Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/LohCG.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief LohCG: Artificial compressibility solver for incompressible flow
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "XystBuildConfig.hpp"
14 : : #include "LohCG.hpp"
15 : : #include "Vector.hpp"
16 : : #include "Reader.hpp"
17 : : #include "ContainerUtil.hpp"
18 : : #include "UnsMesh.hpp"
19 : : #include "ExodusIIMeshWriter.hpp"
20 : : #include "InciterConfig.hpp"
21 : : #include "DerivedData.hpp"
22 : : #include "Discretization.hpp"
23 : : #include "DiagReducer.hpp"
24 : : #include "IntegralReducer.hpp"
25 : : #include "Integrals.hpp"
26 : : #include "Refiner.hpp"
27 : : #include "Reorder.hpp"
28 : : #include "Around.hpp"
29 : : #include "Lohner.hpp"
30 : : #include "Problems.hpp"
31 : : #include "EOS.hpp"
32 : : #include "BC.hpp"
33 : : #include "Print.hpp"
34 : :
35 : : namespace inciter {
36 : :
37 : : extern ctr::Config g_cfg;
38 : :
39 : : static CkReduction::reducerType IntegralsMerger;
40 : :
41 : : //! Runge-Kutta coefficients
42 : : static const std::array< std::vector< tk::real >, 4 > rkcoef{{
43 : : { 1.0 },
44 : : { 1.0/2.0, 1.0 },
45 : : { 1.0/3.0, 1.0/2.0, 1.0 },
46 : : { 1.0/4.0, 1.0/3.0, 1.0/2.0, 1.0 }
47 : : }};
48 : :
49 : : } // inciter::
50 : :
51 : : using inciter::g_cfg;
52 : : using inciter::LohCG;
53 : :
54 : 325 : LohCG::LohCG( const CProxy_Discretization& disc,
55 : : const tk::CProxy_ConjugateGradients& cgpre,
56 : : const std::map< int, std::vector< std::size_t > >& bface,
57 : : const std::map< int, std::vector< std::size_t > >& bnode,
58 : 325 : const std::vector< std::size_t >& triinpoel )
59 : : try :
60 : : m_disc( disc ),
61 [ + - ]: 325 : m_cgpre( cgpre ),
62 : 325 : m_nrhs( 0 ),
63 : 325 : m_nnorm( 0 ),
64 : 325 : m_ngrad( 0 ),
65 : 325 : m_nsgrad( 0 ),
66 : 325 : m_nvgrad( 0 ),
67 : 325 : m_nflux( 0 ),
68 : 325 : m_ndiv( 0 ),
69 : 325 : m_nbpint( 0 ),
70 [ + - ]: 325 : m_np( 0 ),
71 : : m_bnode( bnode ),
72 : : m_bface( bface ),
73 [ + - ][ + - ]: 325 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
74 [ + - ]: 325 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
75 : : m_un( m_u.nunk(), m_u.nprop() ),
76 [ + - ]: 325 : m_grad( m_u.nunk(), ngradcomp() ),
77 : : m_vgrad( m_u.nunk(), 9UL ),
78 : : m_flux( m_u.nunk(), 3UL ),
79 [ + - ][ + - ]: 325 : m_div( m_u.nunk() ),
80 : : m_sgrad( m_u.nunk(), 3UL ),
81 : : m_rhs( m_u.nunk(), m_u.nprop() ),
82 : 325 : m_stage( 0 ),
83 : 325 : m_finished( 0 ),
84 [ + - ][ + - ]: 650 : m_rkcoef( rkcoef[ g_cfg.get< tag::rk >() - 1 ] )
[ + - ]
85 : : // *****************************************************************************
86 : : // Constructor
87 : : //! \param[in] disc Discretization proxy
88 : : //! \param[in] cgpre ConjugateGradients Charm++ proxy for initial pressure solve
89 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
90 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
91 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
92 : : // *****************************************************************************
93 : : {
94 : 325 : usesAtSync = true; // enable migration at AtSync
95 : :
96 [ + - ]: 325 : auto d = Disc();
97 : :
98 : : // Create new local ids based on mesh locality
99 : : std::unordered_map< std::size_t, std::size_t > map;
100 : : std::size_t n = 0;
101 : :
102 [ + - ][ + - ]: 325 : auto psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
103 [ + + ]: 16597 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
104 [ + - ]: 881 : if (!map.count(p)) map[p] = n++;
105 [ + + ][ + + ]: 147156 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
106 [ + - ]: 15391 : if (!map.count(q)) map[q] = n++;
107 : : }
108 : : }
109 : :
110 : : Assert( map.size() == d->Gid().size(),
111 : : "Mesh-locality reorder map size mismatch" );
112 : :
113 : : // Remap data in bound Discretization object
114 [ + - ]: 325 : d->remap( map );
115 : : // Remap boundary triangle face connectivity
116 [ + - ]: 325 : tk::remap( m_triinpoel, map );
117 : : // Recompute points surrounding points
118 [ + - ][ + - ]: 325 : psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
119 : :
120 : : // Compute total box IC volume
121 [ + - ]: 325 : d->boxvol();
122 : :
123 : : // Setup LHS matrix for pressure solve
124 [ + - ][ + - ]: 650 : m_cgpre[ thisIndex ].insert( prelhs( psup ),
[ + - ][ - + ]
[ - - ]
125 : : d->Gid(),
126 : : d->Lid(),
127 : : d->NodeCommMap() );
128 : :
129 : : // Activate SDAG waits for setup
130 [ + - ][ + - ]: 325 : thisProxy[ thisIndex ].wait4int();
131 : :
132 : 325 : } // Catch std::exception
133 [ - - ]: 0 : catch (std::exception& se) {
134 : : // (re-)throw tk::Excpetion
135 [ - - ][ - - ]: 0 : Throw( std::string("RUNTIME ERROR in CSR constructor: ") + se.what() );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
136 : 325 : }
137 : :
138 : :
139 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
140 : 325 : LohCG::prelhs( const std::pair< std::vector< std::size_t >,
141 : : std::vector< std::size_t > >& psup )
142 : : // *****************************************************************************
143 : : // Setup lhs matrix for pressure solve
144 : : //! \param[in] psup Points surrounding points
145 : : //! \return { A, x, b } in linear system A * x = b to solve for pressure
146 : : // *****************************************************************************
147 : : {
148 : 325 : auto d = Disc();
149 : : const auto& inpoel = d->Inpoel();
150 : : const auto& coord = d->Coord();
151 : : const auto& X = coord[0];
152 : : const auto& Y = coord[1];
153 : : const auto& Z = coord[2];
154 : :
155 : : // Matrix with compressed sparse row storage
156 : 325 : tk::CSR A( /*DOF=*/ 1, psup );
157 : :
158 : : // Fill matrix with Laplacian
159 [ + + ]: 34385 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
160 : 34060 : const auto N = inpoel.data() + e*4;
161 : : const std::array< tk::real, 3 >
162 : 34060 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
163 : 34060 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
164 : 34060 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
165 : 34060 : const auto J = tk::triple( ba, ca, da ) * 6.0;
166 : : std::array< std::array< tk::real, 3 >, 4 > grad;
167 : 34060 : grad[1] = tk::cross( ca, da );
168 : 34060 : grad[2] = tk::cross( da, ba );
169 : 34060 : grad[3] = tk::cross( ba, ca );
170 [ + + ]: 136240 : for (std::size_t i=0; i<3; ++i)
171 : 102180 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
172 [ + + ]: 170300 : for (std::size_t a=0; a<4; ++a)
173 [ + + ]: 681200 : for (std::size_t b=0; b<4; ++b)
174 [ + - ]: 544960 : A(N[a],N[b]) -= tk::dot( grad[a], grad[b] ) / J;
175 : : }
176 : :
177 : : auto nunk = X.size();
178 [ + - ][ + - ]: 325 : std::vector< tk::real > x( nunk, 0.0 ), b( nunk, 0.0 );
[ - - ]
179 : :
180 : 325 : return { std::move(A), std::move(x), std::move(b) };
181 : 325 : }
182 : :
183 : : std::size_t
184 : 325 : LohCG::ngradcomp() const
185 : : // *****************************************************************************
186 : : // Compute number of scalar components for gradients
187 : : //! \return Number scalar components required for gradients
188 : : // *****************************************************************************
189 : : {
190 : : std::size_t n = 0;
191 : :
192 [ + + ]: 325 : if (g_cfg.get< tag::flux >() == "damp4") n += 9; // (u,v,w) x 3
193 : :
194 : 325 : return n;
195 : : }
196 : :
197 : : void
198 : 650 : LohCG::setupDirBC( const std::vector< std::vector< int > >& cfgmask,
199 : : const std::vector< std::vector< double > >& cfgval,
200 : : std::size_t ncomp,
201 : : std::vector< std::size_t >& mask,
202 : : std::vector< double >& val )
203 : : // *****************************************************************************
204 : : // Prepare Dirichlet boundary condition data structures
205 : : //! \param[in] cfgmask Boundary condition mask config data to use
206 : : //! \param[in] cfgval Boundary condition values config data to use
207 : : //! \param[in] ncomp Number of scalar component BCs expected per mesh node
208 : : //! \param[in,out] mask Mesh nodes and their Dirichlet BC masks
209 : : //! \param[in,out] val Mesh nodes and their Dirichlet BC values
210 : : // *****************************************************************************
211 : : {
212 : : // Query Dirichlet BC nodes associated to side sets
213 : : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
214 [ + + ]: 733 : for (const auto& s : cfgmask) {
215 : : auto k = m_bface.find(s[0]);
216 [ + + ]: 83 : if (k != end(m_bface)) {
217 [ + - ]: 42 : auto& n = dir[ k->first ];
218 [ + + ]: 12434 : for (auto f : k->second) {
219 [ + - ]: 12392 : n.insert( m_triinpoel[f*3+0] );
220 [ + - ]: 12392 : n.insert( m_triinpoel[f*3+1] );
221 [ + - ]: 12392 : n.insert( m_triinpoel[f*3+2] );
222 : : }
223 : : }
224 : : }
225 : :
226 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
227 [ + - ]: 650 : const auto& lid = Disc()->Lid();
228 [ + + ]: 733 : for (const auto& s : cfgmask) {
229 : : auto k = m_bnode.find(s[0]);
230 [ + + ]: 83 : if (k != end(m_bnode)) {
231 [ + - ]: 42 : auto& n = dir[ k->first ];
232 [ + - ][ + + ]: 8620 : for (auto g : k->second) {
233 : : n.insert( tk::cref_find(lid,g) );
234 : : }
235 : : }
236 : : }
237 : :
238 : : // Associate sidesets to Dirichlet BC values if configured by user
239 : : std::unordered_map< int, std::vector< double > > dirval;
240 [ + + ]: 708 : for (const auto& s : cfgval) {
241 [ + + ]: 58 : auto k = dir.find( static_cast<int>(s[0]) );
242 [ + + ]: 58 : if (k != end(dir)) {
243 [ + - ]: 17 : auto& v = dirval[ k->first ];
244 [ + - ]: 17 : v.resize( s.size()-1 );
245 [ + + ]: 49 : for (std::size_t i=1; i<s.size(); ++i) v[i-1] = s[i];
246 : : }
247 : : }
248 : :
249 : : // Collect unique set of nodes + Dirichlet BC components mask and value
250 : 650 : auto nmask = ncomp + 1;
251 : : std::unordered_map< std::size_t,
252 : : std::pair< std::vector< int >,
253 : : std::vector< double > > > dirbcset;
254 [ + + ]: 733 : for (const auto& vec : cfgmask) {
255 [ - + ][ - - ]: 83 : ErrChk( vec.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
256 : : auto n = dir.find( vec[0] );
257 [ + + ]: 83 : if (n != end(dir)) {
258 [ + - ][ + + ]: 42 : std::vector< double > v( ncomp, 0.0 );
259 : : auto m = dirval.find( vec[0] );
260 [ + + ]: 42 : if (m != end(dirval)) {
261 [ - + ][ - - ]: 17 : ErrChk( m->second.size() == ncomp, "Incorrect Dirichlet BC val ncomp" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
262 [ + - ]: 17 : v = m->second;
263 : : }
264 [ + - ][ + + ]: 7873 : for (auto p : n->second) {
265 : : auto& mv = dirbcset[p]; // mask & value
266 [ + - ][ + - ]: 7831 : if (mv.second.empty()) mv.second = v;
267 [ + - ]: 7831 : auto& mval = mv.first;
268 [ + - ][ + - ]: 7831 : if (mval.empty()) mval.resize( ncomp, 0 );
[ - - ]
269 [ + + ]: 38723 : for (std::size_t c=0; c<ncomp; ++c) {
270 [ + - ]: 30892 : if (!mval[c]) mval[c] = vec[c+1]; // overwrite mask if 0 -> nonzero
271 : : //mval[c] = std::max( mval[c], vec[c+1] );
272 : : //std::cout << p << ": " << mval[c] << '\n';
273 : : }
274 : : }
275 : : }
276 : : }
277 : :
278 : : // Compile streamable list of nodes + Dirichlet BC components mask and values
279 : : tk::destroy( mask );
280 [ + + ]: 8481 : for (const auto& [p,mv] : dirbcset) {
281 [ + - ]: 7831 : mask.push_back( p );
282 [ + - ]: 7831 : mask.insert( end(mask), begin(mv.first), end(mv.first) );
283 [ + - ][ + - ]: 7831 : val.push_back( static_cast< double >( p ) );
284 [ + - ]: 7831 : val.insert( end(val), begin(mv.second), end(mv.second) );
285 : : }
286 : :
287 [ - + ][ - - ]: 650 : ErrChk( mask.size() % nmask == 0, "Dirichlet BC mask incomplete" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
288 [ - + ][ - - ]: 650 : ErrChk( val.size() % nmask == 0, "Dirichlet BC val incomplete" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
289 [ - + ][ - - ]: 650 : ErrChk( mask.size() == val.size(), "Dirichlet BC mask & val size mismatch" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
290 : 650 : }
291 : :
292 : : void
293 : 325 : LohCG::feop()
294 : : // *****************************************************************************
295 : : // Start (re-)computing finite element domain and boundary operators
296 : : // *****************************************************************************
297 : : {
298 : 325 : auto d = Disc();
299 : :
300 : : // Prepare Dirichlet boundary conditions data structures
301 : 325 : setupDirBC( g_cfg.get< tag::bc_dir >(), g_cfg.get< tag::bc_dirval >(),
302 : 325 : m_u.nprop(), m_dirbcmask, m_dirbcval );
303 : 325 : setupDirBC( g_cfg.get< tag::pre_bc_dir >(), g_cfg.get< tag::pre_bc_dirval >(),
304 : 325 : 1, m_dirbcmaskp, m_dirbcvalp );
305 : :
306 : : // Compute local contributions to boundary normals and integrals
307 : 325 : bndint();
308 : : // Compute local contributions to domain edge integrals
309 : 325 : domint();
310 : :
311 : : // Send boundary point normal contributions to neighbor chares
312 [ + + ]: 325 : if (d->NodeCommMap().empty()) {
313 : 7 : comnorm_complete();
314 : : } else {
315 [ + + ]: 3858 : for (const auto& [c,nodes] : d->NodeCommMap()) {
316 : : decltype(m_bnorm) exp;
317 [ + + ]: 17360 : for (auto i : nodes) {
318 [ + + ]: 41640 : for (const auto& [s,b] : m_bnorm) {
319 : : auto k = b.find(i);
320 [ + + ]: 34788 : if (k != end(b)) exp[s][i] = k->second;
321 : : }
322 : : }
323 [ + - ][ + - ]: 7080 : thisProxy[c].comnorm( exp );
324 : : }
325 : : }
326 : 325 : ownnorm_complete();
327 : 325 : }
328 : :
329 : : void
330 : 325 : LohCG::bndint()
331 : : // *****************************************************************************
332 : : // Compute local contributions to boundary normals and integrals
333 : : // *****************************************************************************
334 : : {
335 : 325 : auto d = Disc();
336 : : const auto& coord = d->Coord();
337 : : const auto& gid = d->Gid();
338 : : const auto& x = coord[0];
339 : : const auto& y = coord[1];
340 : : const auto& z = coord[2];
341 : :
342 : : // Lambda to compute the inverse distance squared between boundary face
343 : : // centroid and boundary point. Here p is the global node id and c is the
344 : : // the boundary face centroid.
345 : 70488 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
346 : 70488 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
347 : 70488 : (c[1] - y[p]) * (c[1] - y[p]) +
348 : 70488 : (c[2] - z[p]) * (c[2] - z[p]) );
349 : 325 : };
350 : :
351 : 325 : tk::destroy( m_bnorm );
352 : 325 : tk::destroy( m_bndpoinint );
353 : :
354 [ + + ]: 1060 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
355 [ + + ]: 24231 : for (auto f : faceids) { // for all side set triangles
356 : 23496 : const auto N = m_triinpoel.data() + f*3;
357 : : const std::array< tk::real, 3 >
358 : 23496 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
359 : 23496 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
360 : : auto n = tk::cross( ba, ca );
361 : : auto A2 = tk::length( n );
362 : 23496 : n[0] /= A2;
363 : 23496 : n[1] /= A2;
364 : 23496 : n[2] /= A2;
365 : : const tk::real centroid[3] = {
366 : 23496 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
367 : 23496 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
368 : 23496 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
369 [ + + ]: 93984 : for (const auto& [i,j] : tk::lpoet) {
370 : 70488 : auto p = N[i];
371 : 70488 : tk::real r = invdistsq( centroid, p );
372 : : auto& v = m_bnorm[setid]; // associate side set id
373 : : auto& bpn = v[gid[p]]; // associate global node id of bnd pnt
374 : 70488 : bpn[0] += r * n[0]; // inv.dist.sq-weighted normal
375 : 70488 : bpn[1] += r * n[1];
376 : 70488 : bpn[2] += r * n[2];
377 : 70488 : bpn[3] += r; // inv.dist.sq of node from centroid
378 : : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
379 : 70488 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
380 : 70488 : b[1] += n[1] * A2 / 6.0;
381 : 70488 : b[2] += n[2] * A2 / 6.0;
382 : : }
383 : : }
384 : : }
385 : 325 : }
386 : :
387 : : void
388 : 325 : LohCG::domint()
389 : : // *****************************************************************************
390 : : //! Compute local contributions to domain edge integrals
391 : : // *****************************************************************************
392 : : {
393 : 325 : auto d = Disc();
394 : :
395 : : const auto& gid = d->Gid();
396 : : const auto& inpoel = d->Inpoel();
397 : :
398 : : const auto& coord = d->Coord();
399 : : const auto& x = coord[0];
400 : : const auto& y = coord[1];
401 : : const auto& z = coord[2];
402 : :
403 : 325 : tk::destroy( m_domedgeint );
404 : :
405 [ + + ]: 34385 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
406 : 34060 : const auto N = inpoel.data() + e*4;
407 : : const std::array< tk::real, 3 >
408 : 34060 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
409 : 34060 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
410 : 34060 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
411 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
412 : : Assert( J > 0, "Element Jacobian non-positive" );
413 : : std::array< std::array< tk::real, 3 >, 4 > grad;
414 : 34060 : grad[1] = tk::cross( ca, da );
415 : 34060 : grad[2] = tk::cross( da, ba );
416 : 34060 : grad[3] = tk::cross( ba, ca );
417 [ + + ]: 136240 : for (std::size_t i=0; i<3; ++i)
418 : 102180 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
419 [ + + ]: 238420 : for (const auto& [p,q] : tk::lpoed) {
420 [ + + ]: 204360 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
421 : : tk::real sig = 1.0;
422 [ + + ]: 204360 : if (ed[0] > ed[1]) {
423 : : std::swap( ed[0], ed[1] );
424 : : sig = -1.0;
425 : : }
426 : : auto& n = m_domedgeint[ ed ];
427 : 204360 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
428 : 204360 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
429 : 204360 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
430 : 204360 : n[3] += tk::dot( grad[p], grad[q] ) / J / 6.0;
431 : : }
432 : : }
433 : 325 : }
434 : :
435 : : void
436 : 3540 : LohCG::comnorm( const decltype(m_bnorm)& inbnd )
437 : : // *****************************************************************************
438 : : // Receive contributions to boundary point normals on chare-boundaries
439 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
440 : : // *****************************************************************************
441 : : {
442 : : // Buffer up incoming boundary point normal vector contributions
443 [ + + ]: 6150 : for (const auto& [s,b] : inbnd) {
444 : : auto& bndnorm = m_bnormc[s];
445 [ + + ]: 9578 : for (const auto& [p,n] : b) {
446 : : auto& norm = bndnorm[p];
447 : 6968 : norm[0] += n[0];
448 : 6968 : norm[1] += n[1];
449 : 6968 : norm[2] += n[2];
450 : 6968 : norm[3] += n[3];
451 : : }
452 : : }
453 : :
454 [ + + ]: 3540 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
455 : 318 : m_nnorm = 0;
456 : 318 : comnorm_complete();
457 : : }
458 : 3540 : }
459 : :
460 : : void
461 : 768 : LohCG::registerReducers()
462 : : // *****************************************************************************
463 : : // Configure Charm++ reduction types initiated from this chare array
464 : : //! \details Since this is a [initnode] routine, the runtime system executes the
465 : : //! routine exactly once on every logical node early on in the Charm++ init
466 : : //! sequence. Must be static as it is called without an object. See also:
467 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
468 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
469 : : // *****************************************************************************
470 : : {
471 : 768 : NodeDiagnostics::registerReducers();
472 : 768 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
473 : 768 : }
474 : :
475 : : void
476 : : // cppcheck-suppress unusedFunction
477 : 3065 : LohCG::ResumeFromSync()
478 : : // *****************************************************************************
479 : : // Return from migration
480 : : //! \details This is called when load balancing (LB) completes. The presence of
481 : : //! this function does not affect whether or not we block on LB.
482 : : // *****************************************************************************
483 : : {
484 [ - + ][ - - ]: 3065 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
485 : :
486 [ + - ]: 3065 : if (!g_cfg.get< tag::nonblocking >()) dt();
487 : 3065 : }
488 : :
489 : : void
490 : 325 : LohCG::setup( tk::real v )
491 : : // *****************************************************************************
492 : : // Start setup for solution
493 : : //! \param[in] v Total volume within user-specified box
494 : : // *****************************************************************************
495 : : {
496 : 325 : auto d = Disc();
497 : :
498 : : // Store user-defined box IC volume
499 : 325 : Disc()->Boxvol() = v;
500 : :
501 : : // Set initial conditions
502 : 325 : problems::initialize( d->Coord(), m_u, d->T(), d->BoxNodes() );
503 : :
504 : : // Query time history field output labels from all PDEs integrated
505 [ - + ]: 325 : if (!g_cfg.get< tag::histout >().empty()) {
506 : : std::vector< std::string > var
507 [ - - ][ - - ]: 0 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ][ - - ]
[ - - ]
508 : : auto ncomp = m_u.nprop();
509 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c)
510 [ - - ]: 0 : var.push_back( "c" + std::to_string(c-5) );
511 [ - - ]: 0 : d->histheader( std::move(var) );
512 : 0 : }
513 : :
514 : : // Compute finite element operators
515 : 325 : feop();
516 [ - - ][ - - ]: 325 : }
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
517 : :
518 : : void
519 : 325 : LohCG::bnorm()
520 : : // *****************************************************************************
521 : : // Combine own and communicated portions of the boundary point normals
522 : : // *****************************************************************************
523 : : {
524 : 325 : const auto& lid = Disc()->Lid();
525 : :
526 : : // Combine own and communicated contributions to boundary point normals
527 [ + + ]: 1061 : for (const auto& [s,b] : m_bnormc) {
528 : : auto& bndnorm = m_bnorm[s];
529 [ + + ]: 5675 : for (const auto& [g,n] : b) {
530 : : auto& norm = bndnorm[g];
531 : 4939 : norm[0] += n[0];
532 : 4939 : norm[1] += n[1];
533 : 4939 : norm[2] += n[2];
534 : 4939 : norm[3] += n[3];
535 : : }
536 : : }
537 : 325 : tk::destroy( m_bnormc );
538 : :
539 : : // Divide summed point normals by the sum of the inverse distance squared
540 [ + + ]: 1114 : for (auto& [s,b] : m_bnorm) {
541 [ + + ]: 19394 : for (auto& [g,n] : b) {
542 : 18605 : n[0] /= n[3];
543 : 18605 : n[1] /= n[3];
544 : 18605 : n[2] /= n[3];
545 : : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
546 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
547 : : "Non-unit normal" );
548 : : }
549 : : }
550 : :
551 : : // Replace global->local ids associated to boundary point normals
552 : : decltype(m_bnorm) loc;
553 [ + + ]: 1114 : for (auto& [s,b] : m_bnorm) {
554 : : auto& bnd = loc[s];
555 [ + + ]: 19394 : for (auto&& [g,n] : b) {
556 : 18605 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
557 : : }
558 : : }
559 : : m_bnorm = std::move(loc);
560 : 325 : }
561 : :
562 : : void
563 [ + - ]: 325 : LohCG::streamable()
564 : : // *****************************************************************************
565 : : // Convert integrals into streamable data structures
566 : : // *****************************************************************************
567 : : {
568 : : // Query surface integral output nodes
569 : : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
570 : : const auto& is = g_cfg.get< tag::integout >();
571 [ + - ]: 325 : std::set< int > outsets( begin(is), end(is) );
572 [ - + ]: 325 : for (auto s : outsets) {
573 : : auto m = m_bface.find(s);
574 [ - - ]: 0 : if (m != end(m_bface)) {
575 [ - - ]: 0 : auto& n = surfintnodes[ m->first ]; // associate set id
576 [ - - ]: 0 : for (auto f : m->second) { // face ids on side set
577 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+0] ); // nodes on side set
578 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+1] );
579 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+2] );
580 : : }
581 : : }
582 : : }
583 [ - - ][ - + ]: 325 : for (auto& [s,n] : surfintnodes) tk::unique( n );
584 : : // Prepare surface integral data
585 : 325 : tk::destroy( m_surfint );
586 [ + - ]: 325 : const auto& gid = Disc()->Gid();
587 [ - + ]: 325 : for (auto&& [s,n] : surfintnodes) {
588 [ - - ]: 0 : auto& sint = m_surfint[s]; // associate set id
589 : 0 : auto& nodes = sint.first;
590 : 0 : auto& ndA = sint.second;
591 : : nodes = std::move(n);
592 [ - - ]: 0 : ndA.resize( nodes.size()*3 );
593 : : std::size_t a = 0;
594 [ - - ]: 0 : for (auto p : nodes) {
595 : : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
596 : 0 : ndA[a*3+0] = b[0]; // store ni * dA
597 : 0 : ndA[a*3+1] = b[1];
598 : 0 : ndA[a*3+2] = b[2];
599 : 0 : ++a;
600 : : }
601 : : }
602 : 325 : tk::destroy( m_bndpoinint );
603 : :
604 : : // Generate domain superedges
605 [ + - ]: 325 : domsuped();
606 : :
607 : : // Prepare symmetry boundary condition data structures
608 : :
609 : : // Query symmetry BC nodes associated to side sets
610 : : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
611 [ + + ]: 617 : for (auto s : g_cfg.get< tag::bc_sym >()) {
612 : : auto k = m_bface.find(s);
613 [ + + ]: 292 : if (k != end(m_bface)) {
614 [ + - ]: 110 : auto& n = sym[ k->first ];
615 [ + - ][ + + ]: 1334 : for (auto f : k->second) {
616 [ + - ]: 1224 : const auto& t = m_triinpoel.data() + f*3;
617 : : n.insert( t[0] );
618 [ + - ]: 1224 : n.insert( t[1] );
619 [ + - ]: 1224 : n.insert( t[2] );
620 : : }
621 : : }
622 : : }
623 : :
624 : : // Generate unique set of symmetry BC nodes of all symmetryc BC side sets
625 : : std::set< std::size_t > symbcnodeset;
626 [ + + ]: 435 : for (const auto& [s,n] : sym) symbcnodeset.insert( begin(n), end(n) );
627 : :
628 : : // Generate symmetry BC data as streamable data structures
629 [ - + ]: 325 : tk::destroy( m_symbcnodes );
630 [ - + ]: 325 : tk::destroy( m_symbcnorms );
631 [ + + ]: 1543 : for (auto p : symbcnodeset) {
632 [ + + ]: 2436 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
633 : : auto m = m_bnorm.find( s );
634 [ + - ]: 1218 : if (m != end(m_bnorm)) {
635 : : auto r = m->second.find( p );
636 [ + - ]: 1218 : if (r != end(m->second)) {
637 [ + - ]: 1218 : m_symbcnodes.push_back( p );
638 [ + - ]: 1218 : m_symbcnorms.push_back( r->second[0] );
639 [ + - ]: 1218 : m_symbcnorms.push_back( r->second[1] );
640 [ + - ]: 1218 : m_symbcnorms.push_back( r->second[2] );
641 : : }
642 : : }
643 : : }
644 : : }
645 : 325 : tk::destroy( m_bnorm );
646 : :
647 : : // Prepare noslip boundary condition data structures
648 : :
649 : : // Query noslip BC nodes associated to side sets
650 : : std::unordered_map< int, std::unordered_set< std::size_t > > noslip;
651 [ + + ]: 415 : for (auto s : g_cfg.get< tag::bc_noslip >()) {
652 : : auto k = m_bface.find(s);
653 [ + + ]: 90 : if (k != end(m_bface)) {
654 [ + - ]: 83 : auto& n = noslip[ k->first ];
655 [ + - ][ + + ]: 3843 : for (auto f : k->second) {
656 [ + - ]: 3760 : const auto& t = m_triinpoel.data() + f*3;
657 : : n.insert( t[0] );
658 [ + - ]: 3760 : n.insert( t[1] );
659 [ + - ]: 3760 : n.insert( t[2] );
660 : : }
661 : : }
662 : : }
663 : :
664 : : // Generate unique set of noslip BC nodes of all noslip BC side sets
665 : : std::set< std::size_t > noslipbcnodeset;
666 [ + + ]: 408 : for (const auto& [s,n] : noslip) noslipbcnodeset.insert( begin(n), end(n) );
667 : :
668 : : // Generate noslip BC data as streamable data structures
669 [ - + ]: 325 : tk::destroy( m_noslipbcnodes );
670 [ + - ]: 325 : m_noslipbcnodes.insert( m_noslipbcnodes.end(),
671 : : begin(noslipbcnodeset), end(noslipbcnodeset) );
672 : 325 : }
673 : :
674 : : void
675 : 325 : LohCG::domsuped()
676 : : // *****************************************************************************
677 : : // Generate superedge-groups for domain-edge loops
678 : : //! \see See Lohner, Sec. 15.1.6.2, An Introduction to Applied CFD Techniques,
679 : : //! Wiley, 2008.
680 : : // *****************************************************************************
681 : : {
682 : : Assert( !m_domedgeint.empty(), "No domain edges to group" );
683 : :
684 : : #ifndef NDEBUG
685 : : auto nedge = m_domedgeint.size();
686 : : #endif
687 : :
688 : 325 : const auto& inpoel = Disc()->Inpoel();
689 : 325 : const auto& lid = Disc()->Lid();
690 : 325 : const auto& gid = Disc()->Gid();
691 : :
692 : : tk::destroy( m_dsupedge[0] );
693 : : tk::destroy( m_dsupedge[1] );
694 : : tk::destroy( m_dsupedge[2] );
695 : :
696 : : tk::destroy( m_dsupint[0] );
697 : : tk::destroy( m_dsupint[1] );
698 : : tk::destroy( m_dsupint[2] );
699 : :
700 : : tk::UnsMesh::FaceSet untri;
701 [ + + ]: 34385 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
702 : : std::size_t N[4] = {
703 : 34060 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
704 [ + - ][ + + ]: 170300 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
705 : : }
706 : :
707 [ + + ]: 34385 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
708 : : std::size_t N[4] = {
709 : 34060 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
710 : : int f = 0;
711 : : tk::real sig[6];
712 [ + + ]: 238420 : decltype(m_domedgeint)::const_iterator d[6];
713 [ + + ]: 101235 : for (const auto& [p,q] : tk::lpoed) {
714 [ + + ]: 95556 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
715 [ + + ]: 124495 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
716 : 95556 : d[f] = m_domedgeint.find( ed );
717 [ + + ]: 95556 : if (d[f] == end(m_domedgeint)) break; else ++f;
718 : : }
719 [ + + ]: 34060 : if (f == 6) {
720 [ + - ]: 5679 : m_dsupedge[0].push_back( N[0] );
721 [ + - ]: 5679 : m_dsupedge[0].push_back( N[1] );
722 [ + - ]: 5679 : m_dsupedge[0].push_back( N[2] );
723 [ + - ]: 5679 : m_dsupedge[0].push_back( N[3] );
724 [ + + ]: 28395 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
725 [ + + ]: 39753 : for (int ed=0; ed<6; ++ed) {
726 : : const auto& ded = d[ed]->second;
727 [ + - ]: 34074 : m_dsupint[0].push_back( sig[ed] * ded[0] );
728 [ + - ]: 34074 : m_dsupint[0].push_back( sig[ed] * ded[1] );
729 [ + - ][ + - ]: 34074 : m_dsupint[0].push_back( sig[ed] * ded[2] );
730 [ + - ]: 34074 : m_dsupint[0].push_back( ded[3] );
731 : 34074 : m_domedgeint.erase( d[ed] );
732 : : }
733 : : }
734 : : }
735 : :
736 [ + + ]: 61164 : for (const auto& N : untri) {
737 : : int f = 0;
738 : : tk::real sig[3];
739 [ + + ]: 243356 : decltype(m_domedgeint)::const_iterator d[3];
740 [ + + ]: 99857 : for (const auto& [p,q] : tk::lpoet) {
741 [ + + ]: 94172 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
742 [ + + ]: 144424 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
743 : 94172 : d[f] = m_domedgeint.find( ed );
744 [ + + ]: 94172 : if (d[f] == end(m_domedgeint)) break; else ++f;
745 : : }
746 [ + + ]: 60839 : if (f == 3) {
747 [ + - ]: 5685 : m_dsupedge[1].push_back( N[0] );
748 [ + - ]: 5685 : m_dsupedge[1].push_back( N[1] );
749 [ + - ]: 5685 : m_dsupedge[1].push_back( N[2] );
750 [ + + ]: 22740 : for (int ed=0; ed<3; ++ed) {
751 : : const auto& ded = d[ed]->second;
752 [ + - ]: 17055 : m_dsupint[1].push_back( sig[ed] * ded[0] );
753 [ + - ]: 17055 : m_dsupint[1].push_back( sig[ed] * ded[1] );
754 [ + - ][ + - ]: 17055 : m_dsupint[1].push_back( sig[ed] * ded[2] );
755 [ + - ]: 17055 : m_dsupint[1].push_back( ded[3] );
756 : 17055 : m_domedgeint.erase( d[ed] );
757 : : }
758 : : }
759 : : }
760 : :
761 [ + - ]: 325 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
762 [ + - ]: 325 : m_dsupint[2].resize( m_domedgeint.size()*4 );
763 : : std::size_t k = 0;
764 [ + + ]: 14638 : for (const auto& [ed,d] : m_domedgeint) {
765 : 14313 : auto e = m_dsupedge[2].data() + k*2;
766 : 14313 : e[0] = tk::cref_find( lid, ed[0] );
767 : 14313 : e[1] = tk::cref_find( lid, ed[1] );
768 : 14313 : auto i = m_dsupint[2].data() + k*4;
769 : 14313 : i[0] = d[0];
770 : 14313 : i[1] = d[1];
771 : 14313 : i[2] = d[2];
772 : 14313 : i[3] = d[3];
773 : 14313 : ++k;
774 : : }
775 : :
776 : 325 : tk::destroy( m_domedgeint );
777 : :
778 : : //std::cout << std::setprecision(2)
779 : : // << "superedges: ntet:" << m_dsupedge[0].size()/4 << "(nedge:"
780 : : // << m_dsupedge[0].size()/4*6 << ","
781 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[0].size()/4*6 ) /
782 : : // static_cast< tk::real >( nedge )
783 : : // << "%) + ntri:" << m_dsupedge[1].size()/3
784 : : // << "(nedge:" << m_dsupedge[1].size() << ","
785 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[1].size() ) /
786 : : // static_cast< tk::real >( nedge )
787 : : // << "%) + nedge:"
788 : : // << m_dsupedge[2].size()/2 << "("
789 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[2].size()/2 ) /
790 : : // static_cast< tk::real >( nedge )
791 : : // << "%) = " << m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
792 : : // m_dsupedge[2].size()/2 << " of "<< nedge << " total edges\n";
793 : :
794 : : Assert( m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
795 : : m_dsupedge[2].size()/2 == nedge,
796 : : "Not all edges accounted for in superedge groups" );
797 : 325 : }
798 : :
799 : : void
800 : : // cppcheck-suppress unusedFunction
801 : 325 : LohCG::merge()
802 : : // *****************************************************************************
803 : : // Combine own and communicated portions of the integrals
804 : : // *****************************************************************************
805 : : {
806 : 325 : auto d = Disc();
807 : :
808 : : // Combine own and communicated contributions to boundary point normals
809 : 325 : bnorm();
810 : :
811 : : // Convert integrals into streamable data structures
812 : 325 : streamable();
813 : :
814 : : // Enforce boundary conditions on initial conditions
815 : 325 : auto t = d->T() + d->Dt();
816 : 325 : physics::dirbc( m_u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
817 : 325 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
818 : 325 : physics::noslipbc( m_u, m_noslipbcnodes, /*pos=*/1 );
819 : :
820 : : // Start measuring initial div-free time
821 : 325 : m_timer.emplace_back();
822 : :
823 : : // Compute initial momentum flux
824 [ + - ]: 325 : thisProxy[ thisIndex ].wait4div();
825 [ + - ]: 325 : thisProxy[ thisIndex ].wait4sgrad();
826 : 325 : div( m_u, /*pos=*/1 );
827 : 325 : }
828 : :
829 : : void
830 : 975 : LohCG::fingrad( tk::Fields& grad,
831 : : std::unordered_map< std::size_t, std::vector< tk::real > >& gradc )
832 : : // *****************************************************************************
833 : : // Finalize computing gradient
834 : : //! \param[in,out] grad Gradient to finalize
835 : : //! \param[in,out] gradc Gradient communication buffer to finalize
836 : : // *****************************************************************************
837 : : {
838 : 975 : auto d = Disc();
839 : : const auto lid = d->Lid();
840 : :
841 : : // Combine own and communicated contributions
842 [ + + ]: 18942 : for (const auto& [g,r] : gradc) {
843 : 17967 : auto i = tk::cref_find( lid, g );
844 [ + + ]: 107802 : for (std::size_t c=0; c<r.size(); ++c) grad(i,c) += r[c];
845 : : }
846 : 975 : tk::destroy(gradc);
847 : :
848 : : // Divide weak result by nodal volume
849 : : const auto& vol = d->Vol();
850 [ + + ]: 49791 : for (std::size_t p=0; p<grad.nunk(); ++p) {
851 [ + + ]: 292896 : for (std::size_t c=0; c<grad.nprop(); ++c) {
852 : 244080 : grad(p,c) /= vol[p];
853 : : }
854 : : }
855 : 975 : }
856 : :
857 : : void
858 : 650 : LohCG::div( const tk::Fields& u, std::size_t pos )
859 : : // *****************************************************************************
860 : : // Start computing divergence
861 : : //! \param[in] u Vector field whose divergence to compute
862 : : //! \param[in] pos Position at which the three vector components are in u
863 : : // *****************************************************************************
864 : : {
865 : 650 : auto d = Disc();
866 : : const auto lid = d->Lid();
867 : :
868 : : // Finalize momentum flux communications if needed
869 [ + + ]: 650 : if (m_np == 1) {
870 [ + - ]: 325 : fingrad( m_flux, m_fluxc );
871 [ + - ]: 325 : physics::symbc( m_flux, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
872 : : }
873 : :
874 : : // Compute divergence
875 : : std::fill( begin(m_div), end(m_div), 0.0 );
876 [ + - ]: 650 : lohner::div( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, u, m_div, pos );
877 : :
878 : : // Communicate velocity divergence to other chares on chare-boundary
879 [ + + ]: 650 : if (d->NodeCommMap().empty()) {
880 [ + - ]: 14 : comdiv_complete();
881 : : } else {
882 [ + + ]: 7716 : for (const auto& [c,n] : d->NodeCommMap()) {
883 : : decltype(m_divc) exp;
884 [ + - ][ + + ]: 34720 : for (auto g : n) exp[g] = m_div[ tk::cref_find(lid,g) ];
885 [ + - ][ + - ]: 14160 : thisProxy[c].comdiv( exp );
886 : : }
887 : : }
888 [ + - ]: 650 : owndiv_complete();
889 : 650 : }
890 : :
891 : : void
892 : 7080 : LohCG::comdiv( const std::unordered_map< std::size_t, tk::real >& indiv )
893 : : // *****************************************************************************
894 : : // Receive contributions to velocity divergence on chare-boundaries
895 : : //! \param[in] indiv Partial contributions of velocity divergence to
896 : : //! chare-boundary nodes. Key: global mesh node IDs, value: contribution.
897 : : //! \details This function receives contributions to m_div, which stores the
898 : : //! velocity divergence at mesh nodes. While m_div stores own contributions,
899 : : //! m_divc collects the neighbor chare contributions during communication.
900 : : //! This way work on m_div and m_divc is overlapped.
901 : : // *****************************************************************************
902 : : {
903 [ + + ]: 34720 : for (const auto& [g,r] : indiv) m_divc[g] += r;
904 : :
905 : : // When we have heard from all chares we communicate with, this chare is done
906 [ + + ]: 7080 : if (++m_ndiv == Disc()->NodeCommMap().size()) {
907 : 636 : m_ndiv = 0;
908 : 636 : comdiv_complete();
909 : : }
910 : 7080 : }
911 : :
912 : : void
913 : 325 : LohCG::velgrad()
914 : : // *****************************************************************************
915 : : // Start computing velocity gradient
916 : : // *****************************************************************************
917 : : {
918 : 325 : auto d = Disc();
919 : :
920 : : // Compute momentum flux
921 : : m_vgrad.fill( 0.0 );
922 : 325 : lohner::vgrad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad );
923 : :
924 : : // Communicate velocity divergence to other chares on chare-boundary
925 : : const auto& lid = d->Lid();
926 [ + + ]: 325 : if (d->NodeCommMap().empty()) {
927 : 7 : comvgrad_complete();
928 : : } else {
929 [ + + ]: 3858 : for (const auto& [c,n] : d->NodeCommMap()) {
930 : : decltype(m_vgradc) exp;
931 [ + - ][ + + ]: 17360 : for (auto g : n) exp[g] = m_vgrad[ tk::cref_find(lid,g) ];
932 [ + - ][ + - ]: 7080 : thisProxy[c].comvgrad( exp );
933 : : }
934 : : }
935 : 325 : ownvgrad_complete();
936 : 325 : }
937 : :
938 : : void
939 : 3540 : LohCG::comvgrad(
940 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
941 : : // *****************************************************************************
942 : : // Receive contributions to velocity gradients on chare-boundaries
943 : : //! \param[in] ingrad Partial contributions of momentum flux to
944 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
945 : : //! \details This function receives contributions to m_vgrad, which stores the
946 : : //! velocity gradients at mesh nodes. While m_vgrad stores own contributions,
947 : : //! m_vgradc collects the neighbor chare contributions during communication.
948 : : //! This way work on m_vgrad and m_vgradc is overlapped.
949 : : // *****************************************************************************
950 : : {
951 : : using tk::operator+=;
952 [ + + ]: 31180 : for (const auto& [g,r] : ingrad) m_vgradc[g] += r;
953 : :
954 : : // When we have heard from all chares we communicate with, this chare is done
955 [ + + ]: 3540 : if (++m_nvgrad == Disc()->NodeCommMap().size()) {
956 : 318 : m_nvgrad = 0;
957 : 318 : comvgrad_complete();
958 : : }
959 : 3540 : }
960 : :
961 : : void
962 : 325 : LohCG::flux()
963 : : // *****************************************************************************
964 : : // Start computing momentum flux
965 : : // *****************************************************************************
966 : : {
967 : 325 : auto d = Disc();
968 : :
969 : : // Finalize computing velocity gradients
970 : 325 : fingrad( m_vgrad, m_vgradc );
971 : :
972 : : // Compute momentum flux
973 : : m_flux.fill( 0.0 );
974 : 325 : lohner::flux( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad,
975 : 325 : m_flux );
976 : :
977 : : // Communicate velocity divergence to other chares on chare-boundary
978 : : const auto& lid = d->Lid();
979 [ + + ]: 325 : if (d->NodeCommMap().empty()) {
980 : 7 : comflux_complete();
981 : : } else {
982 [ + + ]: 3858 : for (const auto& [c,n] : d->NodeCommMap()) {
983 : : decltype(m_fluxc) exp;
984 [ + - ][ + + ]: 17360 : for (auto g : n) exp[g] = m_flux[ tk::cref_find(lid,g) ];
985 [ + - ][ + - ]: 7080 : thisProxy[c].comflux( exp );
986 : : }
987 : : }
988 : 325 : ownflux_complete();
989 : 325 : }
990 : :
991 : : void
992 : 3540 : LohCG::comflux(
993 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& influx )
994 : : // *****************************************************************************
995 : : // Receive contributions to momentum flux on chare-boundaries
996 : : //! \param[in] influx Partial contributions of momentum flux to
997 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
998 : : //! \details This function receives contributions to m_flux, which stores the
999 : : //! momentum flux at mesh nodes. While m_flux stores own contributions,
1000 : : //! m_fluxc collects the neighbor chare contributions during communication.
1001 : : //! This way work on m_flux and m_fluxc is overlapped.
1002 : : // *****************************************************************************
1003 : : {
1004 : : using tk::operator+=;
1005 [ + + ]: 31180 : for (const auto& [g,r] : influx) m_fluxc[g] += r;
1006 : :
1007 : : // When we have heard from all chares we communicate with, this chare is done
1008 [ + + ]: 3540 : if (++m_nflux == Disc()->NodeCommMap().size()) {
1009 : 318 : m_nflux = 0;
1010 : 318 : comflux_complete();
1011 : : }
1012 : 3540 : }
1013 : :
1014 : : void
1015 : 650 : LohCG::pinit()
1016 : : // *****************************************************************************
1017 : : // Initialize Poisson solve
1018 : : // *****************************************************************************
1019 : : {
1020 : 650 : auto d = Disc();
1021 : : const auto lid = d->Lid();
1022 : : const auto& coord = d->Coord();
1023 : : const auto& x = coord[0];
1024 : : const auto& y = coord[1];
1025 : : const auto& z = coord[2];
1026 : :
1027 : : // Combine own and communicated contributions to velocity divergence
1028 [ + + ]: 12628 : for (const auto& [g,r] : m_divc) m_div[ tk::cref_find(lid,g) ] += r;
1029 : 650 : tk::destroy(m_divc);
1030 : :
1031 : : // Configure Dirichlet BCs
1032 : : std::unordered_map< std::size_t,
1033 : : std::vector< std::pair< int, tk::real > > > dirbc;
1034 [ + + ]: 650 : if (!g_cfg.get< tag::pre_bc_dir >().empty()) {
1035 [ + - ]: 50 : auto ic = problems::PRESSURE_IC();
1036 : : std::size_t nmask = 1 + 1;
1037 : : Assert( m_dirbcmaskp.size() % nmask == 0, "Size mismatch" );
1038 [ + + ]: 338 : for (std::size_t i=0; i<m_dirbcmaskp.size()/nmask; ++i) {
1039 [ - + ]: 288 : auto p = m_dirbcmaskp[i*nmask+0]; // local node id
1040 : 288 : auto mask = m_dirbcmaskp[i*nmask+1];
1041 [ - + ]: 288 : if (mask == 1) { // mask == 1: IC value
1042 [ - - ][ - - ]: 0 : auto val = ic( x[p], y[p], z[p] );
1043 [ - - ]: 0 : dirbc[p] = {{ { 1, val } }};
1044 [ + - ][ + - ]: 288 : } else if (mask == 2 && !m_dirbcvalp.empty()) { // mask == 2: BC value
1045 [ + - ]: 288 : auto val = m_dirbcvalp[i*nmask+1];
1046 [ + - ][ - - ]: 576 : dirbc[p] = {{ { 1, val } }};
1047 : : }
1048 : : }
1049 : : }
1050 : :
1051 : : // Configure Neumann BCs
1052 : : std::vector< tk::real > neubc;
1053 [ + - ]: 650 : auto pg = problems::PRESSURE_GRAD();
1054 [ - + ]: 650 : if (pg) {
1055 : : // Collect Neumann BC elements
1056 [ - - ][ - - ]: 0 : std::vector< std::uint8_t > besym( m_triinpoel.size(), 0 );
1057 [ - - ]: 0 : for (auto s : g_cfg.get< tag::pre_bc_sym >()) {
1058 : : auto k = m_bface.find(s);
1059 [ - - ][ - - ]: 0 : if (k != end(m_bface)) for (auto f : k->second) besym[f] = 1;
1060 : : }
1061 : : // Setup Neumann BCs
1062 [ - - ][ - - ]: 0 : neubc.resize( x.size(), 0.0 );
1063 [ - - ]: 0 : for (std::size_t e=0; e<m_triinpoel.size()/3; ++e) {
1064 [ - - ]: 0 : if (besym[e]) {
1065 : 0 : const auto N = m_triinpoel.data() + e*3;
1066 : : tk::real n[3];
1067 [ - - ]: 0 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
1068 [ - - ]: 0 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
1069 : : n[0], n[1], n[2] );
1070 : 0 : auto g = pg( x[N[0]], y[N[0]], z[N[0]] );
1071 [ - - ]: 0 : neubc[ N[0] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1072 [ - - ]: 0 : g = pg( x[N[1]], y[N[1]], z[N[1]] );
1073 [ - - ]: 0 : neubc[ N[1] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1074 [ - - ]: 0 : g = pg( x[N[2]], y[N[2]], z[N[2]] );
1075 : 0 : neubc[ N[2] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1076 : : }
1077 : : }
1078 : : }
1079 : :
1080 : : // Set hydrostat
1081 : 650 : auto h = g_cfg.get< tag::pre_hydrostat >();
1082 [ + + ]: 650 : if (h != std::numeric_limits< uint64_t >::max()) {
1083 [ + + ]: 16 : auto g = lid.find( 1 );
1084 [ + + ][ + - ]: 24 : if (g != end(lid)) dirbc[ h ] = {{ { g->second, 0.0 }} };
1085 : : }
1086 : :
1087 : : // Configure right hand side
1088 [ + - ]: 650 : auto pr = problems::PRESSURE_RHS();
1089 [ - + ]: 650 : if (pr) {
1090 : : const auto& vol = d->Vol();
1091 [ - - ]: 0 : for (std::size_t i=0; i<x.size(); ++i) {
1092 [ - - ]: 0 : m_div[i] = pr( x[i], y[i], z[i] ) * vol[i];
1093 : : }
1094 : : }
1095 : :
1096 : : // Initialize Poisson solve
1097 : : const auto& pc = g_cfg.get< tag::pre_pc >();
1098 [ + - ][ + - ]: 2600 : m_cgpre[ thisIndex ].ckLocal()->init( {}, m_div, neubc, dirbc, pc,
[ - + ][ - - ]
1099 [ + - ][ + - ]: 1950 : CkCallback( CkIndex_LohCG::psolve(), thisProxy[thisIndex] ) );
[ - + ][ - - ]
1100 : 650 : }
1101 : :
1102 : : void
1103 : 650 : LohCG::psolve()
1104 : : // *****************************************************************************
1105 : : // Solve Poisson equation
1106 : : // *****************************************************************************
1107 : : {
1108 : 650 : auto iter = g_cfg.get< tag::pre_iter >();
1109 : 650 : auto tol = g_cfg.get< tag::pre_tol >();
1110 : 650 : auto verbose = g_cfg.get< tag::pre_verbose >();
1111 : :
1112 [ + + ]: 650 : auto c = m_np != 1 ?
1113 : 325 : CkCallback( CkIndex_LohCG::sgrad(), thisProxy[thisIndex] ) :
1114 [ + - ][ + - ]: 1625 : CkCallback( CkIndex_LohCG::psolved(), thisProxy[thisIndex] );
[ + - ][ + - ]
[ + + ][ - - ]
1115 : :
1116 [ + - ][ + - ]: 1950 : m_cgpre[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose, c );
[ - + ][ - - ]
1117 : 650 : }
1118 : :
1119 : : void
1120 : 325 : LohCG::sgrad()
1121 : : // *****************************************************************************
1122 : : // Compute recent conjugate gradients solution gradient
1123 : : // *****************************************************************************
1124 : : {
1125 : 325 : auto d = Disc();
1126 : :
1127 [ + - ]: 975 : auto sol = m_cgpre[ thisIndex ].ckLocal()->solution();
1128 : : m_sgrad.fill( 0.0 );
1129 [ + - ]: 325 : lohner::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, sol, m_sgrad );
1130 : :
1131 : : // Send gradient contributions to neighbor chares
1132 [ + + ]: 325 : if (d->NodeCommMap().empty()) {
1133 [ + - ]: 7 : comsgrad_complete();
1134 : : } else {
1135 : : const auto& lid = d->Lid();
1136 [ + + ]: 3858 : for (const auto& [c,n] : d->NodeCommMap()) {
1137 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1138 [ + - ][ + + ]: 17360 : for (auto g : n) exp[g] = m_sgrad[ tk::cref_find(lid,g) ];
1139 [ + - ][ + - ]: 7080 : thisProxy[c].comsgrad( exp );
1140 : : }
1141 : : }
1142 [ + - ]: 325 : ownsgrad_complete();
1143 : 325 : }
1144 : :
1145 : : void
1146 : 3540 : LohCG::comsgrad(
1147 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
1148 : : // *****************************************************************************
1149 : : // Receive contributions to conjugrate gradients solution gradient
1150 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
1151 : : //! global mesh node IDs, value: contributions for all scalar components.
1152 : : //! \details This function receives contributions to m_sgrad, which stores the
1153 : : //! gradients at mesh nodes. While m_sgrad stores own contributions, m_sgradc
1154 : : //! collects the neighbor chare contributions during communication. This way
1155 : : //! work on m_sgrad and m_sgradc is overlapped.
1156 : : // *****************************************************************************
1157 : : {
1158 : : using tk::operator+=;
1159 [ + + ]: 31180 : for (const auto& [g,r] : ingrad) m_sgradc[g] += r;
1160 : :
1161 : : // When we have heard from all chares we communicate with, this chare is done
1162 [ + + ]: 3540 : if (++m_nsgrad == Disc()->NodeCommMap().size()) {
1163 : 318 : m_nsgrad = 0;
1164 : 318 : comsgrad_complete();
1165 : : }
1166 : 3540 : }
1167 : :
1168 : : void
1169 : 650 : LohCG::psolved()
1170 : : // *****************************************************************************
1171 : : // Continue setup after Poisson solve and gradient computation
1172 : : // *****************************************************************************
1173 : : {
1174 : 650 : auto d = Disc();
1175 : :
1176 [ + + ][ + - ]: 762 : if (thisIndex == 0) d->pit( m_cgpre[ thisIndex ].ckLocal()->it() );
1177 : :
1178 [ + + ]: 650 : if (m_np != 1) {
1179 : : // Finalize gradient communications
1180 : 325 : fingrad( m_sgrad, m_sgradc );
1181 : : // Project velocity to divergence-free subspace
1182 [ + + ]: 16597 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1183 : 16272 : m_u(i,1) -= m_sgrad(i,0);
1184 : 16272 : m_u(i,2) -= m_sgrad(i,1);
1185 : 16272 : m_u(i,3) -= m_sgrad(i,2);
1186 : : }
1187 : : // Enforce boundary conditions
1188 : 325 : auto t = d->T() + d->Dt();
1189 : 325 : physics::dirbc( m_u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
1190 : 325 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
1191 : 325 : physics::noslipbc( m_u, m_noslipbcnodes, /*pos=*/1 );
1192 : : }
1193 : :
1194 [ - + ]: 650 : if (g_cfg.get< tag::nstep >() == 1) { // test first Poisson solve only
1195 : :
1196 [ - - ]: 0 : auto p = m_cgpre[ thisIndex ].ckLocal()->solution();
1197 [ - - ]: 0 : for (std::size_t i=0; i<m_u.nunk(); ++i) m_u(i,0) = p[i];
1198 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4step();
1199 [ - - ][ - - ]: 0 : writeFields( CkCallback(CkIndex_LohCG::diag(), thisProxy[thisIndex]) );
[ - - ][ - - ]
[ - - ]
1200 : :
1201 : : } else {
1202 : :
1203 [ + + ]: 650 : if (++m_np < 2) {
1204 : : // Compute momentum flux for initial pressure-Poisson rhs
1205 [ + - ]: 325 : thisProxy[ thisIndex ].wait4vgrad();
1206 [ + - ]: 325 : thisProxy[ thisIndex ].wait4flux();
1207 [ + - ]: 325 : thisProxy[ thisIndex ].wait4div();
1208 : 325 : velgrad();
1209 : : } else {
1210 [ + + ]: 325 : if (thisIndex == 0) {
1211 : : tk::Print() << "Initial div-free time: " << m_timer[0].dsec()
1212 : : << " sec\n";
1213 : : }
1214 : : // Assign initial pressure and start timestepping
1215 [ + - ]: 650 : auto p = m_cgpre[ thisIndex ].ckLocal()->solution();
1216 [ + + ]: 16597 : for (std::size_t i=0; i<m_u.nunk(); ++i) m_u(i,0) = p[i];
1217 [ + - ][ + - ]: 975 : writeFields( CkCallback(CkIndex_LohCG::start(), thisProxy[thisIndex]) );
[ + - ][ + - ]
[ - - ]
1218 : : }
1219 : :
1220 : : }
1221 : 650 : }
1222 : :
1223 : : void
1224 : 3530 : LohCG::diag()
1225 : : // *****************************************************************************
1226 : : // Compute diagnostics
1227 : : // *****************************************************************************
1228 : : {
1229 : 3530 : auto d = Disc();
1230 : :
1231 : : // Increase number of iterations and physical time
1232 : 3530 : d->next();
1233 : :
1234 : : // Compute diagnostics, e.g., residuals
1235 : 3530 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1236 : 3530 : auto diagnostics = m_diag.accompute( *d, m_u, m_un, diag_iter );
1237 : :
1238 : : // Evaluate residuals
1239 [ - + ][ - - ]: 3530 : if (!diagnostics) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
1240 : 3530 : }
1241 : :
1242 : : void
1243 : 325 : LohCG::start()
1244 : : // *****************************************************************************
1245 : : // Start time stepping
1246 : : // *****************************************************************************
1247 : : {
1248 : : // Set flag that indicates that we are now during time stepping
1249 : 325 : Disc()->Initial( 0 );
1250 : : // Start timer measuring time stepping wall clock time
1251 : 325 : Disc()->Timer().zero();
1252 : : // Zero grind-timer
1253 : 325 : Disc()->grindZero();
1254 : : // Continue to first time step
1255 : 325 : dt();
1256 : 325 : }
1257 : :
1258 : : void
1259 : 3530 : LohCG::dt()
1260 : : // *****************************************************************************
1261 : : // Compute time step size
1262 : : // *****************************************************************************
1263 : : {
1264 : 3530 : auto d = Disc();
1265 : : const auto& vol = d->Vol();
1266 : :
1267 : 3530 : tk::real mindt = std::numeric_limits< tk::real >::max();
1268 [ + + ]: 3530 : auto const_dt = g_cfg.get< tag::dt >();
1269 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1270 : :
1271 : : // use constant dt if configured
1272 [ + + ]: 3530 : if (std::abs(const_dt) > eps) {
1273 : :
1274 : : // cppcheck-suppress redundantInitialization
1275 : 2920 : mindt = const_dt;
1276 : :
1277 : : } else {
1278 : :
1279 : 610 : auto cfl = g_cfg.get< tag::cfl >();
1280 : 610 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1281 : : auto large = std::numeric_limits< tk::real >::max();
1282 : 610 : auto c = g_cfg.get< tag::soundspeed >();
1283 : :
1284 [ + + ]: 166520 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1285 : 165910 : auto u = m_u(i,1);
1286 : 165910 : auto v = m_u(i,2);
1287 : 165910 : auto w = m_u(i,3);
1288 : 165910 : auto vel = std::sqrt( u*u + v*v + w*w );
1289 : 165910 : auto L = std::cbrt( vol[i] );
1290 [ - + ][ + + ]: 165910 : auto euler_dt = L / std::max( vel+c, 1.0e-8 );
1291 : 165910 : mindt = std::min( mindt, euler_dt );
1292 [ + - ][ - + ]: 165910 : auto visc_dt = mu > eps ? L * L / mu : large;
1293 : 165910 : mindt = std::min( mindt, visc_dt );
1294 : : }
1295 : 610 : mindt *= cfl;
1296 : :
1297 : : }
1298 : :
1299 : : // Actiavate SDAG waits for next time step stage
1300 [ + - ]: 3530 : thisProxy[ thisIndex ].wait4step();
1301 : :
1302 : : // Contribute to minimum dt across all chares and advance to next step
1303 [ + - ]: 3530 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
1304 : 3530 : CkCallback(CkReductionTarget(LohCG,advance), thisProxy) );
1305 : 3530 : }
1306 : :
1307 : : void
1308 : 3530 : LohCG::advance( tk::real newdt )
1309 : : // *****************************************************************************
1310 : : // Advance equations to next time step
1311 : : //! \param[in] newdt The smallest dt across the whole problem
1312 : : // *****************************************************************************
1313 : : {
1314 : : // Set new time step size
1315 : 3530 : Disc()->setdt( newdt );
1316 : :
1317 : : // Start next time step stage
1318 : 3530 : stage();
1319 : 3530 : }
1320 : :
1321 : : void
1322 : 4660 : LohCG::stage()
1323 : : // *****************************************************************************
1324 : : // Start next time step stage
1325 : : // *****************************************************************************
1326 : : {
1327 : : // Activate SDAG waits for next time step stage
1328 [ + - ]: 4660 : thisProxy[ thisIndex ].wait4grad();
1329 [ + - ]: 4660 : thisProxy[ thisIndex ].wait4rhs();
1330 : :
1331 : : // Compute gradients
1332 : 4660 : grad();
1333 : 4660 : }
1334 : :
1335 : : void
1336 : 4660 : LohCG::grad()
1337 : : // *****************************************************************************
1338 : : // Compute gradients
1339 : : // *****************************************************************************
1340 : : {
1341 : 4660 : auto d = Disc();
1342 : :
1343 [ + + ]: 4660 : if (m_grad.nprop()) {
1344 : : m_grad.fill( 0.0 );
1345 : 140 : lohner::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_grad );
1346 : : }
1347 : :
1348 : : // Send gradient contributions to neighbor chares
1349 [ + + ]: 4660 : if (d->NodeCommMap().empty()) {
1350 : 280 : comgrad_complete();
1351 : : } else {
1352 : : const auto& lid = d->Lid();
1353 [ + + ]: 41860 : for (const auto& [c,n] : d->NodeCommMap()) {
1354 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1355 [ + - ][ + + ]: 202680 : for (auto g : n) exp[g] = m_grad[ tk::cref_find(lid,g) ];
1356 [ + - ][ + - ]: 74960 : thisProxy[c].comgrad( exp );
1357 : : }
1358 : : }
1359 : 4660 : owngrad_complete();
1360 : 4660 : }
1361 : :
1362 : : void
1363 : 37480 : LohCG::comgrad(
1364 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
1365 : : // *****************************************************************************
1366 : : // Receive contributions to gradients on chare-boundaries
1367 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
1368 : : //! global mesh node IDs, value: contributions for all scalar components.
1369 : : //! \details This function receives contributions to m_grad, which stores the
1370 : : //! gradients at mesh nodes. While m_grad stores own contributions, m_gradc
1371 : : //! collects the neighbor chare contributions during communication. This way
1372 : : //! work on m_grad and m_gradc is overlapped.
1373 : : // *****************************************************************************
1374 : : {
1375 : : using tk::operator+=;
1376 [ + + ]: 367880 : for (const auto& [g,r] : ingrad) m_gradc[g] += r;
1377 : :
1378 : : // When we have heard from all chares we communicate with, this chare is done
1379 [ + + ]: 37480 : if (++m_ngrad == Disc()->NodeCommMap().size()) {
1380 : 4380 : m_ngrad = 0;
1381 : 4380 : comgrad_complete();
1382 : : }
1383 : 37480 : }
1384 : :
1385 : : void
1386 : 4660 : LohCG::rhs()
1387 : : // *****************************************************************************
1388 : : // Compute right-hand side of transport equations
1389 : : // *****************************************************************************
1390 : : {
1391 : 4660 : auto d = Disc();
1392 : : const auto& lid = d->Lid();
1393 : :
1394 [ + + ]: 4660 : if (m_grad.nprop()) {
1395 : : // Combine own and communicated contributions to gradients
1396 [ - + ]: 140 : for (const auto& [g,r] : m_gradc) {
1397 : 0 : auto i = tk::cref_find( lid, g );
1398 [ - - ]: 0 : for (std::size_t c=0; c<r.size(); ++c) m_grad(i,c) += r[c];
1399 : : }
1400 : 140 : tk::destroy(m_gradc);
1401 : :
1402 : : // Divide weak result by nodal volume
1403 : : const auto& vol = d->Vol();
1404 [ + + ]: 169820 : for (std::size_t p=0; p<m_grad.nunk(); ++p) {
1405 [ + + ]: 1696800 : for (std::size_t c=0; c<m_grad.nprop(); ++c) {
1406 : 1527120 : m_grad(p,c) /= vol[p];
1407 : : }
1408 : : }
1409 : : }
1410 : :
1411 : : // Compute own portion of right-hand side for all equations
1412 : 4660 : lohner::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_grad,
1413 : 4660 : m_rhs );
1414 : :
1415 : : // Communicate rhs to other chares on chare-boundary
1416 [ + + ]: 4660 : if (d->NodeCommMap().empty()) {
1417 : 280 : comrhs_complete();
1418 : : } else {
1419 [ + + ]: 41860 : for (const auto& [c,n] : d->NodeCommMap()) {
1420 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1421 [ + - ][ + + ]: 202680 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
1422 [ + - ][ + - ]: 74960 : thisProxy[c].comrhs( exp );
1423 : : }
1424 : : }
1425 : 4660 : ownrhs_complete();
1426 : 4660 : }
1427 : :
1428 : : void
1429 : 37480 : LohCG::comrhs(
1430 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1431 : : // *****************************************************************************
1432 : : // Receive contributions to right-hand side vector on chare-boundaries
1433 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1434 : : //! global mesh node IDs, value: contributions for all scalar components.
1435 : : //! \details This function receives contributions to m_rhs, which stores the
1436 : : //! right hand side vector at mesh nodes. While m_rhs stores own
1437 : : //! contributions, m_rhsc collects the neighbor chare contributions during
1438 : : //! communication. This way work on m_rhs and m_rhsc is overlapped.
1439 : : // *****************************************************************************
1440 : : {
1441 : : using tk::operator+=;
1442 [ + + ]: 367880 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
1443 : :
1444 : : // When we have heard from all chares we communicate with, this chare is done
1445 [ + + ]: 37480 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
1446 : 4380 : m_nrhs = 0;
1447 : 4380 : comrhs_complete();
1448 : : }
1449 : 37480 : }
1450 : :
1451 : : void
1452 : : // cppcheck-suppress unusedFunction
1453 : 4660 : LohCG::solve()
1454 : : // *****************************************************************************
1455 : : // Advance systems of equations
1456 : : // *****************************************************************************
1457 : : {
1458 : 4660 : auto d = Disc();
1459 : : const auto npoin = m_u.nunk();
1460 : : const auto ncomp = m_u.nprop();
1461 : : const auto& vol = d->Vol();
1462 : :
1463 : : // Combine own and communicated contributions to rhs
1464 : : const auto lid = d->Lid();
1465 [ + + ]: 91160 : for (const auto& [g,r] : m_rhsc) {
1466 : 86500 : auto i = tk::cref_find( lid, g );
1467 [ + + ]: 432500 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1468 : : }
1469 : 4660 : tk::destroy(m_rhsc);
1470 : :
1471 [ + + ]: 4660 : if (m_stage == 0) m_un = m_u;
1472 : :
1473 : : // Advance system
1474 : 4660 : auto dt = m_rkcoef[m_stage] * d->Dt();
1475 [ + + ]: 587700 : for (std::size_t i=0; i<npoin; ++i) {
1476 [ + + ]: 2915200 : for (std::size_t c=0; c<ncomp; ++c) {
1477 : 2332160 : m_u(i,c) = m_un(i,c) - dt*m_rhs(i,c)/vol[i];
1478 : : }
1479 : : }
1480 : :
1481 : : // Configure and apply scalar source to solution (if defined)
1482 [ + - ]: 4660 : auto src = problems::PHYS_SRC();
1483 [ - + ][ - - ]: 4660 : if (src) src( d->Coord(), d->T(), m_u );
1484 : :
1485 : : // Enforce boundary conditions
1486 [ + - ]: 4660 : auto t = d->T() + m_rkcoef[m_stage] * d->Dt();
1487 [ + - ]: 4660 : physics::dirbc( m_u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
1488 [ + - ]: 4660 : physics::dirbcp( m_u, d->Coord(), m_dirbcmaskp, m_dirbcvalp );
1489 [ + - ]: 4660 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
1490 [ + - ]: 4660 : physics::noslipbc( m_u, m_noslipbcnodes, /*pos=*/1 );
1491 : :
1492 [ + + ]: 4660 : if (++m_stage < m_rkcoef.size()) {
1493 : :
1494 : : // Start next time step stage
1495 [ + - ]: 1130 : stage();
1496 : :
1497 : : } else {
1498 : :
1499 : : // Reset Runge-Kutta stage counter
1500 : 3530 : m_stage = 0;
1501 : : // Compute diagnostics, e.g., residuals
1502 [ + - ]: 3530 : diag();
1503 : :
1504 : : }
1505 : 4660 : }
1506 : :
1507 : : void
1508 : 3530 : LohCG::evalres( const std::vector< tk::real >& )
1509 : : // *****************************************************************************
1510 : : // Evaluate residuals
1511 : : // *****************************************************************************
1512 : : {
1513 : 3530 : refine();
1514 : 3530 : }
1515 : :
1516 : : void
1517 : 3530 : LohCG::refine()
1518 : : // *****************************************************************************
1519 : : // Optionally refine/derefine mesh
1520 : : // *****************************************************************************
1521 : : {
1522 : 3530 : auto d = Disc();
1523 : :
1524 : : // See if this is the last time step
1525 [ + + ]: 3530 : if (d->finished()) m_finished = 1;
1526 : :
1527 : 3530 : auto dtref = g_cfg.get< tag::href_dt >();
1528 : 3530 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
1529 : :
1530 : : // if t>0 refinement enabled and we hit the frequency
1531 [ - + ][ - - ]: 3530 : if (dtref && !(d->It() % dtfreq)) { // refine
1532 : :
1533 : 0 : d->refined() = 1;
1534 : 0 : d->startvol();
1535 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1536 : :
1537 : : // Activate SDAG waits for re-computing the integrals
1538 [ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1539 : :
1540 : : } else { // do not refine
1541 : :
1542 : 3530 : d->refined() = 0;
1543 : 3530 : feop_complete();
1544 : 3530 : resize_complete();
1545 : :
1546 : : }
1547 : 3530 : }
1548 : :
1549 : : void
1550 : 0 : LohCG::resizePostAMR(
1551 : : const std::vector< std::size_t >& /*ginpoel*/,
1552 : : const tk::UnsMesh::Chunk& chunk,
1553 : : const tk::UnsMesh::Coords& coord,
1554 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1555 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1556 : : const std::set< std::size_t >& removedNodes,
1557 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1558 : : nodeCommMap,
1559 : : const std::map< int, std::vector< std::size_t > >& bface,
1560 : : const std::map< int, std::vector< std::size_t > >& bnode,
1561 : : const std::vector< std::size_t >& triinpoel )
1562 : : // *****************************************************************************
1563 : : // Receive new mesh from Refiner
1564 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1565 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1566 : : //! \param[in] coord New mesh node coordinates
1567 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1568 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1569 : : //! \param[in] removedNodes Newly removed mesh node local ids
1570 : : //! \param[in] nodeCommMap New node communication map
1571 : : //! \param[in] bface Boundary-faces mapped to side set ids
1572 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1573 : : //! \param[in] triinpoel Boundary-face connectivity
1574 : : // *****************************************************************************
1575 : : {
1576 : 0 : auto d = Disc();
1577 : :
1578 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1579 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1580 : :
1581 : : // Resize mesh data structures after mesh refinement
1582 : 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1583 : :
1584 : : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
1585 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1586 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1587 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1588 : :
1589 : : // Remove newly removed nodes from solution vectors
1590 : 0 : m_u.rm( removedNodes );
1591 : 0 : m_rhs.rm( removedNodes );
1592 : :
1593 : : // Resize auxiliary solution vectors
1594 : : auto npoin = coord[0].size();
1595 : : m_u.resize( npoin );
1596 : : m_rhs.resize( npoin );
1597 : :
1598 : : // Update solution on new mesh
1599 [ - - ]: 0 : for (const auto& n : addedNodes)
1600 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1601 : : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
1602 : : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
1603 : : "Indices of parent-edge nodes out of bounds post-AMR");
1604 : 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
1605 : : }
1606 : :
1607 : : // Update physical-boundary node-, face-, and element lists
1608 : : m_bnode = bnode;
1609 : : m_bface = bface;
1610 : 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1611 : :
1612 : 0 : auto meshid = d->MeshId();
1613 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1614 : 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1615 : 0 : }
1616 : :
1617 : : void
1618 : 797 : LohCG::writeFields( CkCallback cb )
1619 : : // *****************************************************************************
1620 : : // Output mesh-based fields to file
1621 : : //! \param[in] cb Function to continue with after the write
1622 : : // *****************************************************************************
1623 : : {
1624 [ + + ]: 797 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1625 : :
1626 : 213 : auto d = Disc();
1627 : :
1628 : : // Field output
1629 : :
1630 : : std::vector< std::string > nodefieldnames{
1631 [ - + ][ + + ]: 1065 : "pressure", "velocityx", "velocityy", "velocityz" };
[ - + ][ - - ]
[ - - ]
1632 : :
1633 : : std::vector< std::vector< tk::real > > nodefields;
1634 : :
1635 [ + - ]: 213 : nodefields.push_back( m_u.extract(0) );
1636 [ + - ]: 213 : nodefields.push_back( m_u.extract(1) );
1637 [ + - ]: 213 : nodefields.push_back( m_u.extract(2) );
1638 [ + - ][ - + ]: 426 : nodefields.push_back( m_u.extract(3) );
1639 : :
1640 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
1641 : :
1642 : : // Surface output
1643 : :
1644 : : std::vector< std::string > nodesurfnames;
1645 : : std::vector< std::vector< tk::real > > nodesurfs;
1646 : :
1647 : : const auto& f = g_cfg.get< tag::fieldout >();
1648 : :
1649 [ - + ]: 213 : if (!f.empty()) {
1650 : : std::size_t ncomp = 4;
1651 [ - - ]: 0 : nodesurfnames.push_back( "pressure" );
1652 [ - - ]: 0 : nodesurfnames.push_back( "velocityx" );
1653 [ - - ]: 0 : nodesurfnames.push_back( "velocityy" );
1654 [ - - ]: 0 : nodesurfnames.push_back( "velocityz" );
1655 : :
1656 [ - - ]: 0 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
1657 [ - - ]: 0 : std::set< int > outsets( begin(f), end(f) );
1658 [ - - ]: 0 : for (auto sideset : outsets) {
1659 : : auto b = bnode.find(sideset);
1660 [ - - ]: 0 : if (b == end(bnode)) continue;
1661 : : const auto& nodes = b->second;
1662 : : auto i = nodesurfs.size();
1663 : : nodesurfs.insert( end(nodesurfs), ncomp,
1664 [ - - ]: 0 : std::vector< tk::real >( nodes.size() ) );
1665 : : std::size_t j = 0;
1666 [ - - ]: 0 : for (auto n : nodes) {
1667 [ - - ]: 0 : const auto s = m_u[n];
1668 : : std::size_t p = 0;
1669 : 0 : nodesurfs[i+(p++)][j] = s[0];
1670 : 0 : nodesurfs[i+(p++)][j] = s[1];
1671 : 0 : nodesurfs[i+(p++)][j] = s[2];
1672 : 0 : nodesurfs[i+(p++)][j] = s[3];
1673 : 0 : ++j;
1674 : : }
1675 : : }
1676 : : }
1677 : :
1678 : : // Send mesh and fields data (solution dump) for output to file
1679 [ + - ][ + - ]: 426 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
1680 [ + - ]: 213 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
1681 : : {}, nodefields, {}, nodesurfs, cb );
1682 [ + - ][ + - ]: 639 : }
[ + - ][ + - ]
[ + - ][ - - ]
[ - - ]
1683 : :
1684 : : void
1685 : 3530 : LohCG::out()
1686 : : // *****************************************************************************
1687 : : // Output mesh field data
1688 : : // *****************************************************************************
1689 : : {
1690 : 3530 : auto d = Disc();
1691 : :
1692 : : // Time history
1693 [ + - ][ + - ]: 3530 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
1694 : : auto ncomp = m_u.nprop();
1695 : : const auto& inpoel = d->Inpoel();
1696 : 0 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
1697 : : std::size_t j = 0;
1698 [ - - ]: 0 : for (const auto& p : d->Hist()) {
1699 [ - - ]: 0 : auto e = p.get< tag::elem >(); // host element id
1700 : : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
1701 [ - - ]: 0 : hist[j].resize( ncomp+1, 0.0 );
1702 [ - - ]: 0 : for (std::size_t i=0; i<4; ++i) {
1703 [ - - ]: 0 : const auto u = m_u[ inpoel[e*4+i] ];
1704 : 0 : hist[j][0] += n[i] * u[0];
1705 : 0 : hist[j][1] += n[i] * u[1]/u[0];
1706 : 0 : hist[j][2] += n[i] * u[2]/u[0];
1707 : 0 : hist[j][3] += n[i] * u[3]/u[0];
1708 : 0 : hist[j][4] += n[i] * u[4]/u[0];
1709 : 0 : auto ei = u[4]/u[0] - 0.5*(u[1]*u[1] + u[2]*u[2] + u[3]*u[3])/u[0]/u[0];
1710 : 0 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
1711 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
1712 : : }
1713 : 0 : ++j;
1714 : : }
1715 [ - - ]: 0 : d->history( std::move(hist) );
1716 : 0 : }
1717 : :
1718 : : // Field data
1719 [ + + ][ + - ]: 3530 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
1720 [ + - ][ + - ]: 1416 : writeFields( CkCallback(CkIndex_LohCG::integrals(), thisProxy[thisIndex]) );
1721 : : } else {
1722 : 3058 : integrals();
1723 : : }
1724 : 3530 : }
1725 : :
1726 : : void
1727 : 3530 : LohCG::integrals()
1728 : : // *****************************************************************************
1729 : : // Compute integral quantities for output
1730 : : // *****************************************************************************
1731 : : {
1732 : 3530 : auto d = Disc();
1733 : :
1734 [ + - ][ + - ]: 3530 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ]
1735 : :
1736 : : using namespace integrals;
1737 [ - - ]: 0 : std::vector< std::map< int, tk::real > > ints( NUMINT );
1738 : : // Prepend integral vector with metadata on the current time step:
1739 : : // current iteration count, current physical time, time step size
1740 [ - - ][ - - ]: 0 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
1741 [ - - ][ - - ]: 0 : ints[ TIME ][ 0 ] = d->T();
1742 [ - - ]: 0 : ints[ DT ][ 0 ] = d->Dt();
1743 : : // Compute mass flow rate for surfaces requested
1744 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
1745 : : // cppcheck-suppress unreadVariable
1746 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
1747 : : const auto& nodes = sint.first;
1748 : : const auto& ndA = sint.second;
1749 [ - - ]: 0 : for (std::size_t i=0; i<nodes.size(); ++i) {
1750 : 0 : auto p = nodes[i];
1751 : 0 : mfr += ndA[i*3+0] * m_u(p,1)
1752 : 0 : + ndA[i*3+1] * m_u(p,2)
1753 : 0 : + ndA[i*3+2] * m_u(p,3);
1754 : : }
1755 : : }
1756 [ - - ]: 0 : auto stream = serialize( d->MeshId(), ints );
1757 [ - - ][ - - ]: 0 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
1758 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
1759 : :
1760 : 0 : } else {
1761 : :
1762 : 3530 : step();
1763 : :
1764 : : }
1765 : 3530 : }
1766 : :
1767 : : void
1768 : 3205 : LohCG::evalLB( int nrestart )
1769 : : // *****************************************************************************
1770 : : // Evaluate whether to do load balancing
1771 : : //! \param[in] nrestart Number of times restarted
1772 : : // *****************************************************************************
1773 : : {
1774 : 3205 : auto d = Disc();
1775 : :
1776 : : // Detect if just returned from a checkpoint and if so, zero timers and
1777 : : // finished flag
1778 [ + + ]: 3205 : if (d->restarted( nrestart )) m_finished = 0;
1779 : :
1780 : 3205 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
1781 [ + + ]: 3205 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
1782 : :
1783 : : // Load balancing if user frequency is reached or after the second time-step
1784 [ + + ][ + + ]: 3205 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
1785 : :
1786 : 3065 : AtSync();
1787 [ - + ]: 3065 : if (nonblocking) dt();
1788 : :
1789 : : } else {
1790 : :
1791 : 140 : dt();
1792 : :
1793 : : }
1794 : 3205 : }
1795 : :
1796 : : void
1797 : 3200 : LohCG::evalRestart()
1798 : : // *****************************************************************************
1799 : : // Evaluate whether to save checkpoint/restart
1800 : : // *****************************************************************************
1801 : : {
1802 : 3200 : auto d = Disc();
1803 : :
1804 : 3200 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1805 : 3200 : const auto benchmark = g_cfg.get< tag::benchmark >();
1806 : :
1807 [ + + ][ - + ]: 3200 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
1808 : :
1809 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1810 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1811 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
1812 : :
1813 : : } else {
1814 : :
1815 : 3200 : evalLB( /* nrestart = */ -1 );
1816 : :
1817 : : }
1818 : 3200 : }
1819 : :
1820 : : void
1821 : 3530 : LohCG::step()
1822 : : // *****************************************************************************
1823 : : // Evaluate whether to continue with next time step
1824 : : // *****************************************************************************
1825 : : {
1826 : 3530 : auto d = Disc();
1827 : :
1828 : : // Output one-liner status report to screen
1829 [ + + ]: 3530 : if(thisIndex == 0) d->status();
1830 : :
1831 [ + + ]: 3530 : if (not m_finished) {
1832 : :
1833 : 3200 : evalRestart();
1834 : :
1835 : : } else {
1836 : :
1837 : 330 : auto meshid = d->MeshId();
1838 [ + - ]: 660 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1839 : 330 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1840 : :
1841 : : }
1842 : 3530 : }
1843 : :
1844 : : #include "NoWarning/lohcg.def.h"
|
