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