Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ChoCG.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 ChoCG: Projection-based solver for incompressible flow
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "XystBuildConfig.hpp"
14 : : #include "ChoCG.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 "Chorin.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 : : //! Runge-Kutta coefficients
43 : : static const std::array< std::vector< tk::real >, 4 > rkcoef{{
44 : : { 1.0 },
45 : : { 1.0/2.0, 1.0 },
46 : : { 1.0/3.0, 1.0/2.0, 1.0 },
47 : : { 1.0/4.0, 1.0/3.0, 1.0/2.0, 1.0 }
48 : : }};
49 : :
50 : : } // inciter::
51 : :
52 : : using inciter::g_cfg;
53 : : using inciter::ChoCG;
54 : :
55 : 365 : ChoCG::ChoCG( const CProxy_Discretization& disc,
56 : : const tk::CProxy_ConjugateGradients& cgpre,
57 : : const tk::CProxy_ConjugateGradients& cgmom,
58 : : const std::map< int, std::vector< std::size_t > >& bface,
59 : : const std::map< int, std::vector< std::size_t > >& bnode,
60 : 365 : const std::vector< std::size_t >& triinpoel ) :
61 [ + - ]: 365 : m_disc( disc ),
62 [ + - ]: 365 : m_cgpre( cgpre ),
63 [ + - ]: 365 : m_cgmom( cgmom ),
64 : 365 : m_nrhs( 0 ),
65 : 365 : m_nnorm( 0 ),
66 : 365 : m_naec( 0 ),
67 : 365 : m_nalw( 0 ),
68 : 365 : m_nlim( 0 ),
69 : 365 : m_nsgrad( 0 ),
70 : 365 : m_npgrad( 0 ),
71 : 365 : m_nvgrad( 0 ),
72 : 365 : m_nflux( 0 ),
73 : 365 : m_ndiv( 0 ),
74 : 365 : m_nbpint( 0 ),
75 : 365 : m_np( 0 ),
76 [ + - ]: 365 : m_bnode( bnode ),
77 [ + - ]: 365 : m_bface( bface ),
78 [ + - ][ + - ]: 365 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
79 [ + - ][ + - ]: 365 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
80 [ + - ]: 365 : m_un( m_u.nunk(), m_u.nprop() ),
81 [ + - ]: 365 : m_pr( m_u.nunk(), 0.0 ),
82 [ + - ]: 365 : m_p( m_u.nunk(), m_u.nprop()*2 ),
83 [ + - ]: 365 : m_q( m_u.nunk(), m_u.nprop()*2 ),
84 [ + - ]: 365 : m_a( m_u.nunk(), m_u.nprop() ),
85 [ + - ]: 365 : m_rhs( m_u.nunk(), m_u.nprop() ),
86 [ + - ]: 365 : m_sgrad( m_u.nunk(), 3UL ),
87 [ + - ]: 365 : m_pgrad( m_u.nunk(), 3UL ),
88 [ + - ]: 365 : m_vgrad( m_u.nunk(), 9UL ),
89 [ + - ]: 365 : m_flux( m_u.nunk(), 3UL ),
90 [ + - ]: 365 : m_div( m_u.nunk() ),
91 : 365 : m_stage( 0 ),
92 : 365 : m_finished( 0 ),
93 [ + - ]: 3650 : m_rkcoef( rkcoef[ g_cfg.get< tag::rk >() - 1 ] )
94 : : // *****************************************************************************
95 : : // Constructor
96 : : //! \param[in] disc Discretization proxy
97 : : //! \param[in] cgpre ConjugateGradients Charm++ proxy for pressure solve
98 : : //! \param[in] cgmom ConjugateGradients Charm++ proxy for momentum solve
99 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
100 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
101 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
102 : : // *****************************************************************************
103 : : {
104 : 365 : usesAtSync = true; // enable migration at AtSync
105 : :
106 [ + - ]: 365 : auto d = Disc();
107 : :
108 : : // Create new local ids based on mesh locality
109 : 365 : std::unordered_map< std::size_t, std::size_t > map;
110 : 365 : std::size_t n = 0;
111 : :
112 [ + - ][ + - ]: 365 : auto psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
113 [ + + ]: 25667 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
114 [ + - ][ + + ]: 25302 : if (!map.count(p)) map[p] = n++;
[ + - ]
115 [ + + ]: 240280 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
116 [ + - ][ + + ]: 214978 : if (!map.count(q)) map[q] = n++;
[ + - ]
117 : : }
118 : : }
119 : :
120 [ - + ][ - - ]: 365 : Assert( map.size() == d->Gid().size(),
[ - - ][ - - ]
121 : : "Mesh-locality reorder map size mismatch" );
122 : :
123 : : // Remap data in bound Discretization object
124 [ + - ]: 365 : d->remap( map );
125 : : // Remap boundary triangle face connectivity
126 [ + - ]: 365 : tk::remap( m_triinpoel, map );
127 : : // Recompute points surrounding points
128 [ + - ][ + - ]: 365 : psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
129 : :
130 : : // Compute total box IC volume
131 [ + - ]: 365 : d->boxvol();
132 : :
133 : : // Setup LHS matrix for pressure solve
134 [ + - ][ + - ]: 365 : m_cgpre[ thisIndex ].insert( prelhs( psup ),
[ + - ]
135 : : d->Gid(),
136 : : d->Lid(),
137 : : d->NodeCommMap() );
138 : :
139 : : // Setup empty LHS matrix for momentum solve if needed
140 [ + + ]: 365 : if (g_cfg.get< tag::theta >() > std::numeric_limits< tk::real >::epsilon()) {
141 [ + - ][ + - ]: 2 : m_cgmom[ thisIndex ].insert( momlhs( psup ),
[ + - ]
142 : : d->Gid(),
143 : : d->Lid(),
144 : : d->NodeCommMap() );
145 : : }
146 : :
147 : : // Activate SDAG waits for setup
148 [ + - ][ + - ]: 365 : thisProxy[ thisIndex ].wait4int();
149 : 365 : }
150 : :
151 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
152 : 365 : ChoCG::prelhs( const std::pair< std::vector< std::size_t >,
153 : : std::vector< std::size_t > >& psup )
154 : : // *****************************************************************************
155 : : // Setup lhs matrix for pressure solve
156 : : //! \param[in] psup Points surrounding points
157 : : //! \return { A, x, b } in linear system A * x = b to solve for pressure
158 : : // *****************************************************************************
159 : : {
160 [ + - ]: 365 : auto d = Disc();
161 : 365 : const auto& inpoel = d->Inpoel();
162 : 365 : const auto& coord = d->Coord();
163 : 365 : const auto& X = coord[0];
164 : 365 : const auto& Y = coord[1];
165 : 365 : const auto& Z = coord[2];
166 : :
167 : : // Matrix with compressed sparse row storage
168 [ + - ]: 365 : tk::CSR A( /*DOF=*/ 1, psup );
169 : :
170 : : // Fill matrix with Laplacian
171 [ + + ]: 59289 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
172 : 58924 : const auto N = inpoel.data() + e*4;
173 : : const std::array< tk::real, 3 >
174 : 58924 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
175 : 58924 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
176 : 58924 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
177 : 58924 : const auto J = tk::triple( ba, ca, da ) * 6.0;
178 : : std::array< std::array< tk::real, 3 >, 4 > grad;
179 : 58924 : grad[1] = tk::cross( ca, da );
180 : 58924 : grad[2] = tk::cross( da, ba );
181 : 58924 : grad[3] = tk::cross( ba, ca );
182 [ + + ]: 235696 : for (std::size_t i=0; i<3; ++i)
183 : 176772 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
184 [ + + ]: 294620 : for (std::size_t a=0; a<4; ++a)
185 [ + + ]: 1178480 : for (std::size_t b=0; b<4; ++b)
186 [ + - ]: 942784 : A(N[a],N[b]) -= tk::dot( grad[a], grad[b] ) / J;
187 : : }
188 : :
189 : 365 : auto nunk = X.size();
190 [ + - ][ + - ]: 365 : std::vector< tk::real > x( nunk, 0.0 ), b( nunk, 0.0 );
191 : :
192 : 730 : return { std::move(A), std::move(x), std::move(b) };
193 : 365 : }
194 : :
195 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
196 : 2 : ChoCG::momlhs( const std::pair< std::vector< std::size_t >,
197 : : std::vector< std::size_t > >& psup )
198 : : // *****************************************************************************
199 : : // Setup empty lhs matrix for momentum solve
200 : : //! \param[in] psup Points surrounding points
201 : : //! \return { A, x, b } in linear system A * x = b to solve for momentum
202 : : // *****************************************************************************
203 : : {
204 : 2 : auto ncomp = m_u.nprop();
205 : :
206 : : // Matrix with compressed sparse row storage
207 [ + - ]: 2 : tk::CSR A( /*DOF=*/ ncomp, psup );
208 : :
209 : 2 : auto nunk = (psup.second.size() - 1) * ncomp;
210 [ + - ][ + - ]: 2 : std::vector< tk::real > x( nunk, 0.0 ), b( nunk, 0.0 );
211 : :
212 : 4 : return { std::move(A), std::move(x), std::move(b) };
213 : 2 : }
214 : :
215 : : void
216 : 730 : ChoCG::setupDirBC( const std::vector< std::vector< int > >& cfgmask,
217 : : const std::vector< std::vector< double > >& cfgval,
218 : : std::size_t ncomp,
219 : : std::vector< std::size_t >& mask,
220 : : std::vector< double >& val )
221 : : // *****************************************************************************
222 : : // Prepare Dirichlet boundary condition data structures
223 : : //! \param[in] cfgmask Boundary condition mask config data to use
224 : : //! \param[in] cfgval Boundary condition values config data to use
225 : : //! \param[in] ncomp Number of scalar component BCs expected per mesh node
226 : : //! \param[in,out] mask Mesh nodes and their Dirichlet BC masks
227 : : //! \param[in,out] val Mesh nodes and their Dirichlet BC values
228 : : // *****************************************************************************
229 : : {
230 : : // Query Dirichlet BC nodes associated to side sets
231 : 730 : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
232 [ + + ]: 1052 : for (const auto& s : cfgmask) {
233 [ + - ]: 322 : auto k = m_bface.find(s[0]);
234 [ + + ]: 322 : if (k != end(m_bface)) {
235 [ + - ]: 176 : auto& n = dir[ k->first ];
236 [ + + ]: 23596 : for (auto f : k->second) {
237 [ + - ]: 23420 : n.insert( m_triinpoel[f*3+0] );
238 [ + - ]: 23420 : n.insert( m_triinpoel[f*3+1] );
239 [ + - ]: 23420 : n.insert( m_triinpoel[f*3+2] );
240 : : }
241 : : }
242 : : }
243 : :
244 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
245 [ + - ]: 730 : const auto& lid = Disc()->Lid();
246 [ + + ]: 1052 : for (const auto& s : cfgmask) {
247 [ + - ]: 322 : auto k = m_bnode.find(s[0]);
248 [ + + ]: 322 : if (k != end(m_bnode)) {
249 [ + - ]: 180 : auto& n = dir[ k->first ];
250 [ + + ]: 17746 : for (auto g : k->second) {
251 [ + - ][ + - ]: 17566 : n.insert( tk::cref_find(lid,g) );
252 : : }
253 : : }
254 : : }
255 : :
256 : : // Associate sidesets to Dirichlet BC values if configured by user
257 : 730 : std::unordered_map< int, std::vector< double > > dirval;
258 [ + + ]: 796 : for (const auto& s : cfgval) {
259 [ + - ]: 66 : auto k = dir.find( static_cast<int>(s[0]) );
260 [ + + ]: 66 : if (k != end(dir)) {
261 [ + - ]: 24 : auto& v = dirval[ k->first ];
262 [ + - ]: 24 : v.resize( s.size()-1 );
263 [ + + ]: 56 : for (std::size_t i=1; i<s.size(); ++i) v[i-1] = s[i];
264 : : }
265 : : }
266 : :
267 : : // Collect unique set of nodes + Dirichlet BC components mask and value
268 : 730 : auto nmask = ncomp + 1;
269 : : std::unordered_map< std::size_t,
270 : : std::pair< std::vector< int >,
271 : 730 : std::vector< double > > > dirbcset;
272 [ + + ]: 1052 : for (const auto& vec : cfgmask) {
273 [ - + ][ - - ]: 322 : ErrChk( vec.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
274 [ + - ]: 322 : auto n = dir.find( vec[0] );
275 [ + + ]: 322 : if (n != end(dir)) {
276 [ + - ]: 180 : std::vector< double > v( ncomp, 0.0 );
277 [ + - ]: 180 : auto m = dirval.find( vec[0] );
278 [ + + ]: 180 : if (m != end(dirval)) {
279 [ - + ][ - - ]: 24 : ErrChk( m->second.size() == ncomp, "Incorrect Dirichlet BC val ncomp" );
[ - - ][ - - ]
280 [ + - ]: 24 : v = m->second;
281 : : }
282 [ + + ]: 15341 : for (auto p : n->second) {
283 [ + - ]: 15161 : auto& mv = dirbcset[p]; // mask & value
284 [ + - ]: 15161 : mv.second = v;
285 : 15161 : auto& mval = mv.first;
286 [ + + ][ + - ]: 15161 : if (mval.empty()) mval.resize( ncomp, 0 );
287 [ + + ]: 55758 : for (std::size_t c=0; c<ncomp; ++c)
288 [ + + ]: 40597 : if (!mval[c]) mval[c] = vec[c+1]; // overwrite mask if 0 -> 1
289 : : }
290 : 180 : }
291 : : }
292 : :
293 : : // Compile streamable list of nodes + Dirichlet BC components mask and values
294 : 730 : tk::destroy( mask );
295 [ + + ]: 15329 : for (const auto& [p,mv] : dirbcset) {
296 [ + - ]: 14599 : mask.push_back( p );
297 [ + - ]: 14599 : mask.insert( end(mask), begin(mv.first), end(mv.first) );
298 [ + - ]: 14599 : val.push_back( static_cast< double >( p ) );
299 [ + - ]: 14599 : val.insert( end(val), begin(mv.second), end(mv.second) );
300 : : }
301 : :
302 [ - + ][ - - ]: 730 : ErrChk( mask.size() % nmask == 0, "Dirichlet BC mask incomplete" );
[ - - ][ - - ]
303 [ - + ][ - - ]: 730 : ErrChk( val.size() % nmask == 0, "Dirichlet BC val incomplete" );
[ - - ][ - - ]
304 [ - + ][ - - ]: 730 : ErrChk( mask.size() == val.size(), "Dirichlet BC mask & val size mismatch" );
[ - - ][ - - ]
305 : 730 : }
306 : :
307 : : void
308 : 365 : ChoCG::feop()
309 : : // *****************************************************************************
310 : : // Start (re-)computing finite element domain and boundary operators
311 : : // *****************************************************************************
312 : : {
313 : 365 : auto d = Disc();
314 : :
315 : : // Prepare Dirichlet boundary conditions data structures
316 : 365 : setupDirBC( g_cfg.get< tag::bc_dir >(), g_cfg.get< tag::bc_dirval >(),
317 : 365 : m_u.nprop(), m_dirbcmask, m_dirbcval );
318 : 365 : setupDirBC( g_cfg.get< tag::pre_bc_dir >(), g_cfg.get< tag::pre_bc_dirval >(),
319 : 365 : 1, m_dirbcmaskp, m_dirbcvalp );
320 : :
321 : : // Compute local contributions to boundary normals and integrals
322 : 365 : bndint();
323 : : // Compute local contributions to domain edge integrals
324 : 365 : domint();
325 : :
326 : : // Send boundary point normal contributions to neighbor chares
327 [ + + ]: 365 : if (d->NodeCommMap().empty()) {
328 : 12 : comnorm_complete();
329 : : } else {
330 [ + + ]: 4093 : for (const auto& [c,nodes] : d->NodeCommMap()) {
331 : 3740 : decltype(m_bnorm) exp;
332 [ + + ]: 19386 : for (auto i : nodes) {
333 [ + + ]: 49082 : for (const auto& [s,b] : m_bnorm) {
334 [ + - ]: 33436 : auto k = b.find(i);
335 [ + + ][ + - ]: 33436 : if (k != end(b)) exp[s][i] = k->second;
[ + - ]
336 : : }
337 : : }
338 [ + - ][ + - ]: 3740 : thisProxy[c].comnorm( exp );
339 : 3740 : }
340 : : }
341 : 365 : ownnorm_complete();
342 : 365 : }
343 : :
344 : : void
345 : 365 : ChoCG::bndint()
346 : : // *****************************************************************************
347 : : // Compute local contributions to boundary normals and integrals
348 : : // *****************************************************************************
349 : : {
350 [ + - ]: 365 : auto d = Disc();
351 : 365 : const auto& coord = d->Coord();
352 : 365 : const auto& gid = d->Gid();
353 : 365 : const auto& x = coord[0];
354 : 365 : const auto& y = coord[1];
355 : 365 : const auto& z = coord[2];
356 : :
357 : : // Lambda to compute the inverse distance squared between boundary face
358 : : // centroid and boundary point. Here p is the global node id and c is the
359 : : // the boundary face centroid.
360 : 114324 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
361 : 114324 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
362 : 114324 : (c[1] - y[p]) * (c[1] - y[p]) +
363 : 114324 : (c[2] - z[p]) * (c[2] - z[p]) );
364 : 365 : };
365 : :
366 : 365 : tk::destroy( m_bnorm );
367 : 365 : tk::destroy( m_bndpoinint );
368 : :
369 [ + + ]: 1239 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
370 [ + + ]: 38982 : for (auto f : faceids) { // for all side set triangles
371 : 38108 : const auto N = m_triinpoel.data() + f*3;
372 : : const std::array< tk::real, 3 >
373 : 38108 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
374 : 38108 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
375 : 38108 : auto n = tk::cross( ba, ca );
376 : 38108 : auto A2 = tk::length( n );
377 : 38108 : n[0] /= A2;
378 : 38108 : n[1] /= A2;
379 : 38108 : n[2] /= A2;
380 : : const tk::real centroid[3] = {
381 : 38108 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
382 : 38108 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
383 : 76216 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
384 [ + + ]: 152432 : for (const auto& [i,j] : tk::lpoet) {
385 : 114324 : auto p = N[i];
386 : 114324 : tk::real r = invdistsq( centroid, p );
387 [ + - ]: 114324 : auto& v = m_bnorm[setid]; // associate side set id
388 [ + - ]: 114324 : auto& bpn = v[gid[p]]; // associate global node id of bnd pnt
389 : 114324 : bpn[0] += r * n[0]; // inv.dist.sq-weighted normal
390 : 114324 : bpn[1] += r * n[1];
391 : 114324 : bpn[2] += r * n[2];
392 : 114324 : bpn[3] += r; // inv.dist.sq of node from centroid
393 [ + - ]: 114324 : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
394 : 114324 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
395 : 114324 : b[1] += n[1] * A2 / 6.0;
396 : 114324 : b[2] += n[2] * A2 / 6.0;
397 : : }
398 : : }
399 : : }
400 : 365 : }
401 : :
402 : : void
403 : 365 : ChoCG::domint()
404 : : // *****************************************************************************
405 : : //! Compute local contributions to domain edge integrals
406 : : // *****************************************************************************
407 : : {
408 : 365 : auto d = Disc();
409 : :
410 : 365 : const auto& gid = d->Gid();
411 : 365 : const auto& inpoel = d->Inpoel();
412 : :
413 : 365 : const auto& coord = d->Coord();
414 : 365 : const auto& x = coord[0];
415 : 365 : const auto& y = coord[1];
416 : 365 : const auto& z = coord[2];
417 : :
418 : 365 : tk::destroy( m_domedgeint );
419 : :
420 [ + + ]: 59289 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
421 : 58924 : const auto N = inpoel.data() + e*4;
422 : : const std::array< tk::real, 3 >
423 : 58924 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
424 : 58924 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
425 : 58924 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
426 : 58924 : const auto J = tk::triple( ba, ca, da ); // J = 6V
427 [ - + ][ - - ]: 58924 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
428 : : std::array< std::array< tk::real, 3 >, 4 > grad;
429 : 58924 : grad[1] = tk::cross( ca, da );
430 : 58924 : grad[2] = tk::cross( da, ba );
431 : 58924 : grad[3] = tk::cross( ba, ca );
432 [ + + ]: 235696 : for (std::size_t i=0; i<3; ++i)
433 : 176772 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
434 [ + + ]: 412468 : for (const auto& [p,q] : tk::lpoed) {
435 : 353544 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
436 : 353544 : tk::real sig = 1.0;
437 [ + + ]: 353544 : if (ed[0] > ed[1]) {
438 : 113856 : std::swap( ed[0], ed[1] );
439 : 113856 : sig = -1.0;
440 : : }
441 [ + - ]: 353544 : auto& n = m_domedgeint[ ed ];
442 : 353544 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
443 : 353544 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
444 : 353544 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
445 : 353544 : n[3] += J / 120.0;
446 : 353544 : n[4] += tk::dot( grad[p], grad[q] ) / J / 6.0;
447 : : }
448 : : }
449 : 365 : }
450 : :
451 : : void
452 : 3740 : ChoCG::comnorm( const decltype(m_bnorm)& inbnd )
453 : : // *****************************************************************************
454 : : // Receive contributions to boundary point normals on chare-boundaries
455 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
456 : : // *****************************************************************************
457 : : {
458 : : // Buffer up incoming boundary point normal vector contributions
459 [ + + ]: 6635 : for (const auto& [s,b] : inbnd) {
460 [ + - ]: 2895 : auto& bndnorm = m_bnormc[s];
461 [ + + ]: 10922 : for (const auto& [p,n] : b) {
462 [ + - ]: 8027 : auto& norm = bndnorm[p];
463 : 8027 : norm[0] += n[0];
464 : 8027 : norm[1] += n[1];
465 : 8027 : norm[2] += n[2];
466 : 8027 : norm[3] += n[3];
467 : : }
468 : : }
469 : :
470 [ + + ]: 3740 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
471 : 353 : m_nnorm = 0;
472 : 353 : comnorm_complete();
473 : : }
474 : 3740 : }
475 : :
476 : : void
477 : 783 : ChoCG::registerReducers()
478 : : // *****************************************************************************
479 : : // Configure Charm++ reduction types initiated from this chare array
480 : : //! \details Since this is a [initnode] routine, the runtime system executes the
481 : : //! routine exactly once on every logical node early on in the Charm++ init
482 : : //! sequence. Must be static as it is called without an object. See also:
483 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
484 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
485 : : // *****************************************************************************
486 : : {
487 : 783 : NodeDiagnostics::registerReducers();
488 : 783 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
489 : 783 : }
490 : :
491 : : void
492 : : // cppcheck-suppress unusedFunction
493 : 3157 : ChoCG::ResumeFromSync()
494 : : // *****************************************************************************
495 : : // Return from migration
496 : : //! \details This is called when load balancing (LB) completes. The presence of
497 : : //! this function does not affect whether or not we block on LB.
498 : : // *****************************************************************************
499 : : {
500 [ - + ][ - - ]: 3157 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
501 : :
502 [ + - ]: 3157 : if (!g_cfg.get< tag::nonblocking >()) dt();
503 : 3157 : }
504 : :
505 : : void
506 : 365 : ChoCG::setup( tk::real v )
507 : : // *****************************************************************************
508 : : // Start setup for solution
509 : : //! \param[in] v Total volume within user-specified box
510 : : // *****************************************************************************
511 : : {
512 : 365 : auto d = Disc();
513 : :
514 : : // Store user-defined box IC volume
515 : 365 : Disc()->Boxvol() = v;
516 : :
517 : : // Set initial conditions
518 : 365 : problems::initialize( d->Coord(), m_u, d->T(), d->BoxNodes() );
519 : :
520 : : // Query time history field output labels from all PDEs integrated
521 [ - + ]: 365 : if (!g_cfg.get< tag::histout >().empty()) {
522 : : std::vector< std::string > var
523 [ - - ][ - - ]: 0 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ]
524 : 0 : auto ncomp = m_u.nprop();
525 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c)
526 [ - - ][ - - ]: 0 : var.push_back( "c" + std::to_string(c-5) );
[ - - ]
527 [ - - ]: 0 : d->histheader( std::move(var) );
528 : 0 : }
529 : :
530 : : // Compute finite element operators
531 : 365 : feop();
532 [ - - ][ - - ]: 365 : }
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
533 : :
534 : : void
535 : 365 : ChoCG::bnorm()
536 : : // *****************************************************************************
537 : : // Combine own and communicated portions of the boundary point normals
538 : : // *****************************************************************************
539 : : {
540 [ + - ]: 365 : const auto& lid = Disc()->Lid();
541 : :
542 : : // Combine own and communicated contributions to boundary point normals
543 [ + + ]: 1212 : for (const auto& [s,b] : m_bnormc) {
544 [ + - ]: 847 : auto& bndnorm = m_bnorm[s];
545 [ + + ]: 6653 : for (const auto& [g,n] : b) {
546 [ + - ]: 5806 : auto& norm = bndnorm[g];
547 : 5806 : norm[0] += n[0];
548 : 5806 : norm[1] += n[1];
549 : 5806 : norm[2] += n[2];
550 : 5806 : norm[3] += n[3];
551 : : }
552 : : }
553 : 365 : tk::destroy( m_bnormc );
554 : :
555 : : // Divide summed point normals by the sum of the inverse distance squared
556 [ + + ]: 1298 : for (auto& [s,b] : m_bnorm) {
557 [ + + ]: 29513 : for (auto& [g,n] : b) {
558 : 28580 : n[0] /= n[3];
559 : 28580 : n[1] /= n[3];
560 : 28580 : n[2] /= n[3];
561 [ - + ][ - - ]: 28580 : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
[ - - ][ - - ]
562 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
563 : : "Non-unit normal" );
564 : : }
565 : : }
566 : :
567 : : // Replace global->local ids associated to boundary point normals
568 : 365 : decltype(m_bnorm) loc;
569 [ + + ]: 1298 : for (auto& [s,b] : m_bnorm) {
570 [ + - ]: 933 : auto& bnd = loc[s];
571 [ + + ]: 29513 : for (auto&& [g,n] : b) {
572 [ + - ][ + - ]: 28580 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
573 : : }
574 : : }
575 : 365 : m_bnorm = std::move(loc);
576 : 365 : }
577 : :
578 : : void
579 : 365 : ChoCG::streamable()
580 : : // *****************************************************************************
581 : : // Convert integrals into streamable data structures
582 : : // *****************************************************************************
583 : : {
584 : : // Query surface integral output nodes
585 : 365 : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
586 : 365 : const auto& is = g_cfg.get< tag::integout >();
587 [ + - ]: 365 : std::set< int > outsets( begin(is), end(is) );
588 [ - + ]: 365 : for (auto s : outsets) {
589 [ - - ]: 0 : auto m = m_bface.find(s);
590 [ - - ]: 0 : if (m != end(m_bface)) {
591 [ - - ]: 0 : auto& n = surfintnodes[ m->first ]; // associate set id
592 [ - - ]: 0 : for (auto f : m->second) { // face ids on side set
593 : 0 : auto t = m_triinpoel.data() + f*3;
594 [ - - ]: 0 : n.push_back( t[0] ); // nodes on side set
595 [ - - ]: 0 : n.push_back( t[1] );
596 [ - - ]: 0 : n.push_back( t[2] );
597 : : }
598 : : }
599 : : }
600 [ - - ][ - + ]: 365 : for (auto& [s,n] : surfintnodes) tk::unique( n );
601 : : // Prepare surface integral data
602 : 365 : tk::destroy( m_surfint );
603 [ + - ]: 365 : const auto& gid = Disc()->Gid();
604 [ - + ]: 365 : for (auto&& [s,n] : surfintnodes) {
605 [ - - ]: 0 : auto& sint = m_surfint[s]; // associate set id
606 : 0 : auto& nodes = sint.first;
607 : 0 : auto& ndA = sint.second;
608 : 0 : nodes = std::move(n);
609 [ - - ]: 0 : ndA.resize( nodes.size()*3 );
610 : 0 : auto a = ndA.data();
611 [ - - ]: 0 : for (auto p : nodes) {
612 [ - - ]: 0 : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
613 : 0 : a[0] = b[0]; // store ni * dA
614 : 0 : a[1] = b[1];
615 : 0 : a[2] = b[2];
616 : 0 : a += 3;
617 : : }
618 : : }
619 : 365 : tk::destroy( m_bndpoinint );
620 : :
621 : : // Generate domain superedges
622 [ + - ]: 365 : domsuped();
623 : :
624 : : // Prepare symmetry boundary condition data structures
625 : :
626 : : // Query symmetry BC nodes associated to side sets
627 : 365 : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
628 [ + + ]: 667 : for (auto s : g_cfg.get< tag::bc_sym >()) {
629 [ + - ]: 302 : auto k = m_bface.find(s);
630 [ + + ]: 302 : if (k != end(m_bface)) {
631 [ + - ]: 119 : auto& n = sym[ k->first ];
632 [ + + ]: 3563 : for (auto f : k->second) {
633 : 3444 : const auto& t = m_triinpoel.data() + f*3;
634 [ + - ]: 3444 : n.insert( t[0] );
635 [ + - ]: 3444 : n.insert( t[1] );
636 [ + - ]: 3444 : n.insert( t[2] );
637 : : }
638 : : }
639 : : }
640 : :
641 : : // Generate unique set of symmetry BC nodes of all symmetryc BC side sets
642 : 365 : std::set< std::size_t > symbcnodeset;
643 [ + - ][ + + ]: 484 : for (const auto& [s,n] : sym) symbcnodeset.insert( begin(n), end(n) );
644 : :
645 : : // Generate symmetry BC data as streamable data structures
646 : 365 : tk::destroy( m_symbcnodes );
647 : 365 : tk::destroy( m_symbcnorms );
648 [ + + ]: 2831 : for (auto p : symbcnodeset) {
649 [ + + ]: 6180 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
650 [ + - ]: 3714 : auto m = m_bnorm.find( s );
651 [ + - ]: 3714 : if (m != end(m_bnorm)) {
652 [ + - ]: 3714 : auto r = m->second.find( p );
653 [ + + ]: 3714 : if (r != end(m->second)) {
654 [ + - ]: 2466 : m_symbcnodes.push_back( p );
655 [ + - ]: 2466 : m_symbcnorms.push_back( r->second[0] );
656 [ + - ]: 2466 : m_symbcnorms.push_back( r->second[1] );
657 [ + - ]: 2466 : m_symbcnorms.push_back( r->second[2] );
658 : : }
659 : : }
660 : : }
661 : : }
662 : 365 : tk::destroy( m_bnorm );
663 : :
664 : : // Prepare noslip boundary condition data structures
665 : :
666 : : // Query noslip BC nodes associated to side sets
667 : 365 : std::unordered_map< int, std::unordered_set< std::size_t > > noslip;
668 [ + + ]: 455 : for (auto s : g_cfg.get< tag::bc_noslip >()) {
669 [ + - ]: 90 : auto k = m_bface.find(s);
670 [ + + ]: 90 : if (k != end(m_bface)) {
671 [ + - ]: 84 : auto& n = noslip[ k->first ];
672 [ + + ]: 5104 : for (auto f : k->second) {
673 : 5020 : const auto& t = m_triinpoel.data() + f*3;
674 [ + - ]: 5020 : n.insert( t[0] );
675 [ + - ]: 5020 : n.insert( t[1] );
676 [ + - ]: 5020 : n.insert( t[2] );
677 : : }
678 : : }
679 : : }
680 : :
681 : : // Generate unique set of noslip BC nodes of all noslip BC side sets
682 : 365 : std::set< std::size_t > noslipbcnodeset;
683 [ + - ][ + + ]: 449 : for (const auto& [s,n] : noslip) noslipbcnodeset.insert( begin(n), end(n) );
684 : :
685 : : // Generate noslip BC data as streamable data structures
686 : 365 : tk::destroy( m_noslipbcnodes );
687 [ + - ]: 365 : m_noslipbcnodes.insert( m_noslipbcnodes.end(),
688 : : begin(noslipbcnodeset), end(noslipbcnodeset) );
689 : 365 : }
690 : :
691 : : void
692 : 365 : ChoCG::domsuped()
693 : : // *****************************************************************************
694 : : // Generate superedge-groups for domain-edge loops
695 : : //! \see See Lohner, Sec. 15.1.6.2, An Introduction to Applied CFD Techniques,
696 : : //! Wiley, 2008.
697 : : // *****************************************************************************
698 : : {
699 [ - + ][ - - ]: 365 : Assert( !m_domedgeint.empty(), "No domain edges to group" );
[ - - ][ - - ]
700 : :
701 : : #ifndef NDEBUG
702 : 365 : auto nedge = m_domedgeint.size();
703 : : #endif
704 : :
705 [ + - ]: 365 : const auto& inpoel = Disc()->Inpoel();
706 [ + - ]: 365 : const auto& lid = Disc()->Lid();
707 [ + - ]: 365 : const auto& gid = Disc()->Gid();
708 : :
709 : 365 : tk::destroy( m_dsupedge[0] );
710 : 365 : tk::destroy( m_dsupedge[1] );
711 : 365 : tk::destroy( m_dsupedge[2] );
712 : :
713 : 365 : tk::destroy( m_dsupint[0] );
714 : 365 : tk::destroy( m_dsupint[1] );
715 : 365 : tk::destroy( m_dsupint[2] );
716 : :
717 : 365 : tk::UnsMesh::FaceSet untri;
718 [ + + ]: 59289 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
719 : : std::size_t N[4] = {
720 : 58924 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
721 [ + - ][ + + ]: 294620 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
722 : : }
723 : :
724 [ + + ]: 59289 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
725 : : std::size_t N[4] = {
726 : 58924 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
727 : 58924 : int f = 0;
728 : : tk::real sig[6];
729 [ + - ][ + + ]: 412468 : decltype(m_domedgeint)::const_iterator d[6];
730 [ + + ]: 170622 : for (const auto& [p,q] : tk::lpoed) {
731 : 161150 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
732 [ + + ]: 161150 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
733 [ + - ]: 161150 : d[f] = m_domedgeint.find( ed );
734 [ + + ]: 161150 : if (d[f] == end(m_domedgeint)) break; else ++f;
735 : : }
736 [ + + ]: 58924 : if (f == 6) {
737 [ + - ]: 9472 : m_dsupedge[0].push_back( N[0] );
738 [ + - ]: 9472 : m_dsupedge[0].push_back( N[1] );
739 [ + - ]: 9472 : m_dsupedge[0].push_back( N[2] );
740 [ + - ]: 9472 : m_dsupedge[0].push_back( N[3] );
741 [ + - ][ + + ]: 47360 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
742 [ + + ]: 66304 : for (int ed=0; ed<6; ++ed) {
743 : 56832 : const auto& ded = d[ed]->second;
744 [ + - ]: 56832 : m_dsupint[0].push_back( sig[ed] * ded[0] );
745 [ + - ]: 56832 : m_dsupint[0].push_back( sig[ed] * ded[1] );
746 [ + - ]: 56832 : m_dsupint[0].push_back( sig[ed] * ded[2] );
747 [ + - ]: 56832 : m_dsupint[0].push_back( ded[3] );
748 [ + - ]: 56832 : m_dsupint[0].push_back( ded[4] );
749 [ + - ]: 56832 : m_domedgeint.erase( d[ed] );
750 : : }
751 : : }
752 : : }
753 : :
754 [ + + ]: 103954 : for (const auto& N : untri) {
755 : 103589 : int f = 0;
756 : : tk::real sig[3];
757 [ + - ][ + + ]: 414356 : decltype(m_domedgeint)::const_iterator d[3];
758 [ + + ]: 167483 : for (const auto& [p,q] : tk::lpoet) {
759 : 158289 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
760 [ + + ]: 158289 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
761 [ + - ]: 158289 : d[f] = m_domedgeint.find( ed );
762 [ + + ]: 158289 : if (d[f] == end(m_domedgeint)) break; else ++f;
763 : : }
764 [ + + ]: 103589 : if (f == 3) {
765 [ + - ]: 9194 : m_dsupedge[1].push_back( N[0] );
766 [ + - ]: 9194 : m_dsupedge[1].push_back( N[1] );
767 [ + - ]: 9194 : m_dsupedge[1].push_back( N[2] );
768 [ + + ]: 36776 : for (int ed=0; ed<3; ++ed) {
769 : 27582 : const auto& ded = d[ed]->second;
770 [ + - ]: 27582 : m_dsupint[1].push_back( sig[ed] * ded[0] );
771 [ + - ]: 27582 : m_dsupint[1].push_back( sig[ed] * ded[1] );
772 [ + - ]: 27582 : m_dsupint[1].push_back( sig[ed] * ded[2] );
773 [ + - ]: 27582 : m_dsupint[1].push_back( ded[3] );
774 [ + - ]: 27582 : m_dsupint[1].push_back( ded[4] );
775 [ + - ]: 27582 : m_domedgeint.erase( d[ed] );
776 : : }
777 : : }
778 : : }
779 : :
780 [ + - ]: 365 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
781 [ + - ]: 365 : m_dsupint[2].resize( m_domedgeint.size()*5 );
782 : 365 : std::size_t k = 0;
783 [ + + ]: 23440 : for (const auto& [ed,d] : m_domedgeint) {
784 : 23075 : auto e = m_dsupedge[2].data() + k*2;
785 [ + - ]: 23075 : e[0] = tk::cref_find( lid, ed[0] );
786 [ + - ]: 23075 : e[1] = tk::cref_find( lid, ed[1] );
787 : 23075 : auto i = m_dsupint[2].data() + k*5;
788 : 23075 : i[0] = d[0];
789 : 23075 : i[1] = d[1];
790 : 23075 : i[2] = d[2];
791 : 23075 : i[3] = d[3];
792 : 23075 : i[4] = d[4];
793 : 23075 : ++k;
794 : : }
795 : :
796 [ + + ]: 365 : if (g_cfg.get< tag::fct >()) {
797 : 324 : const auto ncomp = m_u.nprop();
798 [ + - ]: 324 : m_dsuplim[0].resize( m_dsupedge[0].size() * 6 * ncomp );
799 [ + - ]: 324 : m_dsuplim[1].resize( m_dsupedge[1].size() * 3 * ncomp );
800 [ + - ]: 324 : m_dsuplim[2].resize( m_dsupedge[2].size() * ncomp );
801 : : }
802 : :
803 : 365 : tk::destroy( m_domedgeint );
804 : :
805 : : //std::cout << std::setprecision(2)
806 : : // << "superedges: ntet:" << m_dsupedge[0].size()/4 << "(nedge:"
807 : : // << m_dsupedge[0].size()/4*6 << ","
808 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[0].size()/4*6 ) /
809 : : // static_cast< tk::real >( nedge )
810 : : // << "%) + ntri:" << m_dsupedge[1].size()/3
811 : : // << "(nedge:" << m_dsupedge[1].size() << ","
812 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[1].size() ) /
813 : : // static_cast< tk::real >( nedge )
814 : : // << "%) + nedge:"
815 : : // << m_dsupedge[2].size()/2 << "("
816 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[2].size()/2 ) /
817 : : // static_cast< tk::real >( nedge )
818 : : // << "%) = " << m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
819 : : // m_dsupedge[2].size()/2 << " of "<< nedge << " total edges\n";
820 : :
821 [ - + ][ - - ]: 365 : Assert( m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
[ - - ][ - - ]
822 : : m_dsupedge[2].size()/2 == nedge,
823 : : "Not all edges accounted for in superedge groups" );
824 : 365 : }
825 : :
826 : : void
827 : : // cppcheck-suppress unusedFunction
828 : 365 : ChoCG::merge()
829 : : // *****************************************************************************
830 : : // Combine own and communicated portions of the integrals
831 : : // *****************************************************************************
832 : : {
833 : : // Combine own and communicated contributions to boundary point normals
834 : 365 : bnorm();
835 : :
836 : : // Convert integrals into streamable data structures
837 : 365 : streamable();
838 : :
839 : : // Enforce boundary conditions on initial conditions
840 : 365 : BC( m_u, Disc()->T() );
841 : :
842 : : // Start measuring initial div-free time
843 : 365 : m_timer.emplace_back();
844 : :
845 : : // Compute initial momentum flux
846 [ + - ][ + - ]: 365 : thisProxy[ thisIndex ].wait4div();
847 [ + - ][ + - ]: 365 : thisProxy[ thisIndex ].wait4sgrad();
848 : 365 : div( m_u );
849 : 365 : }
850 : :
851 : : void
852 : 11704 : ChoCG::fingrad( tk::Fields& grad,
853 : : std::unordered_map< std::size_t, std::vector< tk::real > >& gradc )
854 : : // *****************************************************************************
855 : : // Finalize computing gradient
856 : : //! \param[in,out] grad Gradient to finalize
857 : : //! \param[in,out] gradc Gradient communication buffer to finalize
858 : : // *****************************************************************************
859 : : {
860 [ + - ]: 11704 : auto d = Disc();
861 [ + - ]: 11704 : const auto lid = d->Lid();
862 : :
863 : : // Combine own and communicated contributions
864 [ + + ]: 231417 : for (const auto& [g,r] : gradc) {
865 [ + - ]: 219713 : auto i = tk::cref_find( lid, g );
866 [ + - ][ + + ]: 1244834 : for (std::size_t c=0; c<r.size(); ++c) grad(i,c) += r[c];
867 : : }
868 : 11704 : tk::destroy(gradc);
869 : :
870 : : // Divide weak result by nodal volume
871 : 11704 : const auto& vol = d->Vol();
872 [ + + ]: 1023026 : for (std::size_t p=0; p<grad.nunk(); ++p) {
873 [ + + ]: 5421088 : for (std::size_t c=0; c<grad.nprop(); ++c) {
874 [ + - ]: 4409766 : grad(p,c) /= vol[p];
875 : : }
876 : : }
877 : 11704 : }
878 : :
879 : : void
880 : 4328 : ChoCG::div( const tk::Fields& u )
881 : : // *****************************************************************************
882 : : // Start computing divergence
883 : : // \para[in] u Vector field whose divergence to compute
884 : : // *****************************************************************************
885 : : {
886 [ + - ]: 4328 : auto d = Disc();
887 [ + - ]: 4328 : const auto lid = d->Lid();
888 : :
889 : : // Finalize momentum flux communications if needed
890 [ + + ]: 4328 : if (m_np == 1) {
891 [ + - ]: 333 : fingrad( m_flux, m_fluxc );
892 [ + - ]: 333 : physics::symbc( m_flux, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
893 : : }
894 : :
895 : : // Compute divergence
896 [ + - ]: 4328 : std::fill( begin(m_div), end(m_div), 0.0 );
897 [ + - ]: 4328 : chorin::div( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel,
898 : 4328 : d->Dt(), m_pr, m_pgrad, u, m_div, m_np>1 );
899 : :
900 : : // Communicate velocity divergence to other chares on chare-boundary
901 [ + + ]: 4328 : if (d->NodeCommMap().empty()) {
902 [ + - ]: 193 : comdiv_complete();
903 : : } else {
904 [ + + ]: 47347 : for (const auto& [c,n] : d->NodeCommMap()) {
905 : 43212 : decltype(m_divc) exp;
906 [ + - ][ + - ]: 221816 : for (auto g : n) exp[g] = m_div[ tk::cref_find(lid,g) ];
[ + + ]
907 [ + - ][ + - ]: 43212 : thisProxy[c].comdiv( exp );
908 : 43212 : }
909 : : }
910 [ + - ]: 4328 : owndiv_complete();
911 : 4328 : }
912 : :
913 : : void
914 : 43212 : ChoCG::comdiv( const std::unordered_map< std::size_t, tk::real >& indiv )
915 : : // *****************************************************************************
916 : : // Receive contributions to velocity divergence on chare-boundaries
917 : : //! \param[in] indiv Partial contributions of velocity divergence to
918 : : //! chare-boundary nodes. Key: global mesh node IDs, value: contribution.
919 : : //! \details This function receives contributions to m_div, which stores the
920 : : //! velocity divergence at mesh nodes. While m_div stores own contributions,
921 : : //! m_divc collects the neighbor chare contributions during communication.
922 : : //! This way work on m_div and m_divc is overlapped.
923 : : // *****************************************************************************
924 : : {
925 [ + - ][ + + ]: 221816 : for (const auto& [g,r] : indiv) m_divc[g] += r;
926 : :
927 : : // When we have heard from all chares we communicate with, this chare is done
928 [ + + ]: 43212 : if (++m_ndiv == Disc()->NodeCommMap().size()) {
929 : 4135 : m_ndiv = 0;
930 : 4135 : comdiv_complete();
931 : : }
932 : 43212 : }
933 : :
934 : : void
935 : 3413 : ChoCG::velgrad()
936 : : // *****************************************************************************
937 : : // Start computing velocity gradient
938 : : // *****************************************************************************
939 : : {
940 : 3413 : auto d = Disc();
941 : :
942 : : // Compute momentum flux
943 : 3413 : m_vgrad.fill( 0.0 );
944 : 3413 : chorin::vgrad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad );
945 : :
946 : : // Communicate velocity divergence to other chares on chare-boundary
947 : 3413 : const auto& lid = d->Lid();
948 [ + + ]: 3413 : if (d->NodeCommMap().empty()) {
949 : 151 : comvgrad_complete();
950 : : } else {
951 [ + + ]: 41814 : for (const auto& [c,n] : d->NodeCommMap()) {
952 : 38552 : decltype(m_vgradc) exp;
953 [ + - ][ + - ]: 185790 : for (auto g : n) exp[g] = m_vgrad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
954 [ + - ][ + - ]: 38552 : thisProxy[c].comvgrad( exp );
955 : 38552 : }
956 : : }
957 : 3413 : ownvgrad_complete();
958 : 3413 : }
959 : :
960 : : void
961 : 38552 : ChoCG::comvgrad(
962 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
963 : : // *****************************************************************************
964 : : // Receive contributions to velocity gradients on chare-boundaries
965 : : //! \param[in] ingrad Partial contributions of momentum flux to
966 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
967 : : //! \details This function receives contributions to m_vgrad, which stores the
968 : : //! velocity gradients at mesh nodes. While m_vgrad stores own contributions,
969 : : //! m_vgradc collects the neighbor chare contributions during communication.
970 : : //! This way work on m_vgrad and m_vgradc is overlapped.
971 : : // *****************************************************************************
972 : : {
973 : : using tk::operator+=;
974 [ + - ][ + - ]: 185790 : for (const auto& [g,r] : ingrad) m_vgradc[g] += r;
[ + + ]
975 : :
976 : : // When we have heard from all chares we communicate with, this chare is done
977 [ + + ]: 38552 : if (++m_nvgrad == Disc()->NodeCommMap().size()) {
978 : 3262 : m_nvgrad = 0;
979 : 3262 : comvgrad_complete();
980 : : }
981 : 38552 : }
982 : :
983 : : void
984 : 333 : ChoCG::flux()
985 : : // *****************************************************************************
986 : : // Start computing momentum flux
987 : : // *****************************************************************************
988 : : {
989 : 333 : auto d = Disc();
990 : :
991 : : // Finalize computing velocity gradients
992 : 333 : fingrad( m_vgrad, m_vgradc );
993 : :
994 : : // Compute momentum flux
995 : 333 : m_flux.fill( 0.0 );
996 : 333 : chorin::flux( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad,
997 : 333 : m_flux );
998 : :
999 : : // Communicate velocity divergence to other chares on chare-boundary
1000 : 333 : const auto& lid = d->Lid();
1001 [ + + ]: 333 : if (d->NodeCommMap().empty()) {
1002 : 11 : comflux_complete();
1003 : : } else {
1004 [ + + ]: 3874 : for (const auto& [c,n] : d->NodeCommMap()) {
1005 : 3552 : decltype(m_fluxc) exp;
1006 [ + - ][ + - ]: 17930 : for (auto g : n) exp[g] = m_flux[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1007 [ + - ][ + - ]: 3552 : thisProxy[c].comflux( exp );
1008 : 3552 : }
1009 : : }
1010 : 333 : ownflux_complete();
1011 : 333 : }
1012 : :
1013 : : void
1014 : 3552 : ChoCG::comflux(
1015 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& influx )
1016 : : // *****************************************************************************
1017 : : // Receive contributions to momentum flux on chare-boundaries
1018 : : //! \param[in] influx Partial contributions of momentum flux to
1019 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
1020 : : //! \details This function receives contributions to m_flux, which stores the
1021 : : //! momentum flux at mesh nodes. While m_flux stores own contributions,
1022 : : //! m_fluxc collects the neighbor chare contributions during communication.
1023 : : //! This way work on m_flux and m_fluxc is overlapped.
1024 : : // *****************************************************************************
1025 : : {
1026 : : using tk::operator+=;
1027 [ + - ][ + - ]: 17930 : for (const auto& [g,r] : influx) m_fluxc[g] += r;
[ + + ]
1028 : :
1029 : : // When we have heard from all chares we communicate with, this chare is done
1030 [ + + ]: 3552 : if (++m_nflux == Disc()->NodeCommMap().size()) {
1031 : 322 : m_nflux = 0;
1032 : 322 : comflux_complete();
1033 : : }
1034 : 3552 : }
1035 : :
1036 : : void
1037 : 4328 : ChoCG::pinit()
1038 : : // *****************************************************************************
1039 : : // Initialize Poisson solve
1040 : : // *****************************************************************************
1041 : : {
1042 [ + - ]: 4328 : auto d = Disc();
1043 [ + - ]: 4328 : const auto lid = d->Lid();
1044 : 4328 : const auto& coord = d->Coord();
1045 : 4328 : const auto& x = coord[0];
1046 : 4328 : const auto& y = coord[1];
1047 : 4328 : const auto& z = coord[2];
1048 : :
1049 : : // Combine own and communicated contributions to velocity divergence
1050 [ + - ][ + + ]: 87327 : for (const auto& [g,r] : m_divc) m_div[ tk::cref_find(lid,g) ] += r;
1051 : 4328 : tk::destroy(m_divc);
1052 : :
1053 [ + + ][ + + ]: 359528 : if (m_np > 1) for (auto& div : m_div) div /= d->Dt();
1054 : :
1055 : : // Configure Dirichlet BCs
1056 : : std::unordered_map< std::size_t,
1057 : 4328 : std::vector< std::pair< int, tk::real > > > dirbc;
1058 [ + + ]: 4328 : if (!g_cfg.get< tag::pre_bc_dir >().empty()) {
1059 [ + - ]: 752 : auto ic = problems::PRESSURE_IC();
1060 : 752 : std::size_t nmask = 1 + 1;
1061 [ - + ][ - - ]: 752 : Assert( m_dirbcmaskp.size() % nmask == 0, "Size mismatch" );
[ - - ][ - - ]
1062 [ + + ]: 8816 : for (std::size_t i=0; i<m_dirbcmaskp.size()/nmask; ++i) {
1063 : 8064 : auto p = m_dirbcmaskp[i*nmask+0]; // local node id
1064 : 8064 : auto mask = m_dirbcmaskp[i*nmask+1];
1065 [ + + ]: 8064 : if (mask == 1) { // mask == 1: IC value
1066 [ + + ][ + - ]: 2784 : auto val = m_np>1 ? 0.0 : ic( x[p], y[p], z[p] );
1067 [ + - ][ + - ]: 2784 : dirbc[p] = {{ { 1, val } }};
1068 [ + - ][ + - ]: 5280 : } else if (mask == 2 && !m_dirbcvalp.empty()) { // mask == 2: BC value
[ + - ]
1069 [ + + ]: 5280 : auto val = m_np>1 ? 0.0 : m_dirbcvalp[i*nmask+1];
1070 [ + - ][ + - ]: 5280 : dirbc[p] = {{ { 1, val } }};
1071 : : }
1072 : : }
1073 : 752 : }
1074 : :
1075 : : // Configure Neumann BCs
1076 : 4328 : std::vector< tk::real > neubc;
1077 [ + - ]: 4328 : auto pg = problems::PRESSURE_GRAD();
1078 [ + + ]: 4328 : if (pg) {
1079 : : // Collect Neumann BC elements
1080 [ + - ]: 2 : std::vector< std::uint8_t > besym( m_triinpoel.size(), 0 );
1081 [ + + ]: 8 : for (auto s : g_cfg.get< tag::pre_bc_sym >()) {
1082 [ + - ]: 6 : auto k = m_bface.find(s);
1083 [ + + ][ + + ]: 210 : if (k != end(m_bface)) for (auto f : k->second) besym[f] = 1;
1084 : : }
1085 : : // Setup Neumann BCs
1086 [ + - ]: 2 : neubc.resize( x.size(), 0.0 );
1087 [ + + ]: 410 : for (std::size_t e=0; e<m_triinpoel.size()/3; ++e) {
1088 [ + + ]: 408 : if (besym[e]) {
1089 : 204 : const auto N = m_triinpoel.data() + e*3;
1090 : : tk::real n[3];
1091 : 204 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
1092 : 204 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
1093 : : n[0], n[1], n[2] );
1094 [ + - ]: 204 : auto g = pg( x[N[0]], y[N[0]], z[N[0]] );
1095 : 204 : neubc[ N[0] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1096 [ + - ]: 204 : g = pg( x[N[1]], y[N[1]], z[N[1]] );
1097 : 204 : neubc[ N[1] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1098 [ + - ]: 204 : g = pg( x[N[2]], y[N[2]], z[N[2]] );
1099 : 204 : neubc[ N[2] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1100 : : }
1101 : : }
1102 : 2 : }
1103 : :
1104 : : // Set hydrostat
1105 : 4328 : auto h = g_cfg.get< tag::pre_hydrostat >();
1106 [ + + ]: 4328 : if (h != std::numeric_limits< uint64_t >::max()) {
1107 [ + - ]: 72 : auto pi = lid.find( h );
1108 [ + + ]: 72 : if (pi != end(lid)) {
1109 : 36 : auto p = pi->second;
1110 [ + - ]: 36 : auto ic = problems::PRESSURE_IC();
1111 [ + + ][ + - ]: 36 : auto val = m_np>1 ? 0.0 : ic( x[p], y[p], z[p] );
1112 [ + - ]: 36 : auto& b = dirbc[p];
1113 [ + - ][ + - ]: 36 : if (b.empty()) b = {{ { 1, val }} };
1114 : 36 : }
1115 : : }
1116 : :
1117 : : // Configure right hand side
1118 [ + - ]: 4328 : auto pr = problems::PRESSURE_RHS();
1119 [ + + ]: 4328 : if (pr) {
1120 : 32 : const auto& vol = d->Vol();
1121 [ + + ]: 2174 : for (std::size_t i=0; i<x.size(); ++i) {
1122 [ + - ]: 2142 : m_div[i] = pr( x[i], y[i], z[i] ) * vol[i];
1123 : : }
1124 : : }
1125 : :
1126 : : // Initialize Poisson solve
1127 : 4328 : const auto& pc = g_cfg.get< tag::pre_pc >();
1128 [ + - ][ + - ]: 8656 : m_cgpre[ thisIndex ].ckLocal()->init( {}, m_div, neubc, dirbc, pc,
[ + - ]
1129 [ + - ][ + - ]: 8656 : CkCallback( CkIndex_ChoCG::psolve(), thisProxy[thisIndex] ) );
[ + - ]
1130 : 4328 : }
1131 : :
1132 : : void
1133 : 4328 : ChoCG::psolve()
1134 : : // *****************************************************************************
1135 : : // Solve Poisson equation
1136 : : // *****************************************************************************
1137 : : {
1138 : 4328 : auto iter = g_cfg.get< tag::pre_iter >();
1139 : 4328 : auto tol = g_cfg.get< tag::pre_tol >();
1140 : 4328 : auto verbose = g_cfg.get< tag::pre_verbose >();
1141 : :
1142 : 4328 : auto c = m_np != 1 ?
1143 : 3995 : CkCallback( CkIndex_ChoCG::sgrad(), thisProxy[thisIndex] ) :
1144 [ + + ][ + - ]: 4328 : CkCallback( CkIndex_ChoCG::psolved(), thisProxy[thisIndex] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ][ - - ]
[ - - ]
1145 : :
1146 [ + - ][ + - ]: 4328 : m_cgpre[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose, c );
[ + - ]
1147 : 4328 : }
1148 : :
1149 : : void
1150 : 3995 : ChoCG::sgrad()
1151 : : // *****************************************************************************
1152 : : // Compute recent conjugate gradients solution gradient
1153 : : // *****************************************************************************
1154 : : {
1155 [ + - ]: 3995 : auto d = Disc();
1156 : :
1157 [ + - ][ + - ]: 3995 : auto sol = m_cgpre[ thisIndex ].ckLocal()->solution();
[ + - ]
1158 [ + - ]: 3995 : m_sgrad.fill( 0.0 );
1159 [ + - ]: 3995 : chorin::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, sol, m_sgrad );
1160 : :
1161 : : // Send gradient contributions to neighbor chares
1162 [ + + ]: 3995 : if (d->NodeCommMap().empty()) {
1163 [ + - ]: 182 : comsgrad_complete();
1164 : : } else {
1165 : 3813 : const auto& lid = d->Lid();
1166 [ + + ]: 43473 : for (const auto& [c,n] : d->NodeCommMap()) {
1167 : 39660 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1168 [ + - ][ + - ]: 203886 : for (auto g : n) exp[g] = m_sgrad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1169 [ + - ][ + - ]: 39660 : thisProxy[c].comsgrad( exp );
1170 : 39660 : }
1171 : : }
1172 [ + - ]: 3995 : ownsgrad_complete();
1173 : 3995 : }
1174 : :
1175 : : void
1176 : 39660 : ChoCG::comsgrad(
1177 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
1178 : : // *****************************************************************************
1179 : : // Receive contributions to conjugrate gradients solution gradient
1180 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
1181 : : //! global mesh node IDs, value: contributions for all scalar components.
1182 : : //! \details This function receives contributions to m_sgrad, which stores the
1183 : : //! gradients at mesh nodes. While m_sgrad stores own contributions, m_sgradc
1184 : : //! collects the neighbor chare contributions during communication. This way
1185 : : //! work on m_sgrad and m_sgradc is overlapped.
1186 : : // *****************************************************************************
1187 : : {
1188 : : using tk::operator+=;
1189 [ + - ][ + - ]: 203886 : for (const auto& [g,r] : ingrad) m_sgradc[g] += r;
[ + + ]
1190 : :
1191 : : // When we have heard from all chares we communicate with, this chare is done
1192 [ + + ]: 39660 : if (++m_nsgrad == Disc()->NodeCommMap().size()) {
1193 : 3813 : m_nsgrad = 0;
1194 : 3813 : comsgrad_complete();
1195 : : }
1196 : 39660 : }
1197 : :
1198 : : void
1199 : 4328 : ChoCG::psolved()
1200 : : // *****************************************************************************
1201 : : // Continue setup after Poisson solve and gradient computation
1202 : : // *****************************************************************************
1203 : : {
1204 : 4328 : auto d = Disc();
1205 : :
1206 [ + + ][ + - ]: 4328 : if (thisIndex == 0) d->pit( m_cgpre[ thisIndex ].ckLocal()->it() );
[ + - ][ + - ]
1207 : :
1208 [ + + ]: 4328 : if (m_np != 1) {
1209 : : // Finalize gradient communications
1210 : 3995 : fingrad( m_sgrad, m_sgradc );
1211 : : // Project velocity to divergence-free subspace
1212 [ + + ]: 3995 : auto dt = m_np > 1 ? d->Dt() : 1.0;
1213 [ + + ]: 384497 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1214 : 380502 : m_u(i,0) -= dt * m_sgrad(i,0);
1215 : 380502 : m_u(i,1) -= dt * m_sgrad(i,1);
1216 : 380502 : m_u(i,2) -= dt * m_sgrad(i,2);
1217 : : }
1218 : : // Enforce boundary conditions
1219 : 3995 : BC( m_u, d->T() + d->Dt() );
1220 : : }
1221 : :
1222 [ + + ]: 4328 : if (d->Initial()) {
1223 : :
1224 [ + + ]: 698 : if (g_cfg.get< tag::nstep >() == 1) { // test first Poisson solve only
1225 : :
1226 [ + - ][ + - ]: 32 : m_pr = m_cgpre[ thisIndex ].ckLocal()->solution();
[ + - ]
1227 [ + - ][ + - ]: 32 : thisProxy[ thisIndex ].wait4step();
1228 [ + - ][ + - ]: 32 : writeFields( CkCallback(CkIndex_ChoCG::diag(), thisProxy[thisIndex]) );
[ + - ][ + - ]
1229 : :
1230 : : } else {
1231 : :
1232 [ + + ]: 666 : if (++m_np < 2) {
1233 : : // Compute momentum flux for initial pressure-Poisson rhs
1234 [ + - ][ + - ]: 333 : thisProxy[ thisIndex ].wait4vgrad();
1235 [ + - ][ + - ]: 333 : thisProxy[ thisIndex ].wait4flux();
1236 [ + - ][ + - ]: 333 : thisProxy[ thisIndex ].wait4div();
1237 : 333 : velgrad();
1238 : : } else {
1239 [ + + ]: 333 : if (thisIndex == 0) {
1240 [ + - ]: 99 : tk::Print() << "Initial div-free time: " << m_timer[0].dsec()
1241 [ + - ][ + - ]: 66 : << " sec\n";
[ + - ]
1242 : : }
1243 : : // Assign initial pressure and compute its gradient
1244 [ + - ][ + - ]: 333 : m_pr = m_cgpre[ thisIndex ].ckLocal()->solution();
[ + - ]
1245 : 333 : pgrad();
1246 : : }
1247 : :
1248 : : }
1249 : :
1250 : : } else {
1251 : :
1252 : : // Update pressure and compute its gradient
1253 : : using tk::operator+=;
1254 [ + - ][ + - ]: 3630 : m_pr += m_cgpre[ thisIndex ].ckLocal()->solution();
[ + - ]
1255 : 3630 : pgrad();
1256 : :
1257 : : }
1258 : 4328 : }
1259 : :
1260 : : void
1261 : 3963 : ChoCG::pgrad()
1262 : : // *****************************************************************************
1263 : : // Compute pressure gradient
1264 : : // *****************************************************************************
1265 : : {
1266 : 3963 : auto d = Disc();
1267 : :
1268 [ + - ][ + - ]: 3963 : thisProxy[ thisIndex ].wait4pgrad();
1269 : :
1270 : 3963 : m_pgrad.fill( 0.0 );
1271 : 3963 : chorin::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_pr, m_pgrad );
1272 : :
1273 : : // Send gradient contributions to neighbor chares
1274 [ + + ]: 3963 : if (d->NodeCommMap().empty()) {
1275 : 181 : compgrad_complete();
1276 : : } else {
1277 : 3782 : const auto& lid = d->Lid();
1278 [ + + ]: 43254 : for (const auto& [c,n] : d->NodeCommMap()) {
1279 : 39472 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1280 [ + - ][ + - ]: 202430 : for (auto g : n) exp[g] = m_pgrad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1281 [ + - ][ + - ]: 39472 : thisProxy[c].compgrad( exp );
1282 : 39472 : }
1283 : : }
1284 : 3963 : ownpgrad_complete();
1285 : 3963 : }
1286 : :
1287 : : void
1288 : 39472 : ChoCG::compgrad(
1289 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
1290 : : // *****************************************************************************
1291 : : // Receive contributions to pressure gradient on chare-boundaries
1292 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
1293 : : //! global mesh node IDs, value: contributions for all scalar components.
1294 : : //! \details This function receives contributions to m_pgrad, which stores the
1295 : : //! gradients at mesh nodes. While m_pgrad stores own contributions, m_pgradc
1296 : : //! collects the neighbor chare contributions during communication. This way
1297 : : //! work on m_pgrad and m_pgradc is overlapped.
1298 : : // *****************************************************************************
1299 : : {
1300 : : using tk::operator+=;
1301 [ + - ][ + - ]: 202430 : for (const auto& [g,r] : ingrad) m_pgradc[g] += r;
[ + + ]
1302 : :
1303 : : // When we have heard from all chares we communicate with, this chare is done
1304 [ + + ]: 39472 : if (++m_npgrad == Disc()->NodeCommMap().size()) {
1305 : 3782 : m_npgrad = 0;
1306 : 3782 : compgrad_complete();
1307 : : }
1308 : 39472 : }
1309 : :
1310 : : void
1311 : 3963 : ChoCG::finpgrad()
1312 : : // *****************************************************************************
1313 : : // Compute pressure gradient
1314 : : // *****************************************************************************
1315 : : {
1316 : 3963 : auto d = Disc();
1317 : :
1318 : : // Finalize pressure gradient communications
1319 : 3963 : fingrad( m_pgrad, m_pgradc );
1320 : :
1321 [ + + ]: 3963 : if (d->Initial()) {
1322 [ + - ][ + - ]: 333 : writeFields( CkCallback(CkIndex_ChoCG::start(), thisProxy[thisIndex]) );
[ + - ][ + - ]
1323 : : } else {
1324 : 3630 : diag();
1325 : : }
1326 : 3963 : }
1327 : :
1328 : : void
1329 : 333 : ChoCG::start()
1330 : : // *****************************************************************************
1331 : : // Start time stepping
1332 : : // *****************************************************************************
1333 : : {
1334 : : // Set flag that indicates that we are now during time stepping
1335 : 333 : Disc()->Initial( 0 );
1336 : : // Start timer measuring time stepping wall clock time
1337 : 333 : Disc()->Timer().zero();
1338 : : // Zero grind-timer
1339 : 333 : Disc()->grindZero();
1340 : : // Continue to first time step
1341 : 333 : dt();
1342 : 333 : }
1343 : :
1344 : : void
1345 : : // cppcheck-suppress unusedFunction
1346 : 2920 : ChoCG::aec()
1347 : : // *****************************************************************************
1348 : : // Compute antidiffusive contributions: P+/-
1349 : : // *****************************************************************************
1350 : : {
1351 : 2920 : auto d = Disc();
1352 : 2920 : const auto ncomp = m_u.nprop();
1353 : 2920 : const auto& lid = d->Lid();
1354 : :
1355 : : // Antidiffusive contributions: P+/-
1356 : :
1357 : 2920 : auto ctau = g_cfg.get< tag::fctdif >();
1358 : 2920 : m_p.fill( 0.0 );
1359 : :
1360 : : // tetrahedron superedges
1361 [ + + ]: 25760 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1362 : 22840 : const auto N = m_dsupedge[0].data() + e*4;
1363 : 22840 : const auto D = m_dsupint[0].data();
1364 : 22840 : std::size_t i = 0;
1365 [ + + ]: 159880 : for (const auto& [p,q] : tk::lpoed) {
1366 : 137040 : auto dif = D[(e*6+i)*5+3];
1367 [ + + ]: 548160 : for (std::size_t c=0; c<ncomp; ++c) {
1368 [ + - ][ + - ]: 411120 : auto aec = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1369 : 411120 : auto a = c*2;
1370 : 411120 : auto b = a+1;
1371 [ - + ]: 411120 : if (aec > 0.0) std::swap(a,b);
1372 [ + - ]: 411120 : m_p(N[p],a) -= aec;
1373 [ + - ]: 411120 : m_p(N[q],b) += aec;
1374 : : }
1375 : 137040 : ++i;
1376 : : }
1377 : : }
1378 : :
1379 : : // triangle superedges
1380 [ + + ]: 25460 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1381 : 22540 : const auto N = m_dsupedge[1].data() + e*3;
1382 : 22540 : const auto D = m_dsupint[1].data();
1383 : 22540 : std::size_t i = 0;
1384 [ + + ]: 90160 : for (const auto& [p,q] : tk::lpoet) {
1385 : 67620 : auto dif = D[(e*3+i)*5+3];
1386 [ + + ]: 270480 : for (std::size_t c=0; c<ncomp; ++c) {
1387 [ + - ][ + - ]: 202860 : auto aec = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1388 : 202860 : auto a = c*2;
1389 : 202860 : auto b = a+1;
1390 [ - + ]: 202860 : if (aec > 0.0) std::swap(a,b);
1391 [ + - ]: 202860 : m_p(N[p],a) -= aec;
1392 [ + - ]: 202860 : m_p(N[q],b) += aec;
1393 : : }
1394 : 67620 : ++i;
1395 : : }
1396 : : }
1397 : :
1398 : : // edges
1399 [ + + ]: 74000 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1400 : 71080 : const auto N = m_dsupedge[2].data() + e*2;
1401 : 71080 : const auto dif = m_dsupint[2][e*5+3];
1402 [ + + ]: 284320 : for (std::size_t c=0; c<ncomp; ++c) {
1403 [ + - ][ + - ]: 213240 : auto aec = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1404 : 213240 : auto a = c*2;
1405 : 213240 : auto b = a+1;
1406 [ - + ]: 213240 : if (aec > 0.0) std::swap(a,b);
1407 [ + - ]: 213240 : m_p(N[0],a) -= aec;
1408 [ + - ]: 213240 : m_p(N[1],b) += aec;
1409 : : }
1410 : : }
1411 : :
1412 : : // Apply symmetry BCs on AEC
1413 [ + + ]: 15100 : for (std::size_t i=0; i<m_symbcnodes.size(); ++i) {
1414 : 12180 : auto p = m_symbcnodes[i];
1415 : 12180 : auto n = m_symbcnorms.data() + i*3;
1416 : 12180 : auto rvnp = m_p(p,0)*n[0] + m_p(p,2)*n[1] + m_p(p,4)*n[2];
1417 : 12180 : auto rvnn = m_p(p,1)*n[0] + m_p(p,3)*n[1] + m_p(p,5)*n[2];
1418 : 12180 : m_p(p,0) -= rvnp * n[0];
1419 : 12180 : m_p(p,1) -= rvnn * n[0];
1420 : 12180 : m_p(p,2) -= rvnp * n[1];
1421 : 12180 : m_p(p,3) -= rvnn * n[1];
1422 : 12180 : m_p(p,4) -= rvnp * n[2];
1423 : 12180 : m_p(p,5) -= rvnn * n[2];
1424 : : }
1425 : :
1426 : : // Communicate antidiffusive edge and low-order solution contributions
1427 [ + + ]: 2920 : if (d->NodeCommMap().empty()) {
1428 : 20 : comaec_complete();
1429 : : } else {
1430 [ + + ]: 37780 : for (const auto& [c,n] : d->NodeCommMap()) {
1431 : 34880 : decltype(m_pc) exp;
1432 [ + - ][ + - ]: 165880 : for (auto g : n) exp[g] = m_p[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1433 [ + - ][ + - ]: 34880 : thisProxy[c].comaec( exp );
1434 : 34880 : }
1435 : : }
1436 : 2920 : ownaec_complete();
1437 : 2920 : }
1438 : :
1439 : : void
1440 : 34880 : ChoCG::comaec( const std::unordered_map< std::size_t,
1441 : : std::vector< tk::real > >& inaec )
1442 : : // *****************************************************************************
1443 : : // Receive antidiffusive and low-order contributions on chare-boundaries
1444 : : //! \param[in] inaec Partial contributions of antidiffusive edge and low-order
1445 : : //! solution contributions on chare-boundary nodes. Key: global mesh node IDs,
1446 : : //! value: 0: antidiffusive contributions, 1: low-order solution.
1447 : : // *****************************************************************************
1448 : : {
1449 : : using tk::operator+=;
1450 [ + - ][ + - ]: 165880 : for (const auto& [g,a] : inaec) m_pc[g] += a;
[ + + ]
1451 : :
1452 : : // When we have heard from all chares we communicate with, this chare is done
1453 [ + + ]: 34880 : if (++m_naec == Disc()->NodeCommMap().size()) {
1454 : 2900 : m_naec = 0;
1455 : 2900 : comaec_complete();
1456 : : }
1457 : 34880 : }
1458 : :
1459 : : void
1460 : 2920 : ChoCG::alw()
1461 : : // *****************************************************************************
1462 : : // Compute allowed limits, Q+/-
1463 : : // *****************************************************************************
1464 : : {
1465 : 2920 : auto d = Disc();
1466 : 2920 : const auto npoin = m_u.nunk();
1467 : 2920 : const auto ncomp = m_u.nprop();
1468 : 2920 : const auto& lid = d->Lid();
1469 : 2920 : const auto& vol = d->Vol();
1470 : :
1471 : : // Combine own and communicated contributions to antidiffusive contributions
1472 : : // and low-order solution
1473 [ + + ]: 56000 : for (const auto& [g,p] : m_pc) {
1474 [ + - ]: 53080 : auto i = tk::cref_find( lid, g );
1475 [ + - ][ + + ]: 371560 : for (std::size_t c=0; c<p.size(); ++c) m_p(i,c) += p[c];
1476 : : }
1477 : 2920 : tk::destroy(m_pc);
1478 : :
1479 : : // Finish computing antidiffusive contributions and low-order solution
1480 : 2920 : auto dt = m_rkcoef[m_stage] * d->Dt();
1481 [ + + ]: 79900 : for (std::size_t i=0; i<npoin; ++i) {
1482 [ + + ]: 307920 : for (std::size_t c=0; c<ncomp; ++c) {
1483 : 230940 : auto a = c*2;
1484 : 230940 : auto b = a+1;
1485 : 230940 : m_p(i,a) /= vol[i];
1486 : 230940 : m_p(i,b) /= vol[i];
1487 : : // low-order solution
1488 : 230940 : m_rhs(i,c) = m_u(i,c) - dt*m_rhs(i,c)/vol[i] - m_p(i,a) - m_p(i,b);
1489 : : }
1490 : : }
1491 : :
1492 : : // Allowed limits: Q+/-
1493 : :
1494 : : using std::max;
1495 : : using std::min;
1496 : :
1497 : 2920 : auto large = std::numeric_limits< tk::real >::max();
1498 [ + + ]: 79900 : for (std::size_t i=0; i<m_q.nunk(); ++i) {
1499 [ + + ]: 307920 : for (std::size_t c=0; c<m_q.nprop()/2; ++c) {
1500 : 230940 : m_q(i,c*2+0) = -large;
1501 : 230940 : m_q(i,c*2+1) = +large;
1502 : : }
1503 : : }
1504 : :
1505 : : // tetrahedron superedges
1506 [ + + ]: 25760 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1507 : 22840 : const auto N = m_dsupedge[0].data() + e*4;
1508 [ + + ]: 91360 : for (std::size_t c=0; c<ncomp; ++c) {
1509 : 68520 : auto a = c*2;
1510 : 68520 : auto b = a+1;
1511 [ + + ]: 479640 : for (const auto& [p,q] : tk::lpoed) {
1512 : : tk::real alwp, alwn;
1513 [ - + ]: 411120 : if (g_cfg.get< tag::fctclip >()) {
1514 [ - - ][ - - ]: 0 : alwp = max( m_rhs(N[p],c), m_rhs(N[q],c) );
1515 [ - - ][ - - ]: 0 : alwn = min( m_rhs(N[p],c), m_rhs(N[q],c) );
1516 : : } else {
1517 [ + - ][ + - ]: 411120 : alwp = max( max(m_rhs(N[p],c), m_u(N[p],c)),
1518 [ + - ][ + - ]: 411120 : max(m_rhs(N[q],c), m_u(N[q],c)) );
1519 [ + - ][ + - ]: 411120 : alwn = min( min(m_rhs(N[p],c), m_u(N[p],c)),
1520 [ + - ][ + - ]: 411120 : min(m_rhs(N[q],c), m_u(N[q],c)) );
1521 : : }
1522 [ + - ][ + - ]: 411120 : m_q(N[p],a) = max(m_q(N[p],a), alwp);
1523 [ + - ][ + - ]: 411120 : m_q(N[p],b) = min(m_q(N[p],b), alwn);
1524 [ + - ][ + - ]: 411120 : m_q(N[q],a) = max(m_q(N[q],a), alwp);
1525 [ + - ][ + - ]: 411120 : m_q(N[q],b) = min(m_q(N[q],b), alwn);
1526 : : }
1527 : : }
1528 : : }
1529 : :
1530 : : // triangle superedges
1531 [ + + ]: 25460 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1532 : 22540 : const auto N = m_dsupedge[1].data() + e*3;
1533 [ + + ]: 90160 : for (std::size_t c=0; c<ncomp; ++c) {
1534 : 67620 : auto a = c*2;
1535 : 67620 : auto b = a+1;
1536 [ + + ]: 270480 : for (const auto& [p,q] : tk::lpoet) {
1537 : : tk::real alwp, alwn;
1538 [ - + ]: 202860 : if (g_cfg.get< tag::fctclip >()) {
1539 [ - - ][ - - ]: 0 : alwp = max( m_rhs(N[p],c), m_rhs(N[q],c) );
1540 [ - - ][ - - ]: 0 : alwn = min( m_rhs(N[p],c), m_rhs(N[q],c) );
1541 : : } else {
1542 [ + - ][ + - ]: 202860 : alwp = max( max(m_rhs(N[p],c), m_u(N[p],c)),
1543 [ + - ][ + - ]: 202860 : max(m_rhs(N[q],c), m_u(N[q],c)) );
1544 [ + - ][ + - ]: 202860 : alwn = min( min(m_rhs(N[p],c), m_u(N[p],c)),
1545 [ + - ][ + - ]: 202860 : min(m_rhs(N[q],c), m_u(N[q],c)) );
1546 : : }
1547 [ + - ][ + - ]: 202860 : m_q(N[p],a) = max(m_q(N[p],a), alwp);
1548 [ + - ][ + - ]: 202860 : m_q(N[p],b) = min(m_q(N[p],b), alwn);
1549 [ + - ][ + - ]: 202860 : m_q(N[q],a) = max(m_q(N[q],a), alwp);
1550 [ + - ][ + - ]: 202860 : m_q(N[q],b) = min(m_q(N[q],b), alwn);
1551 : : }
1552 : : }
1553 : : }
1554 : :
1555 : : // edges
1556 [ + + ]: 74000 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1557 : 71080 : const auto N = m_dsupedge[2].data() + e*2;
1558 [ + + ]: 284320 : for (std::size_t c=0; c<ncomp; ++c) {
1559 : 213240 : auto a = c*2;
1560 : 213240 : auto b = a+1;
1561 : : tk::real alwp, alwn;
1562 [ - + ]: 213240 : if (g_cfg.get< tag::fctclip >()) {
1563 [ - - ][ - - ]: 0 : alwp = max( m_rhs(N[0],c), m_rhs(N[1],c) );
1564 [ - - ][ - - ]: 0 : alwn = min( m_rhs(N[0],c), m_rhs(N[1],c) );
1565 : : } else {
1566 [ + - ][ + - ]: 213240 : alwp = max( max(m_rhs(N[0],c), m_u(N[0],c)),
1567 [ + - ][ + - ]: 213240 : max(m_rhs(N[1],c), m_u(N[1],c)) );
1568 [ + - ][ + - ]: 213240 : alwn = min( min(m_rhs(N[0],c), m_u(N[0],c)),
1569 [ + - ][ + - ]: 213240 : min(m_rhs(N[1],c), m_u(N[1],c)) );
1570 : : }
1571 [ + - ][ + - ]: 213240 : m_q(N[0],a) = max(m_q(N[0],a), alwp);
1572 [ + - ][ + - ]: 213240 : m_q(N[0],b) = min(m_q(N[0],b), alwn);
1573 [ + - ][ + - ]: 213240 : m_q(N[1],a) = max(m_q(N[1],a), alwp);
1574 [ + - ][ + - ]: 213240 : m_q(N[1],b) = min(m_q(N[1],b), alwn);
1575 : : }
1576 : : }
1577 : :
1578 : : // Communicate allowed limits contributions
1579 [ + + ]: 2920 : if (d->NodeCommMap().empty()) {
1580 : 20 : comalw_complete();
1581 : : } else {
1582 [ + + ]: 37780 : for (const auto& [c,n] : d->NodeCommMap()) {
1583 : 34880 : decltype(m_qc) exp;
1584 [ + - ][ + - ]: 165880 : for (auto g : n) exp[g] = m_q[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1585 [ + - ][ + - ]: 34880 : thisProxy[c].comalw( exp );
1586 : 34880 : }
1587 : : }
1588 : 2920 : ownalw_complete();
1589 : 2920 : }
1590 : :
1591 : : void
1592 : 34880 : ChoCG::comalw( const std::unordered_map< std::size_t,
1593 : : std::vector< tk::real > >& inalw )
1594 : : // *****************************************************************************
1595 : : // Receive allowed limits contributions on chare-boundaries
1596 : : //! \param[in] inalw Partial contributions of allowed limits contributions on
1597 : : //! chare-boundary nodes. Key: global mesh node IDs, value: allowed limit
1598 : : //! contributions.
1599 : : // *****************************************************************************
1600 : : {
1601 [ + + ]: 165880 : for (const auto& [g,alw] : inalw) {
1602 [ + - ]: 131000 : auto& q = m_qc[g];
1603 [ + - ]: 131000 : q.resize( alw.size() );
1604 [ + + ]: 524000 : for (std::size_t c=0; c<alw.size()/2; ++c) {
1605 : 393000 : auto a = c*2;
1606 : 393000 : auto b = a+1;
1607 : 393000 : q[a] = std::max( q[a], alw[a] );
1608 : 393000 : q[b] = std::min( q[b], alw[b] );
1609 : : }
1610 : : }
1611 : :
1612 : : // When we have heard from all chares we communicate with, this chare is done
1613 [ + + ]: 34880 : if (++m_nalw == Disc()->NodeCommMap().size()) {
1614 : 2900 : m_nalw = 0;
1615 : 2900 : comalw_complete();
1616 : : }
1617 : 34880 : }
1618 : :
1619 : : void
1620 : 2920 : ChoCG::lim()
1621 : : // *****************************************************************************
1622 : : // Compute limit coefficients
1623 : : // *****************************************************************************
1624 : : {
1625 [ + - ]: 2920 : auto d = Disc();
1626 : 2920 : const auto npoin = m_u.nunk();
1627 : 2920 : const auto ncomp = m_u.nprop();
1628 : 2920 : const auto& lid = d->Lid();
1629 : :
1630 : : using std::max;
1631 : : using std::min;
1632 : :
1633 : : // Combine own and communicated contributions to allowed limits
1634 [ + + ]: 56000 : for (const auto& [g,alw] : m_qc) {
1635 [ + - ]: 53080 : auto i = tk::cref_find( lid, g );
1636 [ + + ]: 212320 : for (std::size_t c=0; c<alw.size()/2; ++c) {
1637 : 159240 : auto a = c*2;
1638 : 159240 : auto b = a+1;
1639 [ + - ][ + - ]: 159240 : m_q(i,a) = max( m_q(i,a), alw[a] );
1640 [ + - ][ + - ]: 159240 : m_q(i,b) = min( m_q(i,b), alw[b] );
1641 : : }
1642 : : }
1643 : 2920 : tk::destroy(m_qc);
1644 : :
1645 : : // Finish computing allowed limits
1646 [ + + ]: 79900 : for (std::size_t i=0; i<npoin; ++i) {
1647 [ + + ]: 307920 : for (std::size_t c=0; c<ncomp; ++c) {
1648 : 230940 : auto a = c*2;
1649 : 230940 : auto b = a+1;
1650 [ + - ][ + - ]: 230940 : m_q(i,a) -= m_rhs(i,c);
1651 [ + - ][ + - ]: 230940 : m_q(i,b) -= m_rhs(i,c);
1652 : : }
1653 : : }
1654 : :
1655 : : // Limit coefficients, C
1656 : :
1657 [ + + ]: 79900 : for (std::size_t i=0; i<npoin; ++i) {
1658 [ + + ]: 307920 : for (std::size_t c=0; c<ncomp; ++c) {
1659 : 230940 : auto a = c*2;
1660 : 230940 : auto b = a+1;
1661 : 230940 : auto eps = std::numeric_limits< tk::real >::epsilon();
1662 [ + - ][ + - ]: 230940 : m_q(i,a) = m_p(i,a) < eps ? 0.0 : min(1.0, m_q(i,a)/m_p(i,a));
[ - - ][ - - ]
[ + - ]
1663 [ + - ][ + - ]: 230940 : m_q(i,b) = m_p(i,b) > -eps ? 0.0 : min(1.0, m_q(i,b)/m_p(i,b));
[ - - ][ - - ]
[ + - ]
1664 : : }
1665 : : }
1666 : :
1667 : : // Limited antidiffusive contributions
1668 : :
1669 : 2920 : auto ctau = g_cfg.get< tag::fctdif >();
1670 [ + - ]: 2920 : m_a.fill( 0.0 );
1671 : :
1672 [ + - ]: 2920 : auto fctsys = g_cfg.get< tag::fctsys >();
1673 [ - + ]: 2920 : for (auto& c : fctsys) --c;
1674 : :
1675 : : #if defined(__clang__)
1676 : : #pragma clang diagnostic push
1677 : : #pragma clang diagnostic ignored "-Wvla"
1678 : : #pragma clang diagnostic ignored "-Wvla-extension"
1679 : : #elif defined(STRICT_GNUC)
1680 : : #pragma GCC diagnostic push
1681 : : #pragma GCC diagnostic ignored "-Wvla"
1682 : : #endif
1683 : :
1684 : : // tetrahedron superedges
1685 [ + + ]: 25760 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1686 : 22840 : const auto N = m_dsupedge[0].data() + e*4;
1687 : 22840 : const auto D = m_dsupint[0].data();
1688 : 22840 : auto C = m_dsuplim[0].data();
1689 : 22840 : std::size_t i = 0;
1690 [ + + ]: 296920 : for (const auto& [p,q] : tk::lpoed) {
1691 : 137040 : auto dif = D[(e*6+i)*5+3];
1692 : 137040 : auto coef = C + (e*6+i)*ncomp;
1693 : 137040 : tk::real aec[ncomp];
1694 [ + + ]: 548160 : for (std::size_t c=0; c<ncomp; ++c) {
1695 [ + - ][ + - ]: 411120 : aec[c] = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1696 : 411120 : auto a = c*2;
1697 : 411120 : auto b = a+1;
1698 [ - + ][ - - ]: 411120 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
[ + - ]
1699 [ - + ][ - - ]: 411120 : aec[c] > 0.0 ? m_q(N[q],a) : m_q(N[q],b) );
[ + - ]
1700 : : }
1701 : 137040 : tk::real cs = 1.0;
1702 [ - + ]: 137040 : for (auto c : fctsys) cs = min( cs, coef[c] );
1703 [ - + ]: 137040 : for (auto c : fctsys) coef[c] = cs;
1704 [ + + ]: 548160 : for (std::size_t c=0; c<ncomp; ++c) {
1705 : 411120 : aec[c] *= coef[c];
1706 [ + - ]: 411120 : m_a(N[p],c) -= aec[c];
1707 [ + - ]: 411120 : m_a(N[q],c) += aec[c];
1708 : : }
1709 : 137040 : ++i;
1710 : 137040 : }
1711 : : }
1712 : :
1713 : : // triangle superedges
1714 [ + + ]: 25460 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1715 : 22540 : const auto N = m_dsupedge[1].data() + e*3;
1716 : 22540 : const auto D = m_dsupint[1].data();
1717 : 22540 : auto C = m_dsuplim[0].data();
1718 : 22540 : std::size_t i = 0;
1719 [ + + ]: 157780 : for (const auto& [p,q] : tk::lpoet) {
1720 : 67620 : auto dif = D[(e*3+i)*5+3];
1721 : 67620 : auto coef = C + (e*3+i)*ncomp;
1722 : 67620 : tk::real aec[ncomp];
1723 [ + + ]: 270480 : for (std::size_t c=0; c<ncomp; ++c) {
1724 [ + - ][ + - ]: 202860 : aec[c] = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1725 : 202860 : auto a = c*2;
1726 : 202860 : auto b = a+1;
1727 [ - + ][ - - ]: 202860 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
[ + - ]
1728 [ - + ][ - - ]: 202860 : aec[c] > 0.0 ? m_q(N[q],a) : m_q(N[q],b) );
[ + - ]
1729 : : }
1730 : 67620 : tk::real cs = 1.0;
1731 [ - + ]: 67620 : for (auto c : fctsys) cs = min( cs, coef[c] );
1732 [ - + ]: 67620 : for (auto c : fctsys) coef[c] = cs;
1733 [ + + ]: 270480 : for (std::size_t c=0; c<ncomp; ++c) {
1734 : 202860 : aec[c] *= coef[c];
1735 [ + - ]: 202860 : m_a(N[p],c) -= aec[c];
1736 [ + - ]: 202860 : m_a(N[q],c) += aec[c];
1737 : : }
1738 : 67620 : ++i;
1739 : 67620 : }
1740 : : }
1741 : :
1742 : : // edges
1743 [ + + ]: 74000 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1744 : 71080 : const auto N = m_dsupedge[2].data() + e*2;
1745 : 71080 : const auto dif = m_dsupint[2][e*5+3];
1746 : 71080 : auto coef = m_dsuplim[2].data() + e*ncomp;
1747 : 71080 : tk::real aec[ncomp];
1748 [ + + ]: 284320 : for (std::size_t c=0; c<ncomp; ++c) {
1749 [ + - ][ + - ]: 213240 : aec[c] = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1750 : 213240 : auto a = c*2;
1751 : 213240 : auto b = a+1;
1752 [ - + ][ - - ]: 213240 : coef[c] = min( aec[c] < 0.0 ? m_q(N[0],a) : m_q(N[0],b),
[ + - ]
1753 [ - + ][ - - ]: 213240 : aec[c] > 0.0 ? m_q(N[1],a) : m_q(N[1],b) );
[ + - ]
1754 : : }
1755 : 71080 : tk::real cs = 1.0;
1756 [ - + ]: 71080 : for (auto c : fctsys) cs = min( cs, coef[c] );
1757 [ - + ]: 71080 : for (auto c : fctsys) coef[c] = cs;
1758 [ + + ]: 284320 : for (std::size_t c=0; c<ncomp; ++c) {
1759 : 213240 : aec[c] *= coef[c];
1760 [ + - ]: 213240 : m_a(N[0],c) -= aec[c];
1761 [ + - ]: 213240 : m_a(N[1],c) += aec[c];
1762 : : }
1763 : 71080 : }
1764 : :
1765 : : #if defined(__clang__)
1766 : : #pragma clang diagnostic pop
1767 : : #elif defined(STRICT_GNUC)
1768 : : #pragma GCC diagnostic pop
1769 : : #endif
1770 : :
1771 : : // Communicate limited antidiffusive contributions
1772 [ + + ]: 2920 : if (d->NodeCommMap().empty()){
1773 [ + - ]: 20 : comlim_complete();
1774 : : } else {
1775 [ + + ]: 37780 : for (const auto& [c,n] : d->NodeCommMap()) {
1776 : 34880 : decltype(m_ac) exp;
1777 [ + - ][ + - ]: 165880 : for (auto g : n) exp[g] = m_a[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1778 [ + - ][ + - ]: 34880 : thisProxy[c].comlim( exp );
1779 : 34880 : }
1780 : : }
1781 [ + - ]: 2920 : ownlim_complete();
1782 : 2920 : }
1783 : :
1784 : : void
1785 : 34880 : ChoCG::comlim( const std::unordered_map< std::size_t,
1786 : : std::vector< tk::real > >& inlim )
1787 : : // *****************************************************************************
1788 : : // Receive limited antidiffusive contributions on chare-boundaries
1789 : : //! \param[in] inlim Partial contributions of limited contributions on
1790 : : //! chare-boundary nodes. Key: global mesh node IDs, value: limited
1791 : : //! contributions.
1792 : : // *****************************************************************************
1793 : : {
1794 : : using tk::operator+=;
1795 [ + - ][ + - ]: 165880 : for (const auto& [g,a] : inlim) m_ac[g] += a;
[ + + ]
1796 : :
1797 : : // When we have heard from all chares we communicate with, this chare is done
1798 [ + + ]: 34880 : if (++m_nlim == Disc()->NodeCommMap().size()) {
1799 : 2900 : m_nlim = 0;
1800 : 2900 : comlim_complete();
1801 : : }
1802 : 34880 : }
1803 : :
1804 : : void
1805 : 8110 : ChoCG::BC( tk::Fields& u, tk::real t )
1806 : : // *****************************************************************************
1807 : : // Apply boundary conditions
1808 : : //! \param[in,out] u Solution to apply BCs to
1809 : : //! \param[in] t Physical time
1810 : : // *****************************************************************************
1811 : : {
1812 : 8110 : auto d = Disc();
1813 : :
1814 : 8110 : physics::dirbc( u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
1815 : 8110 : physics::symbc( u, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
1816 : 8110 : physics::noslipbc( u, m_noslipbcnodes, /*pos=*/0 );
1817 : 8110 : }
1818 : :
1819 : : void
1820 : 3630 : ChoCG::dt()
1821 : : // *****************************************************************************
1822 : : // Compute time step size
1823 : : // *****************************************************************************
1824 : : {
1825 [ + - ]: 3630 : auto d = Disc();
1826 : 3630 : const auto& vol = d->Vol();
1827 : :
1828 : 3630 : tk::real mindt = std::numeric_limits< tk::real >::max();
1829 : 3630 : auto const_dt = g_cfg.get< tag::dt >();
1830 : 3630 : auto eps = std::numeric_limits< tk::real >::epsilon();
1831 : :
1832 : : // use constant dt if configured
1833 [ + + ]: 3630 : if (std::abs(const_dt) > eps) {
1834 : :
1835 : : // cppcheck-suppress redundantInitialization
1836 : 2920 : mindt = const_dt;
1837 : :
1838 : : } else {
1839 : :
1840 : 710 : auto cfl = g_cfg.get< tag::cfl >();
1841 : 710 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1842 : 710 : auto large = std::numeric_limits< tk::real >::max();
1843 : :
1844 [ + + ]: 278930 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1845 [ + - ]: 278220 : auto u = m_u(i,0);
1846 [ + - ]: 278220 : auto v = m_u(i,1);
1847 [ + - ]: 278220 : auto w = m_u(i,2);
1848 : 278220 : auto vel = std::sqrt( u*u + v*v + w*w );
1849 : 278220 : auto L = std::cbrt( vol[i] );
1850 : 278220 : auto euler_dt = L / std::max( vel, 1.0e-8 );
1851 : 278220 : mindt = std::min( mindt, euler_dt );
1852 [ + + ]: 278220 : auto visc_dt = mu > eps ? L * L / mu : large;
1853 : 278220 : mindt = std::min( mindt, visc_dt );
1854 : : }
1855 : 710 : mindt *= cfl;
1856 : :
1857 : : }
1858 : :
1859 : : // Actiavate SDAG waits for next time step stage
1860 [ + - ][ + - ]: 3630 : thisProxy[ thisIndex ].wait4rhs();
1861 [ + - ][ + - ]: 3630 : thisProxy[ thisIndex ].wait4aec();
1862 [ + - ][ + - ]: 3630 : thisProxy[ thisIndex ].wait4alw();
1863 [ + - ][ + - ]: 3630 : thisProxy[ thisIndex ].wait4sol();
1864 [ + - ][ + - ]: 3630 : thisProxy[ thisIndex ].wait4div();
1865 [ + - ][ + - ]: 3630 : thisProxy[ thisIndex ].wait4sgrad();
1866 [ + - ][ + - ]: 3630 : thisProxy[ thisIndex ].wait4step();
1867 : :
1868 : : // Contribute to minimum dt across all chares and advance to next step
1869 [ + - ]: 3630 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
1870 [ + - ][ + - ]: 7260 : CkCallback(CkReductionTarget(ChoCG,advance), thisProxy) );
1871 : 3630 : }
1872 : :
1873 : : void
1874 : 3630 : ChoCG::advance( tk::real newdt )
1875 : : // *****************************************************************************
1876 : : // Advance equations to next time step
1877 : : //! \param[in] newdt The smallest dt across the whole problem
1878 : : // *****************************************************************************
1879 : : {
1880 : : // Set new time step size
1881 : 3630 : Disc()->setdt( newdt );
1882 : :
1883 : : // Compute lhs and rhs of transport equations
1884 : 3630 : lhs();
1885 : 3630 : rhs();
1886 : 3630 : }
1887 : :
1888 : : void
1889 : 3630 : ChoCG::lhs()
1890 : : // *****************************************************************************
1891 : : // Fill lhs matrix of transport equations
1892 : : // *****************************************************************************
1893 : : {
1894 : 3630 : auto theta = g_cfg.get< tag::theta >();
1895 [ + + ]: 3630 : if (theta < std::numeric_limits< tk::real >::epsilon()) return;
1896 : :
1897 : 40 : auto d = Disc();
1898 : 40 : const auto& inpoel = d->Inpoel();
1899 : 40 : const auto& coord = d->Coord();
1900 : 40 : const auto& X = coord[0];
1901 : 40 : const auto& Y = coord[1];
1902 : 40 : const auto& Z = coord[2];
1903 : 40 : const auto ncomp = m_u.nprop();
1904 : 40 : const auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1905 : :
1906 : 40 : auto dt = d->Dt();
1907 : 40 : auto& A = Lhs();
1908 : 40 : A.zero();
1909 : :
1910 [ + + ]: 60040 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1911 : 60000 : const auto N = inpoel.data() + e*4;
1912 : : const std::array< tk::real, 3 >
1913 : 60000 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
1914 : 60000 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
1915 : 60000 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
1916 : 60000 : const auto J = tk::triple( ba, ca, da ); // J = 6V
1917 [ - + ][ - - ]: 60000 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
1918 : : std::array< std::array< tk::real, 3 >, 4 > grad;
1919 : 60000 : grad[1] = tk::cross( ca, da );
1920 : 60000 : grad[2] = tk::cross( da, ba );
1921 : 60000 : grad[3] = tk::cross( ba, ca );
1922 [ + + ]: 240000 : for (std::size_t i=0; i<3; ++i)
1923 : 180000 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
1924 [ + + ]: 300000 : for (std::size_t a=0; a<4; ++a) {
1925 [ + + ]: 1200000 : for (std::size_t b=0; b<4; ++b) {
1926 [ + + ]: 960000 : auto v = J/dt/120.0 * ((a == b) ? 2.0 : 1.0);
1927 : 960000 : v += theta * mu * tk::dot(grad[a],grad[b]) / J / 6.0;
1928 [ + - ][ + + ]: 3840000 : for (std::size_t c=0; c<ncomp; ++c) A(N[a],N[b],c) -= v;
1929 : : }
1930 : : //for (std::size_t c=0; c<ncomp; ++c) A(N[a],N[a],c) -= J/dt/24.0;
1931 : : }
1932 : : }
1933 : : }
1934 : :
1935 : : void
1936 : 3750 : ChoCG::rhs()
1937 : : // *****************************************************************************
1938 : : // Compute right-hand side of transport equations
1939 : : // *****************************************************************************
1940 : : {
1941 : 3750 : auto d = Disc();
1942 : 3750 : const auto& lid = d->Lid();
1943 : :
1944 : : // Compute own portion of right-hand side for all equations
1945 : 3750 : auto dt = m_rkcoef[m_stage] * d->Dt();
1946 : 3750 : chorin::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel,
1947 : 3750 : dt, m_pr, m_u, m_vgrad, m_pgrad, m_rhs );
1948 : :
1949 : : // Communicate rhs to other chares on chare-boundary
1950 [ + + ]: 3750 : if (d->NodeCommMap().empty()) {
1951 : 290 : comrhs_complete();
1952 : : } else {
1953 [ + + ]: 39380 : for (const auto& [c,n] : d->NodeCommMap()) {
1954 : 35920 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1955 [ + - ][ + - ]: 184500 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1956 [ + - ][ + - ]: 35920 : thisProxy[c].comrhs( exp );
1957 : 35920 : }
1958 : : }
1959 : 3750 : ownrhs_complete();
1960 : 3750 : }
1961 : :
1962 : : void
1963 : 35920 : ChoCG::comrhs(
1964 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1965 : : // *****************************************************************************
1966 : : // Receive contributions to right-hand side vector on chare-boundaries
1967 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1968 : : //! global mesh node IDs, value: contributions for all scalar components.
1969 : : //! \details This function receives contributions to m_rhs, which stores the
1970 : : //! right hand side vector at mesh nodes. While m_rhs stores own
1971 : : //! contributions, m_rhsc collects the neighbor chare contributions during
1972 : : //! communication. This way work on m_rhs and m_rhsc is overlapped.
1973 : : // *****************************************************************************
1974 : : {
1975 : : using tk::operator+=;
1976 [ + - ][ + - ]: 184500 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
[ + + ]
1977 : :
1978 : : // When we have heard from all chares we communicate with, this chare is done
1979 [ + + ]: 35920 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
1980 : 3460 : m_nrhs = 0;
1981 : 3460 : comrhs_complete();
1982 : : }
1983 : 35920 : }
1984 : :
1985 : : void
1986 : 3750 : ChoCG::fct()
1987 : : // *****************************************************************************
1988 : : // Continue with flux-corrected transport if enabled
1989 : : // *****************************************************************************
1990 : : {
1991 : : // Combine own and communicated contributions to rhs
1992 [ + - ][ + - ]: 3750 : const auto lid = Disc()->Lid();
1993 [ + + ]: 73020 : for (const auto& [g,r] : m_rhsc) {
1994 [ + - ]: 69270 : auto i = tk::cref_find( lid, g );
1995 [ + - ][ + + ]: 277080 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1996 : : }
1997 : 3750 : tk::destroy(m_rhsc);
1998 : :
1999 : 3750 : auto eps = std::numeric_limits< tk::real >::epsilon();
2000 [ + + ][ + + ]: 3750 : if (g_cfg.get< tag::theta >() < eps and g_cfg.get< tag::fct >())
[ + + ]
2001 [ + - ]: 2920 : aec();
2002 : : else
2003 [ + - ]: 830 : solve();
2004 : 3750 : }
2005 : :
2006 : : void
2007 : : // cppcheck-suppress unusedFunction
2008 : 3750 : ChoCG::solve()
2009 : : // *****************************************************************************
2010 : : // Advance systems of equations
2011 : : // *****************************************************************************
2012 : : {
2013 : 3750 : auto d = Disc();
2014 : 3750 : const auto npoin = m_u.nunk();
2015 : 3750 : const auto ncomp = m_u.nprop();
2016 : 3750 : const auto& vol = d->Vol();
2017 : :
2018 : : // Combine own and communicated contributions to limited antidiffusive
2019 : : // contributions
2020 [ + + ]: 56830 : for (const auto& [g,a] : m_ac) {
2021 [ + - ]: 53080 : auto i = tk::cref_find( d->Lid(), g );
2022 [ + - ][ + + ]: 212320 : for (std::size_t c=0; c<a.size(); ++c) m_a(i,c) += a[c];
2023 : : }
2024 : 3750 : tk::destroy(m_ac);
2025 : :
2026 [ + + ]: 3750 : if (m_stage == 0) m_un = m_u;
2027 : :
2028 [ + + ]: 3750 : if (g_cfg.get< tag::fct >()) {
2029 : :
2030 : : // Apply limited antidiffusive contributions to low-order solution
2031 [ + + ]: 79900 : for (std::size_t i=0; i<npoin; ++i) {
2032 [ + + ]: 307920 : for (std::size_t c=0; c<ncomp; ++c) {
2033 : 230940 : m_a(i,c) = m_rhs(i,c) + m_a(i,c)/vol[i];
2034 : : }
2035 : : }
2036 : : // Continue to advective-diffusive prediction
2037 : 2920 : pred();
2038 : :
2039 : : } else {
2040 : :
2041 : 830 : auto eps = std::numeric_limits<tk::real>::epsilon();
2042 [ + + ][ - + ]: 830 : if (g_cfg.get< tag::theta >() < eps || m_stage+1 < m_rkcoef.size()) {
[ + + ]
2043 : :
2044 : : // Apply rhs in explicit solve
2045 : 790 : auto dt = m_rkcoef[m_stage] * d->Dt();
2046 [ + + ]: 399970 : for (std::size_t i=0; i<npoin; ++i) {
2047 [ + + ]: 1596720 : for (std::size_t c=0; c<ncomp; ++c) {
2048 : 1197540 : m_a(i,c) = m_un(i,c) - dt*m_rhs(i,c)/vol[i];
2049 : : }
2050 : : }
2051 : : // Continue to advective-diffusive prediction
2052 : 790 : pred();
2053 : :
2054 : : } else {
2055 : :
2056 : : // Configure Dirichlet BCs
2057 : : std::unordered_map< std::size_t,
2058 : 40 : std::vector< std::pair< int, tk::real > > > dirbc;
2059 : 40 : std::size_t nmask = ncomp + 1;
2060 [ - + ][ - - ]: 40 : Assert( m_dirbcmask.size() % nmask == 0, "Size mismatch" );
[ - - ][ - - ]
2061 [ + + ]: 24520 : for (std::size_t i=0; i<m_dirbcmask.size()/nmask; ++i) {
2062 : 24480 : auto p = m_dirbcmask[i*nmask+0]; // local node id
2063 [ + - ]: 24480 : auto& bc = dirbc[p];
2064 [ + - ]: 24480 : bc.resize( ncomp );
2065 [ + + ]: 97920 : for (std::size_t c=0; c<ncomp; ++c) {
2066 : 73440 : bc[c] = { m_dirbcmask[i*nmask+1+c], 0.0 };
2067 : : }
2068 : : }
2069 [ + + ]: 8200 : for (auto p : m_noslipbcnodes) {
2070 [ + - ]: 8160 : auto& bc = dirbc[p];
2071 [ + - ]: 8160 : bc.resize( ncomp );
2072 [ + + ]: 32640 : for (std::size_t c=0; c<ncomp; ++c) {
2073 : 24480 : bc[c] = { 1, 0.0 };
2074 : : }
2075 : : }
2076 : :
2077 : : // Initialize semi-implicit momentum/transport solve
2078 : 40 : const auto& pc = g_cfg.get< tag::mom_pc >();
2079 [ + - ][ + - ]: 80 : m_cgmom[ thisIndex ].ckLocal()->init( {}, m_rhs.vec(), {}, dirbc, pc,
[ + - ]
2080 [ + - ][ + - ]: 80 : CkCallback( CkIndex_ChoCG::msolve(), thisProxy[thisIndex] ) );
[ + - ]
2081 : :
2082 : 40 : }
2083 : :
2084 : : }
2085 : 3750 : }
2086 : :
2087 : : void
2088 : 40 : ChoCG::msolve()
2089 : : // *****************************************************************************
2090 : : // Solve for momentum/transport system of equations
2091 : : // *****************************************************************************
2092 : : {
2093 : 40 : auto iter = g_cfg.get< tag::mom_iter >();
2094 : 40 : auto tol = g_cfg.get< tag::mom_tol >();
2095 : 40 : auto verbose = g_cfg.get< tag::mom_verbose >();
2096 : :
2097 [ + - ][ + - ]: 80 : m_cgmom[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose,
[ + - ]
2098 [ + - ][ + - ]: 80 : CkCallback( CkIndex_ChoCG::msolved(), thisProxy[thisIndex] ) );
[ + - ]
2099 : 40 : }
2100 : :
2101 : : void
2102 : 40 : ChoCG::msolved()
2103 : : // *****************************************************************************
2104 : : // Continue after momentum/transport solve in semi-implcit solve
2105 : : // *****************************************************************************
2106 : : {
2107 : 40 : auto d = Disc();
2108 : 40 : const auto npoin = m_u.nunk();
2109 : 40 : const auto ncomp = m_u.nprop();
2110 : :
2111 [ + + ][ + - ]: 40 : if (thisIndex == 0) d->mit( m_cgmom[ thisIndex ].ckLocal()->it() );
[ + - ][ + - ]
2112 : :
2113 : : // Update momentum/transport solution in semi-implicit solve
2114 [ + - ][ + - ]: 40 : auto& du = m_cgmom[ thisIndex ].ckLocal()->solution();
2115 [ + + ]: 24520 : for (std::size_t i=0; i<npoin; ++i) {
2116 [ + + ]: 97920 : for (std::size_t c=0; c<ncomp; ++c) {
2117 : 73440 : m_a(i,c) = m_un(i,c) + du[i*ncomp+c];
2118 : : }
2119 : : }
2120 : :
2121 : : // Continue to advective-diffusive prediction
2122 : 40 : pred();
2123 : 40 : }
2124 : :
2125 : : void
2126 : 3750 : ChoCG::pred()
2127 : : // *****************************************************************************
2128 : : // Compute advective-diffusive prediction of momentum/transport
2129 : : // *****************************************************************************
2130 : : {
2131 [ + - ]: 3750 : auto d = Disc();
2132 : :
2133 : : // Configure and apply scalar source to solution (if defined)
2134 [ + - ]: 3750 : auto src = problems::PHYS_SRC();
2135 [ - + ][ - - ]: 3750 : if (src) src( d->Coord(), d->T(), m_a );
2136 : :
2137 : : // Enforce boundary conditions
2138 [ + - ]: 3750 : BC( m_a, d->T() + m_rkcoef[m_stage] * d->Dt() );
2139 : :
2140 : : // Update momentum/transport solution
2141 [ + - ]: 3750 : m_u = m_a;
2142 [ + - ]: 3750 : m_a.fill( 0.0 );
2143 : :
2144 : : // Compute velocity gradients if needed
2145 [ + + ]: 3750 : if (g_cfg.get< tag::flux >() == "damp4") {
2146 [ + - ][ + - ]: 3080 : thisProxy[ thisIndex ].wait4vgrad();
2147 [ + - ]: 3080 : velgrad();
2148 : : } else {
2149 [ + - ]: 670 : corr();
2150 : : }
2151 : 3750 : }
2152 : :
2153 : : void
2154 : 3750 : ChoCG::corr()
2155 : : // *****************************************************************************
2156 : : // Compute pressure correction
2157 : : // *****************************************************************************
2158 : : {
2159 : : // Finalize computing velocity gradients
2160 [ + + ]: 3750 : if (g_cfg.get< tag::flux >() == "damp4") fingrad( m_vgrad, m_vgradc );
2161 : :
2162 [ + + ]: 3750 : if (++m_stage < m_rkcoef.size()) {
2163 : :
2164 : : // Activate SDAG wait for next time step stage
2165 [ + - ][ + - ]: 120 : thisProxy[ thisIndex ].wait4rhs();
2166 : : // Continue to next time stage of momentum/transport prediction
2167 : 120 : rhs();
2168 : :
2169 : : } else {
2170 : :
2171 : : // Reset Runge-Kutta stage counter
2172 : 3630 : m_stage = 0;
2173 : : // Continue to pressure correction and projection
2174 : 3630 : div( m_u );
2175 : :
2176 : : }
2177 : 3750 : }
2178 : :
2179 : : void
2180 : 3662 : ChoCG::diag()
2181 : : // *****************************************************************************
2182 : : // Compute diagnostics
2183 : : // *****************************************************************************
2184 : : {
2185 : 3662 : auto d = Disc();
2186 : :
2187 : : // Increase number of iterations and physical time
2188 : 3662 : d->next();
2189 : :
2190 : : // Compute diagnostics, e.g., residuals
2191 : 3662 : auto diag_iter = g_cfg.get< tag::diag_iter >();
2192 [ + - ][ + - ]: 3662 : const auto& dp = m_cgpre[ thisIndex ].ckLocal()->solution();
2193 : 3662 : auto diagnostics = m_diag.precompute( *d, m_u, m_un, m_pr, dp, diag_iter );
2194 : :
2195 : : // Evaluate residuals
2196 [ - + ][ - - ]: 3662 : if (!diagnostics) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ - - ]
2197 : 3662 : }
2198 : :
2199 : : void
2200 : 3662 : ChoCG::evalres( const std::vector< tk::real >& )
2201 : : // *****************************************************************************
2202 : : // Evaluate residuals
2203 : : // *****************************************************************************
2204 : : {
2205 : 3662 : refine();
2206 : 3662 : }
2207 : :
2208 : : void
2209 : 3662 : ChoCG::refine()
2210 : : // *****************************************************************************
2211 : : // Optionally refine/derefine mesh
2212 : : // *****************************************************************************
2213 : : {
2214 : 3662 : auto d = Disc();
2215 : :
2216 : : // See if this is the last time step
2217 [ + + ]: 3662 : if (d->finished()) m_finished = 1;
2218 : :
2219 : 3662 : auto dtref = g_cfg.get< tag::href_dt >();
2220 : 3662 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
2221 : :
2222 : : // if t>0 refinement enabled and we hit the frequency
2223 [ - + ][ - - ]: 3662 : if (dtref && !(d->It() % dtfreq)) { // refine
[ - + ]
2224 : :
2225 : 0 : d->refined() = 1;
2226 : 0 : d->startvol();
2227 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
2228 : :
2229 : : // Activate SDAG waits for re-computing the integrals
2230 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
2231 : :
2232 : : } else { // do not refine
2233 : :
2234 : 3662 : d->refined() = 0;
2235 : 3662 : feop_complete();
2236 : 3662 : resize_complete();
2237 : :
2238 : : }
2239 : 3662 : }
2240 : :
2241 : : void
2242 : 0 : ChoCG::resizePostAMR(
2243 : : const std::vector< std::size_t >& /*ginpoel*/,
2244 : : const tk::UnsMesh::Chunk& chunk,
2245 : : const tk::UnsMesh::Coords& coord,
2246 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
2247 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
2248 : : const std::set< std::size_t >& removedNodes,
2249 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
2250 : : nodeCommMap,
2251 : : const std::map< int, std::vector< std::size_t > >& bface,
2252 : : const std::map< int, std::vector< std::size_t > >& bnode,
2253 : : const std::vector< std::size_t >& triinpoel )
2254 : : // *****************************************************************************
2255 : : // Receive new mesh from Refiner
2256 : : //! \param[in] ginpoel Mesh connectivity with global node ids
2257 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
2258 : : //! \param[in] coord New mesh node coordinates
2259 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
2260 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
2261 : : //! \param[in] removedNodes Newly removed mesh node local ids
2262 : : //! \param[in] nodeCommMap New node communication map
2263 : : //! \param[in] bface Boundary-faces mapped to side set ids
2264 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
2265 : : //! \param[in] triinpoel Boundary-face connectivity
2266 : : // *****************************************************************************
2267 : : {
2268 [ - - ]: 0 : auto d = Disc();
2269 : :
2270 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
2271 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
2272 : :
2273 : : // Resize mesh data structures after mesh refinement
2274 [ - - ]: 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
2275 : :
2276 [ - - ][ - - ]: 0 : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
2277 : : "Incorrect vector length post-AMR: expected length after resizing = " +
2278 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
2279 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
2280 : :
2281 : : // Remove newly removed nodes from solution vectors
2282 [ - - ]: 0 : m_u.rm( removedNodes );
2283 : : //m_pr.rm( removedNodes );
2284 [ - - ]: 0 : m_rhs.rm( removedNodes );
2285 : :
2286 : : // Resize auxiliary solution vectors
2287 : 0 : auto npoin = coord[0].size();
2288 [ - - ]: 0 : m_u.resize( npoin );
2289 [ - - ]: 0 : m_pr.resize( npoin );
2290 [ - - ]: 0 : m_rhs.resize( npoin );
2291 : :
2292 : : // Update solution on new mesh
2293 [ - - ]: 0 : for (const auto& n : addedNodes)
2294 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
2295 [ - - ][ - - ]: 0 : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
[ - - ][ - - ]
2296 [ - - ][ - - ]: 0 : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
[ - - ][ - - ]
[ - - ]
2297 : : "Indices of parent-edge nodes out of bounds post-AMR");
2298 [ - - ][ - - ]: 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
[ - - ]
2299 : : }
2300 : :
2301 : : // Update physical-boundary node-, face-, and element lists
2302 [ - - ]: 0 : m_bnode = bnode;
2303 [ - - ]: 0 : m_bface = bface;
2304 [ - - ]: 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
2305 : :
2306 : 0 : auto meshid = d->MeshId();
2307 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2308 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
2309 : 0 : }
2310 : :
2311 : : void
2312 : 874 : ChoCG::writeFields( CkCallback cb )
2313 : : // *****************************************************************************
2314 : : // Output mesh-based fields to file
2315 : : //! \param[in] cb Function to continue with after the write
2316 : : // *****************************************************************************
2317 : : {
2318 [ + + ][ + - ]: 874 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
2319 : :
2320 [ + - ]: 290 : auto d = Disc();
2321 : :
2322 : : // Field output
2323 : :
2324 : : std::vector< std::string > nodefieldnames{
2325 [ + - ][ + + ]: 1740 : "velocityx", "velocityy", "velocityz", "divergence", "pressure" };
[ - - ]
2326 : :
2327 : 290 : std::vector< std::vector< tk::real > > nodefields;
2328 : :
2329 [ + - ][ + - ]: 290 : nodefields.push_back( m_u.extract(0) );
2330 [ + - ][ + - ]: 290 : nodefields.push_back( m_u.extract(1) );
2331 [ + - ][ + - ]: 290 : nodefields.push_back( m_u.extract(2) );
2332 : :
2333 : : // Divide weak result by nodal volume
2334 : 290 : const auto& vol = d->Vol();
2335 [ + + ]: 80853 : for (std::size_t i=0; i<m_div.size(); ++i) m_div[i] /= vol[i];
2336 [ + - ]: 290 : nodefields.push_back( m_div );
2337 : :
2338 [ + - ]: 290 : nodefields.push_back( m_pr ) ;
2339 : :
2340 : : //nodefieldnames.push_back( "dp/dx" );
2341 : : //nodefieldnames.push_back( "dp/dy" );
2342 : : //nodefieldnames.push_back( "dp/dz" );
2343 : : //nodefields.push_back( m_pgrad.extract(0) );
2344 : : //nodefields.push_back( m_pgrad.extract(1) );
2345 : : //nodefields.push_back( m_pgrad.extract(2) );
2346 : :
2347 : : //nodefieldnames.push_back( "fx" );
2348 : : //nodefieldnames.push_back( "fy" );
2349 : : //nodefieldnames.push_back( "fz" );
2350 : : //nodefields.push_back( m_flux.extract(0) );
2351 : : //nodefields.push_back( m_flux.extract(1) );
2352 : : //nodefields.push_back( m_flux.extract(2) );
2353 : :
2354 : : //nodefieldnames.push_back( "du/dx" );
2355 : : //nodefieldnames.push_back( "du/dy" );
2356 : : //nodefieldnames.push_back( "du/dz" );
2357 : : //nodefieldnames.push_back( "dv/dx" );
2358 : : //nodefieldnames.push_back( "dv/dy" );
2359 : : //nodefieldnames.push_back( "dv/dz" );
2360 : : //nodefieldnames.push_back( "dw/dx" );
2361 : : //nodefieldnames.push_back( "dw/dy" );
2362 : : //nodefieldnames.push_back( "dw/dz" );
2363 : : //nodefields.push_back( m_vgrad.extract(0) );
2364 : : //nodefields.push_back( m_vgrad.extract(1) );
2365 : : //nodefields.push_back( m_vgrad.extract(2) );
2366 : : //nodefields.push_back( m_vgrad.extract(3) );
2367 : : //nodefields.push_back( m_vgrad.extract(4) );
2368 : : //nodefields.push_back( m_vgrad.extract(5) );
2369 : : //nodefields.push_back( m_vgrad.extract(6) );
2370 : : //nodefields.push_back( m_vgrad.extract(7) );
2371 : : //nodefields.push_back( m_vgrad.extract(8) );
2372 : :
2373 : : // also output analytic pressure (if defined)
2374 [ + - ]: 290 : auto pressure_sol = problems::PRESSURE_SOL();
2375 [ + + ]: 290 : if (pressure_sol) {
2376 : 64 : const auto& coord = d->Coord();
2377 : 64 : const auto& x = coord[0];
2378 : 64 : const auto& y = coord[1];
2379 : 64 : const auto& z = coord[2];
2380 [ + - ]: 64 : auto ap = m_pr;
2381 [ + + ]: 4348 : for (std::size_t i=0; i<ap.size(); ++i) {
2382 [ + - ]: 4284 : ap[i] = pressure_sol( x[i], y[i], z[i] );
2383 : : }
2384 [ + - ][ + - ]: 64 : nodefieldnames.push_back( "analytic" );
2385 [ + - ]: 64 : nodefields.push_back( ap );
2386 : 64 : }
2387 : :
2388 [ - + ][ - - ]: 290 : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
2389 : :
2390 : : // Surface output
2391 : :
2392 : 290 : std::vector< std::string > nodesurfnames;
2393 : 290 : std::vector< std::vector< tk::real > > nodesurfs;
2394 : :
2395 : 290 : const auto& f = g_cfg.get< tag::fieldout >();
2396 : :
2397 [ - + ]: 290 : if (!f.empty()) {
2398 : 0 : std::size_t ncomp = 5;
2399 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "velocityx" );
2400 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "velocityy" );
2401 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "velocityz" );
2402 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "divergence" );
2403 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "pressure" );
2404 : :
2405 [ - - ]: 0 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
2406 [ - - ]: 0 : std::set< int > outsets( begin(f), end(f) );
2407 [ - - ]: 0 : for (auto sideset : outsets) {
2408 [ - - ]: 0 : auto b = bnode.find(sideset);
2409 [ - - ]: 0 : if (b == end(bnode)) continue;
2410 : 0 : const auto& nodes = b->second;
2411 : 0 : auto i = nodesurfs.size();
2412 [ - - ]: 0 : nodesurfs.insert( end(nodesurfs), ncomp,
2413 [ - - ]: 0 : std::vector< tk::real >( nodes.size() ) );
2414 : 0 : std::size_t j = 0;
2415 [ - - ]: 0 : for (auto n : nodes) {
2416 [ - - ]: 0 : const auto s = m_u[n];
2417 : 0 : std::size_t p = 0;
2418 : 0 : nodesurfs[i+(p++)][j] = s[0];
2419 : 0 : nodesurfs[i+(p++)][j] = s[1];
2420 : 0 : nodesurfs[i+(p++)][j] = s[2];
2421 : 0 : nodesurfs[i+(p++)][j] = m_div[n];
2422 : 0 : nodesurfs[i+(p++)][j] = m_pr[n];
2423 : 0 : ++j;
2424 : 0 : }
2425 : : }
2426 : 0 : }
2427 : :
2428 : : // Send mesh and fields data (solution dump) for output to file
2429 [ + - ][ + - ]: 580 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
2430 : 290 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
2431 : : {}, nodefields, {}, nodesurfs, cb );
2432 [ + - ][ + - ]: 870 : }
[ + - ][ + - ]
[ + - ][ - - ]
[ - - ]
2433 : :
2434 : : void
2435 : 3662 : ChoCG::out()
2436 : : // *****************************************************************************
2437 : : // Output mesh field data
2438 : : // *****************************************************************************
2439 : : {
2440 : 3662 : auto d = Disc();
2441 : :
2442 : : // Time history
2443 [ + - ][ + - ]: 3662 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ][ - + ]
2444 : 0 : auto ncomp = m_u.nprop();
2445 : 0 : const auto& inpoel = d->Inpoel();
2446 [ - - ]: 0 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
2447 : 0 : std::size_t j = 0;
2448 [ - - ]: 0 : for (const auto& p : d->Hist()) {
2449 : 0 : auto e = p.get< tag::elem >(); // host element id
2450 : 0 : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
2451 [ - - ]: 0 : hist[j].resize( ncomp+1, 0.0 );
2452 [ - - ]: 0 : for (std::size_t i=0; i<4; ++i) {
2453 [ - - ]: 0 : const auto u = m_u[ inpoel[e*4+i] ];
2454 : 0 : hist[j][0] += n[i] * u[0];
2455 : 0 : hist[j][1] += n[i] * u[1]/u[0];
2456 : 0 : hist[j][2] += n[i] * u[2]/u[0];
2457 : 0 : hist[j][3] += n[i] * u[3]/u[0];
2458 : 0 : hist[j][4] += n[i] * u[4]/u[0];
2459 : 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];
2460 : 0 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
2461 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
2462 : 0 : }
2463 : 0 : ++j;
2464 : : }
2465 [ - - ]: 0 : d->history( std::move(hist) );
2466 : 0 : }
2467 : :
2468 : : // Field data
2469 [ + + ][ + - ]: 3662 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
[ + + ]
2470 [ + - ][ + - ]: 509 : writeFields( CkCallback(CkIndex_ChoCG::integrals(), thisProxy[thisIndex]) );
[ + - ][ + - ]
2471 : : } else {
2472 : 3153 : integrals();
2473 : : }
2474 : 3662 : }
2475 : :
2476 : : void
2477 : 3662 : ChoCG::integrals()
2478 : : // *****************************************************************************
2479 : : // Compute integral quantities for output
2480 : : // *****************************************************************************
2481 : : {
2482 : 3662 : auto d = Disc();
2483 : :
2484 [ + - ][ + - ]: 3662 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ][ - + ]
2485 : :
2486 : : using namespace integrals;
2487 [ - - ]: 0 : std::vector< std::map< int, tk::real > > ints( NUMINT );
2488 : : // Prepend integral vector with metadata on the current time step:
2489 : : // current iteration count, current physical time, time step size
2490 [ - - ]: 0 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
2491 [ - - ]: 0 : ints[ TIME ][ 0 ] = d->T();
2492 [ - - ]: 0 : ints[ DT ][ 0 ] = d->Dt();
2493 : :
2494 : : // Compute integrals requested for surfaces requested
2495 : 0 : const auto& reqv = g_cfg.get< tag::integout_integrals >();
2496 [ - - ]: 0 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
2497 [ - - ][ - - ]: 0 : if (req.count("mass_flow_rate")) {
[ - - ]
2498 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
2499 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
2500 : 0 : const auto& nodes = sint.first;
2501 : 0 : const auto& ndA = sint.second;
2502 : 0 : auto n = ndA.data();
2503 [ - - ]: 0 : for (auto p : nodes) {
2504 [ - - ][ - - ]: 0 : mfr += n[0]*m_u(p,0) + n[1]*m_u(p,1) + n[2]*m_u(p,2);
[ - - ]
2505 : 0 : n += 3;
2506 : : }
2507 : : }
2508 : : }
2509 [ - - ][ - - ]: 0 : if (req.count("force")) {
[ - - ]
2510 : 0 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
2511 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
2512 [ - - ]: 0 : auto& fx = ints[ FORCE_X ][ s ];
2513 [ - - ]: 0 : auto& fy = ints[ FORCE_Y ][ s ];
2514 [ - - ]: 0 : auto& fz = ints[ FORCE_Z ][ s ];
2515 : 0 : const auto& nodes = sint.first;
2516 : 0 : const auto& ndA = sint.second;
2517 : 0 : auto n = ndA.data();
2518 [ - - ]: 0 : for (auto p : nodes) {
2519 : : // pressure force
2520 : 0 : fx -= n[0]*m_pr[p];
2521 : 0 : fy -= n[1]*m_pr[p];
2522 : 0 : fz -= n[2]*m_pr[p];
2523 : : // viscous force
2524 [ - - ][ - - ]: 0 : fx += mu*(m_vgrad(p,0)*n[0] + m_vgrad(p,1)*n[1] + m_vgrad(p,2)*n[2]);
[ - - ]
2525 [ - - ][ - - ]: 0 : fy += mu*(m_vgrad(p,3)*n[0] + m_vgrad(p,4)*n[1] + m_vgrad(p,5)*n[2]);
[ - - ]
2526 [ - - ][ - - ]: 0 : fz += mu*(m_vgrad(p,6)*n[0] + m_vgrad(p,7)*n[1] + m_vgrad(p,8)*n[2]);
[ - - ]
2527 : 0 : n += 3;
2528 : : }
2529 : : }
2530 : : }
2531 : :
2532 [ - - ]: 0 : auto stream = serialize( d->MeshId(), ints );
2533 [ - - ]: 0 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
2534 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
2535 : :
2536 : 0 : } else {
2537 : :
2538 : 3662 : step();
2539 : :
2540 : : }
2541 : 3662 : }
2542 : :
2543 : : void
2544 : 3297 : ChoCG::evalLB( int nrestart )
2545 : : // *****************************************************************************
2546 : : // Evaluate whether to do load balancing
2547 : : //! \param[in] nrestart Number of times restarted
2548 : : // *****************************************************************************
2549 : : {
2550 : 3297 : auto d = Disc();
2551 : :
2552 : : // Detect if just returned from a checkpoint and if so, zero timers and
2553 : : // finished flag
2554 [ + + ]: 3297 : if (d->restarted( nrestart )) m_finished = 0;
2555 : :
2556 : 3297 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
2557 : 3297 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
2558 : :
2559 : : // Load balancing if user frequency is reached or after the second time-step
2560 [ + + ][ + + ]: 3297 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
[ + + ]
2561 : :
2562 : 3157 : AtSync();
2563 [ - + ]: 3157 : if (nonblocking) dt();
2564 : :
2565 : : } else {
2566 : :
2567 : 140 : dt();
2568 : :
2569 : : }
2570 : 3297 : }
2571 : :
2572 : : void
2573 : 3292 : ChoCG::evalRestart()
2574 : : // *****************************************************************************
2575 : : // Evaluate whether to save checkpoint/restart
2576 : : // *****************************************************************************
2577 : : {
2578 : 3292 : auto d = Disc();
2579 : :
2580 : 3292 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
2581 : 3292 : const auto benchmark = g_cfg.get< tag::benchmark >();
2582 : :
2583 [ + + ][ - + ]: 3292 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
2584 : :
2585 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
2586 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
2587 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
2588 : :
2589 : 0 : } else {
2590 : :
2591 : 3292 : evalLB( /* nrestart = */ -1 );
2592 : :
2593 : : }
2594 : 3292 : }
2595 : :
2596 : : void
2597 : 3662 : ChoCG::step()
2598 : : // *****************************************************************************
2599 : : // Evaluate whether to continue with next time step
2600 : : // *****************************************************************************
2601 : : {
2602 : 3662 : auto d = Disc();
2603 : :
2604 : : // Output one-liner status report to screen
2605 [ + + ]: 3662 : if(thisIndex == 0) d->status();
2606 : :
2607 [ + + ]: 3662 : if (not m_finished) {
2608 : :
2609 : 3292 : evalRestart();
2610 : :
2611 : : } else {
2612 : :
2613 : 370 : auto meshid = d->MeshId();
2614 [ + - ]: 370 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2615 [ + - ][ + - ]: 740 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
2616 : :
2617 : : }
2618 : 3662 : }
2619 : :
2620 : : #include "NoWarning/chocg.def.h"
|
