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-2025 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 : 366 : 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 : 366 : const std::vector< std::size_t >& triinpoel ) :
61 : : m_disc( disc ),
62 [ + - ]: 366 : m_cgpre( cgpre ),
63 [ + - ]: 366 : m_cgmom( cgmom ),
64 : 366 : m_nrhs( 0 ),
65 : 366 : m_nnorm( 0 ),
66 : 366 : m_naec( 0 ),
67 : 366 : m_nalw( 0 ),
68 : 366 : m_nlim( 0 ),
69 : 366 : m_nsgrad( 0 ),
70 : 366 : m_npgrad( 0 ),
71 : 366 : m_nvgrad( 0 ),
72 : 366 : m_nflux( 0 ),
73 : 366 : m_ndiv( 0 ),
74 : 366 : m_nbpint( 0 ),
75 [ + - ]: 366 : m_np( 0 ),
76 : : m_bnode( bnode ),
77 : : m_bface( bface ),
78 [ + - ][ + - ]: 366 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
79 [ + - ]: 366 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
80 : : m_un( m_u.nunk(), m_u.nprop() ),
81 [ + - ][ + - ]: 366 : m_pr( m_u.nunk(), 0.0 ),
[ - - ]
82 [ + - ]: 366 : m_p( m_u.nunk(), m_u.nprop()*2 ),
83 [ + - ]: 366 : m_q( m_u.nunk(), m_u.nprop()*2 ),
84 : : m_a( m_u.nunk(), m_u.nprop() ),
85 : : m_rhs( m_u.nunk(), m_u.nprop() ),
86 : : m_sgrad( m_u.nunk(), 3UL ),
87 : : m_pgrad( m_u.nunk(), 3UL ),
88 : : m_vgrad( m_u.nunk(), 9UL ),
89 : : m_flux( m_u.nunk(), 3UL ),
90 [ + - ][ + - ]: 732 : m_div( m_u.nunk() ),
91 : 366 : m_stage( 0 ),
92 : 366 : m_finished( 0 ),
93 [ + - ][ + - ]: 732 : 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 : 366 : usesAtSync = true; // enable migration at AtSync
105 : :
106 [ + - ]: 366 : auto d = Disc();
107 : :
108 : : // Create new local ids based on mesh locality
109 : : std::unordered_map< std::size_t, std::size_t > map;
110 : : std::size_t n = 0;
111 : :
112 [ + - ][ + - ]: 366 : auto psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
113 [ + + ]: 26635 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
114 [ + - ]: 1076 : if (!map.count(p)) map[p] = n++;
115 [ + + ][ + + ]: 252769 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
116 [ + - ]: 25193 : if (!map.count(q)) map[q] = n++;
117 : : }
118 : : }
119 : :
120 : : Assert( map.size() == d->Gid().size(),
121 : : "Mesh-locality reorder map size mismatch" );
122 : :
123 : : // Remap data in bound Discretization object
124 [ + - ]: 366 : d->remap( map );
125 : : // Remap boundary triangle face connectivity
126 [ + - ]: 366 : tk::remap( m_triinpoel, map );
127 : : // Recompute points surrounding points
128 [ + - ][ + - ]: 366 : psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
129 : :
130 : : // Compute total box IC volume
131 [ + - ]: 366 : d->boxvol();
132 : :
133 : : // Setup LHS matrix for pressure solve
134 [ + - ][ + - ]: 1098 : 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 [ + + ]: 366 : if (g_cfg.get< tag::theta >() > std::numeric_limits< tk::real >::epsilon()) {
141 [ + - ][ + - ]: 4 : m_cgmom[ thisIndex ].insert( momlhs( psup ),
[ + - ][ - + ]
[ - - ]
142 : : d->Gid(),
143 : : d->Lid(),
144 : : d->NodeCommMap() );
145 : : }
146 : :
147 : : // Activate SDAG waits for setup
148 [ + - ][ + - ]: 366 : thisProxy[ thisIndex ].wait4int();
149 : 732 : }
150 : :
151 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
152 : 366 : 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 : 366 : auto d = Disc();
161 : : const auto& inpoel = d->Inpoel();
162 : : const auto& coord = d->Coord();
163 : : const auto& X = coord[0];
164 : : const auto& Y = coord[1];
165 : : const auto& Z = coord[2];
166 : :
167 : : // Matrix with compressed sparse row storage
168 : 366 : tk::CSR A( /*DOF=*/ 1, psup );
169 : :
170 : : // Fill matrix with Laplacian
171 [ + + ]: 63532 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
172 : 63166 : const auto N = inpoel.data() + e*4;
173 : : const std::array< tk::real, 3 >
174 : 63166 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
175 : 63166 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
176 : 63166 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
177 : 63166 : const auto J = tk::triple( ba, ca, da ) * 6.0;
178 : : std::array< std::array< tk::real, 3 >, 4 > grad;
179 : 63166 : grad[1] = tk::cross( ca, da );
180 : 63166 : grad[2] = tk::cross( da, ba );
181 : 63166 : grad[3] = tk::cross( ba, ca );
182 [ + + ]: 252664 : for (std::size_t i=0; i<3; ++i)
183 : 189498 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
184 [ + + ]: 315830 : for (std::size_t a=0; a<4; ++a)
185 [ + + ]: 1263320 : for (std::size_t b=0; b<4; ++b)
186 [ + - ]: 1010656 : A(N[a],N[b]) -= tk::dot( grad[a], grad[b] ) / J;
187 : : }
188 : :
189 : : auto nunk = X.size();
190 [ + - ][ + - ]: 366 : std::vector< tk::real > x( nunk, 0.0 ), b( nunk, 0.0 );
[ - - ]
191 : :
192 : 366 : return { std::move(A), std::move(x), std::move(b) };
193 : 366 : }
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 : : 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 : 2 : return { std::move(A), std::move(x), std::move(b) };
213 : 2 : }
214 : :
215 : : void
216 : 732 : 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 : : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
232 [ + + ]: 1057 : for (const auto& s : cfgmask) {
233 : : auto k = m_bface.find(s[0]);
234 [ + + ]: 325 : if (k != end(m_bface)) {
235 [ + - ]: 177 : auto& n = dir[ k->first ];
236 [ + + ]: 23993 : for (auto f : k->second) {
237 [ + - ]: 23816 : n.insert( m_triinpoel[f*3+0] );
238 [ + - ]: 23816 : n.insert( m_triinpoel[f*3+1] );
239 [ + - ]: 23816 : n.insert( m_triinpoel[f*3+2] );
240 : : }
241 : : }
242 : : }
243 : :
244 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
245 [ + - ]: 732 : const auto& lid = Disc()->Lid();
246 [ + + ]: 1057 : for (const auto& s : cfgmask) {
247 : : auto k = m_bnode.find(s[0]);
248 [ + + ]: 325 : if (k != end(m_bnode)) {
249 [ + - ]: 178 : auto& n = dir[ k->first ];
250 [ + - ][ + + ]: 17940 : for (auto g : k->second) {
251 : : n.insert( tk::cref_find(lid,g) );
252 : : }
253 : : }
254 : : }
255 : :
256 : : // Associate sidesets to Dirichlet BC values if configured by user
257 : : std::unordered_map< int, std::vector< double > > dirval;
258 [ + + ]: 798 : 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 : 732 : auto nmask = ncomp + 1;
269 : : std::unordered_map< std::size_t,
270 : : std::pair< std::vector< int >,
271 : : std::vector< double > > > dirbcset;
272 [ + + ]: 1057 : for (const auto& vec : cfgmask) {
273 [ - + ][ - - ]: 325 : ErrChk( vec.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
274 : : auto n = dir.find( vec[0] );
275 [ + + ]: 325 : if (n != end(dir)) {
276 [ + - ][ + + ]: 178 : std::vector< double > v( ncomp, 0.0 );
277 : : auto m = dirval.find( vec[0] );
278 [ + + ]: 178 : if (m != end(dirval)) {
279 [ - + ][ - - ]: 24 : ErrChk( m->second.size() == ncomp, "Incorrect Dirichlet BC val ncomp" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
280 [ + - ]: 24 : v = m->second;
281 : : }
282 [ + - ][ + + ]: 15562 : for (auto p : n->second) {
283 : : auto& mv = dirbcset[p]; // mask & value
284 [ + - ]: 15384 : mv.second = v;
285 [ + + ]: 15384 : auto& mval = mv.first;
286 [ + + ][ + - ]: 15384 : if (mval.empty()) mval.resize( ncomp, 0 );
[ - - ]
287 [ + + ]: 56580 : for (std::size_t c=0; c<ncomp; ++c)
288 [ + + ]: 41196 : if (!mval[c]) mval[c] = vec[c+1]; // overwrite mask if 0 -> 1
289 : : }
290 : : }
291 : : }
292 : :
293 : : // Compile streamable list of nodes + Dirichlet BC components mask and values
294 : : tk::destroy( mask );
295 [ + + ]: 15536 : for (const auto& [p,mv] : dirbcset) {
296 [ + - ]: 14804 : mask.push_back( p );
297 [ + - ]: 14804 : mask.insert( end(mask), begin(mv.first), end(mv.first) );
298 [ + - ][ + - ]: 14804 : val.push_back( static_cast< double >( p ) );
299 [ + - ]: 14804 : val.insert( end(val), begin(mv.second), end(mv.second) );
300 : : }
301 : :
302 [ - + ][ - - ]: 732 : ErrChk( mask.size() % nmask == 0, "Dirichlet BC mask incomplete" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
303 [ - + ][ - - ]: 732 : ErrChk( val.size() % nmask == 0, "Dirichlet BC val incomplete" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
304 [ - + ][ - - ]: 732 : ErrChk( mask.size() == val.size(), "Dirichlet BC mask & val size mismatch" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
305 : 732 : }
306 : :
307 : : void
308 : 366 : ChoCG::feop()
309 : : // *****************************************************************************
310 : : // Start (re-)computing finite element domain and boundary operators
311 : : // *****************************************************************************
312 : : {
313 : 366 : auto d = Disc();
314 : :
315 : : // Prepare Dirichlet boundary conditions data structures
316 : 366 : setupDirBC( g_cfg.get< tag::bc_dir >(), g_cfg.get< tag::bc_dirval >(),
317 : 366 : m_u.nprop(), m_dirbcmask, m_dirbcval );
318 : 366 : setupDirBC( g_cfg.get< tag::pre_bc_dir >(), g_cfg.get< tag::pre_bc_dirval >(),
319 : 366 : 1, m_dirbcmaskp, m_dirbcvalp );
320 : :
321 : : // Compute local contributions to boundary normals and integrals
322 : 366 : bndint();
323 : : // Compute local contributions to domain edge integrals
324 : 366 : domint();
325 : :
326 : : // Send boundary point normal contributions to neighbor chares
327 [ + + ]: 366 : if (d->NodeCommMap().empty()) {
328 : 13 : comnorm_complete();
329 : : } else {
330 [ + + ]: 4091 : for (const auto& [c,nodes] : d->NodeCommMap()) {
331 : : decltype(m_bnorm) exp;
332 [ + + ]: 19274 : for (auto i : nodes) {
333 [ + + ]: 48170 : for (const auto& [s,b] : m_bnorm) {
334 : : auto k = b.find(i);
335 [ + + ]: 40557 : if (k != end(b)) exp[s][i] = k->second;
336 : : }
337 : : }
338 [ + - ][ + - ]: 7476 : thisProxy[c].comnorm( exp );
339 : : }
340 : : }
341 : 366 : ownnorm_complete();
342 : 366 : }
343 : :
344 : : void
345 : 366 : ChoCG::bndint()
346 : : // *****************************************************************************
347 : : // Compute local contributions to boundary normals and integrals
348 : : // *****************************************************************************
349 : : {
350 : 366 : auto d = Disc();
351 : : const auto& coord = d->Coord();
352 : : const auto& gid = d->Gid();
353 : : const auto& x = coord[0];
354 : : const auto& y = coord[1];
355 : : 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 : 118050 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
361 : 118050 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
362 : 118050 : (c[1] - y[p]) * (c[1] - y[p]) +
363 : 118050 : (c[2] - z[p]) * (c[2] - z[p]) );
364 : 366 : };
365 : :
366 : 366 : tk::destroy( m_bnorm );
367 : 366 : tk::destroy( m_bndpoinint );
368 : :
369 [ + + ]: 1240 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
370 [ + + ]: 40224 : for (auto f : faceids) { // for all side set triangles
371 : 39350 : const auto N = m_triinpoel.data() + f*3;
372 : : const std::array< tk::real, 3 >
373 : 39350 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
374 : 39350 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
375 : : auto n = tk::cross( ba, ca );
376 : : auto A2 = tk::length( n );
377 : 39350 : n[0] /= A2;
378 : 39350 : n[1] /= A2;
379 : 39350 : n[2] /= A2;
380 : : const tk::real centroid[3] = {
381 : 39350 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
382 : 39350 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
383 : 39350 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
384 [ + + ]: 157400 : for (const auto& [i,j] : tk::lpoet) {
385 : 118050 : auto p = N[i];
386 : 118050 : tk::real r = invdistsq( centroid, p );
387 : : auto& v = m_bnorm[setid]; // associate side set id
388 : : auto& bpn = v[gid[p]]; // associate global node id of bnd pnt
389 : 118050 : bpn[0] += r * n[0]; // inv.dist.sq-weighted normal
390 : 118050 : bpn[1] += r * n[1];
391 : 118050 : bpn[2] += r * n[2];
392 : 118050 : bpn[3] += r; // inv.dist.sq of node from centroid
393 : : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
394 : 118050 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
395 : 118050 : b[1] += n[1] * A2 / 6.0;
396 : 118050 : b[2] += n[2] * A2 / 6.0;
397 : : }
398 : : }
399 : : }
400 : 366 : }
401 : :
402 : : void
403 : 366 : ChoCG::domint()
404 : : // *****************************************************************************
405 : : //! Compute local contributions to domain edge integrals
406 : : // *****************************************************************************
407 : : {
408 : 366 : auto d = Disc();
409 : :
410 : : const auto& gid = d->Gid();
411 : : const auto& inpoel = d->Inpoel();
412 : :
413 : : const auto& coord = d->Coord();
414 : : const auto& x = coord[0];
415 : : const auto& y = coord[1];
416 : : const auto& z = coord[2];
417 : :
418 : 366 : tk::destroy( m_domedgeint );
419 : :
420 [ + + ]: 63532 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
421 : 63166 : const auto N = inpoel.data() + e*4;
422 : : const std::array< tk::real, 3 >
423 : 63166 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
424 : 63166 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
425 : 63166 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
426 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
427 : : Assert( J > 0, "Element Jacobian non-positive" );
428 : : std::array< std::array< tk::real, 3 >, 4 > grad;
429 : 63166 : grad[1] = tk::cross( ca, da );
430 : 63166 : grad[2] = tk::cross( da, ba );
431 : 63166 : grad[3] = tk::cross( ba, ca );
432 [ + + ]: 252664 : for (std::size_t i=0; i<3; ++i)
433 : 189498 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
434 [ + + ]: 442162 : for (const auto& [p,q] : tk::lpoed) {
435 [ + + ]: 378996 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
436 : : tk::real sig = 1.0;
437 [ + + ]: 378996 : if (ed[0] > ed[1]) {
438 : : std::swap( ed[0], ed[1] );
439 : : sig = -1.0;
440 : : }
441 : : auto& n = m_domedgeint[ ed ];
442 : 378996 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
443 : 378996 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
444 : 378996 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
445 : 378996 : n[3] += J / 120.0;
446 : 378996 : n[4] += tk::dot( grad[p], grad[q] ) / J / 6.0;
447 : : }
448 : : }
449 : 366 : }
450 : :
451 : : void
452 : 3738 : 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 [ + + ]: 6612 : for (const auto& [s,b] : inbnd) {
460 : : auto& bndnorm = m_bnormc[s];
461 [ + + ]: 10797 : for (const auto& [p,n] : b) {
462 : : auto& norm = bndnorm[p];
463 : 7923 : norm[0] += n[0];
464 : 7923 : norm[1] += n[1];
465 : 7923 : norm[2] += n[2];
466 : 7923 : norm[3] += n[3];
467 : : }
468 : : }
469 : :
470 [ + + ]: 3738 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
471 : 353 : m_nnorm = 0;
472 : 353 : comnorm_complete();
473 : : }
474 : 3738 : }
475 : :
476 : : void
477 : 784 : 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 : 784 : NodeDiagnostics::registerReducers();
488 : 784 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
489 : 784 : }
490 : :
491 : : void
492 : : // cppcheck-suppress unusedFunction
493 : 3176 : 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 [ - + ][ - - ]: 3176 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
501 : :
502 [ + - ]: 3176 : if (!g_cfg.get< tag::nonblocking >()) dt();
503 : 3176 : }
504 : :
505 : : void
506 : 366 : ChoCG::setup( tk::real v )
507 : : // *****************************************************************************
508 : : // Start setup for solution
509 : : //! \param[in] v Total volume within user-specified box
510 : : // *****************************************************************************
511 : : {
512 : 366 : auto d = Disc();
513 : :
514 : : // Store user-defined box IC volume
515 : 366 : Disc()->Boxvol() = v;
516 : :
517 : : // Set initial conditions
518 : 366 : problems::initialize( d->Coord(), m_u, d->T(), d->BoxNodes() );
519 : :
520 : : // Query time history field output labels from all PDEs integrated
521 [ - + ]: 366 : if (!g_cfg.get< tag::histout >().empty()) {
522 : : std::vector< std::string > var
523 [ - - ][ - - ]: 0 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ][ - - ]
[ - - ]
524 : : 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 : 366 : feop();
532 [ - - ][ - - ]: 366 : }
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
533 : :
534 : : void
535 : 366 : ChoCG::bnorm()
536 : : // *****************************************************************************
537 : : // Combine own and communicated portions of the boundary point normals
538 : : // *****************************************************************************
539 : : {
540 : 366 : const auto& lid = Disc()->Lid();
541 : :
542 : : // Combine own and communicated contributions to boundary point normals
543 [ + + ]: 1203 : for (const auto& [s,b] : m_bnormc) {
544 : : auto& bndnorm = m_bnorm[s];
545 [ + + ]: 6552 : for (const auto& [g,n] : b) {
546 : : auto& norm = bndnorm[g];
547 : 5715 : norm[0] += n[0];
548 : 5715 : norm[1] += n[1];
549 : 5715 : norm[2] += n[2];
550 : 5715 : norm[3] += n[3];
551 : : }
552 : : }
553 : 366 : tk::destroy( m_bnormc );
554 : :
555 : : // Divide summed point normals by the sum of the inverse distance squared
556 [ + + ]: 1295 : for (auto& [s,b] : m_bnorm) {
557 [ + + ]: 30126 : for (auto& [g,n] : b) {
558 : 29197 : n[0] /= n[3];
559 : 29197 : n[1] /= n[3];
560 : 29197 : n[2] /= n[3];
561 : : 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 : : decltype(m_bnorm) loc;
569 [ + + ]: 1295 : for (auto& [s,b] : m_bnorm) {
570 : : auto& bnd = loc[s];
571 [ + + ]: 30126 : for (auto&& [g,n] : b) {
572 : 29197 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
573 : : }
574 : : }
575 : : m_bnorm = std::move(loc);
576 : 366 : }
577 : :
578 : : void
579 [ + - ]: 366 : ChoCG::streamable()
580 : : // *****************************************************************************
581 : : // Convert integrals into streamable data structures
582 : : // *****************************************************************************
583 : : {
584 : : // Query surface integral output nodes
585 : : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
586 : : const auto& is = g_cfg.get< tag::integout >();
587 [ + - ]: 366 : std::set< int > outsets( begin(is), end(is) );
588 [ + + ]: 367 : for (auto s : outsets) {
589 : : auto m = m_bface.find(s);
590 [ + - ]: 1 : if (m != end(m_bface)) {
591 [ + - ]: 1 : auto& n = surfintnodes[ m->first ]; // associate set id
592 [ + - ][ + + ]: 847 : for (auto f : m->second) { // face ids on side set
593 : 846 : auto t = m_triinpoel.data() + f*3;
594 [ + - ]: 846 : n.push_back( t[0] ); // nodes on side set
595 [ + - ]: 846 : n.push_back( t[1] );
596 [ + - ]: 846 : n.push_back( t[2] );
597 : : }
598 : : }
599 : : }
600 [ + - ][ + + ]: 367 : for (auto& [s,n] : surfintnodes) tk::unique( n );
601 : : // Prepare surface integral data
602 : 366 : tk::destroy( m_surfint );
603 [ + - ]: 366 : const auto& gid = Disc()->Gid();
604 [ + + ]: 367 : for (auto&& [s,n] : surfintnodes) {
605 [ + - ]: 1 : auto& sint = m_surfint[s]; // associate set id
606 : 1 : auto& nodes = sint.first;
607 : 1 : auto& ndA = sint.second;
608 : : nodes = std::move(n);
609 [ + - ]: 1 : ndA.resize( nodes.size()*3 );
610 : : auto a = ndA.data();
611 [ + + ]: 426 : for (auto p : nodes) {
612 : : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
613 : 425 : a[0] = b[0]; // store ni * dA
614 : 425 : a[1] = b[1];
615 : 425 : a[2] = b[2];
616 : 425 : a += 3;
617 : : }
618 : : }
619 : 366 : tk::destroy( m_bndpoinint );
620 : :
621 : : // Generate domain superedges
622 [ + - ]: 366 : domsuped();
623 : :
624 : : // Prepare symmetry boundary condition data structures
625 : :
626 : : // Query symmetry BC nodes associated to side sets
627 : : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
628 [ + + ]: 668 : for (auto s : g_cfg.get< tag::bc_sym >()) {
629 : : auto k = m_bface.find(s);
630 [ + + ]: 302 : if (k != end(m_bface)) {
631 [ + - ]: 117 : auto& n = sym[ k->first ];
632 [ + - ][ + + ]: 3561 : for (auto f : k->second) {
633 [ + - ]: 3444 : const auto& t = m_triinpoel.data() + f*3;
634 : : 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 : : std::set< std::size_t > symbcnodeset;
643 [ + + ]: 483 : for (const auto& [s,n] : sym) symbcnodeset.insert( begin(n), end(n) );
644 : :
645 : : // Generate symmetry BC data as streamable data structures
646 [ - + ]: 366 : tk::destroy( m_symbcnodes );
647 [ - + ]: 366 : tk::destroy( m_symbcnorms );
648 [ + + ]: 2803 : for (auto p : symbcnodeset) {
649 [ + + ]: 6093 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
650 : : auto m = m_bnorm.find( s );
651 [ + + ]: 3656 : if (m != end(m_bnorm)) {
652 : : auto r = m->second.find( p );
653 [ + + ]: 3138 : if (r != end(m->second)) {
654 [ + - ]: 2437 : m_symbcnodes.push_back( p );
655 [ + - ]: 2437 : m_symbcnorms.push_back( r->second[0] );
656 [ + - ]: 2437 : m_symbcnorms.push_back( r->second[1] );
657 [ + - ]: 2437 : m_symbcnorms.push_back( r->second[2] );
658 : : }
659 : : }
660 : : }
661 : : }
662 : 366 : tk::destroy( m_bnorm );
663 : :
664 : : // Prepare noslip boundary condition data structures
665 : :
666 : : // Query noslip BC nodes associated to side sets
667 : : std::unordered_map< int, std::unordered_set< std::size_t > > noslip;
668 [ + + ]: 457 : for (auto s : g_cfg.get< tag::bc_noslip >()) {
669 : : auto k = m_bface.find(s);
670 [ + + ]: 91 : if (k != end(m_bface)) {
671 [ + - ]: 85 : auto& n = noslip[ k->first ];
672 [ + - ][ + + ]: 5951 : for (auto f : k->second) {
673 [ + - ]: 5866 : const auto& t = m_triinpoel.data() + f*3;
674 : : n.insert( t[0] );
675 [ + - ]: 5866 : n.insert( t[1] );
676 [ + - ]: 5866 : n.insert( t[2] );
677 : : }
678 : : }
679 : : }
680 : :
681 : : // Generate unique set of noslip BC nodes of all noslip BC side sets
682 : : std::set< std::size_t > noslipbcnodeset;
683 [ + + ]: 451 : for (const auto& [s,n] : noslip) noslipbcnodeset.insert( begin(n), end(n) );
684 : :
685 : : // Generate noslip BC data as streamable data structures
686 [ - + ]: 366 : tk::destroy( m_noslipbcnodes );
687 [ + - ]: 366 : m_noslipbcnodes.insert( m_noslipbcnodes.end(),
688 : : begin(noslipbcnodeset), end(noslipbcnodeset) );
689 : 366 : }
690 : :
691 : : void
692 : 366 : 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 : : Assert( !m_domedgeint.empty(), "No domain edges to group" );
700 : :
701 : : #ifndef NDEBUG
702 : : auto nedge = m_domedgeint.size();
703 : : #endif
704 : :
705 : 366 : const auto& inpoel = Disc()->Inpoel();
706 : 366 : const auto& lid = Disc()->Lid();
707 : 366 : const auto& gid = Disc()->Gid();
708 : :
709 : : tk::destroy( m_dsupedge[0] );
710 : : tk::destroy( m_dsupedge[1] );
711 : : tk::destroy( m_dsupedge[2] );
712 : :
713 : : tk::destroy( m_dsupint[0] );
714 : : tk::destroy( m_dsupint[1] );
715 : : tk::destroy( m_dsupint[2] );
716 : :
717 : : tk::UnsMesh::FaceSet untri;
718 [ + + ]: 63532 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
719 : : std::size_t N[4] = {
720 : 63166 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
721 [ + - ][ + + ]: 315830 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
722 : : }
723 : :
724 [ + + ]: 63532 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
725 : : std::size_t N[4] = {
726 : 63166 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
727 : : int f = 0;
728 : : tk::real sig[6];
729 [ + + ]: 442162 : decltype(m_domedgeint)::const_iterator d[6];
730 [ + + ]: 183152 : for (const auto& [p,q] : tk::lpoed) {
731 [ + + ]: 173128 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
732 [ + + ]: 227818 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
733 : 173128 : d[f] = m_domedgeint.find( ed );
734 [ + + ]: 173128 : if (d[f] == end(m_domedgeint)) break; else ++f;
735 : : }
736 [ + + ]: 63166 : if (f == 6) {
737 [ + - ]: 10024 : m_dsupedge[0].push_back( N[0] );
738 [ + - ]: 10024 : m_dsupedge[0].push_back( N[1] );
739 [ + - ]: 10024 : m_dsupedge[0].push_back( N[2] );
740 [ + - ]: 10024 : m_dsupedge[0].push_back( N[3] );
741 [ + + ]: 50120 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
742 [ + + ]: 70168 : for (int ed=0; ed<6; ++ed) {
743 : : const auto& ded = d[ed]->second;
744 [ + - ]: 60144 : m_dsupint[0].push_back( sig[ed] * ded[0] );
745 [ + - ]: 60144 : m_dsupint[0].push_back( sig[ed] * ded[1] );
746 [ + - ][ + - ]: 60144 : m_dsupint[0].push_back( sig[ed] * ded[2] );
747 [ + - ]: 60144 : m_dsupint[0].push_back( ded[3] );
748 [ + - ]: 60144 : m_dsupint[0].push_back( ded[4] );
749 : 60144 : m_domedgeint.erase( d[ed] );
750 : : }
751 : : }
752 : : }
753 : :
754 [ + + ]: 110786 : for (const auto& N : untri) {
755 : : int f = 0;
756 : : tk::real sig[3];
757 [ + + ]: 441680 : decltype(m_domedgeint)::const_iterator d[3];
758 [ + + ]: 178546 : for (const auto& [p,q] : tk::lpoet) {
759 [ + + ]: 168961 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
760 [ + + ]: 259411 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
761 : 168961 : d[f] = m_domedgeint.find( ed );
762 [ + + ]: 168961 : if (d[f] == end(m_domedgeint)) break; else ++f;
763 : : }
764 [ + + ]: 110420 : if (f == 3) {
765 [ + - ]: 9585 : m_dsupedge[1].push_back( N[0] );
766 [ + - ]: 9585 : m_dsupedge[1].push_back( N[1] );
767 [ + - ]: 9585 : m_dsupedge[1].push_back( N[2] );
768 [ + + ]: 38340 : for (int ed=0; ed<3; ++ed) {
769 : : const auto& ded = d[ed]->second;
770 [ + - ]: 28755 : m_dsupint[1].push_back( sig[ed] * ded[0] );
771 [ + - ]: 28755 : m_dsupint[1].push_back( sig[ed] * ded[1] );
772 [ + - ][ + - ]: 28755 : m_dsupint[1].push_back( sig[ed] * ded[2] );
773 [ + - ]: 28755 : m_dsupint[1].push_back( ded[3] );
774 [ + - ]: 28755 : m_dsupint[1].push_back( ded[4] );
775 : 28755 : m_domedgeint.erase( d[ed] );
776 : : }
777 : : }
778 : : }
779 : :
780 [ + - ]: 366 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
781 [ + - ]: 366 : m_dsupint[2].resize( m_domedgeint.size()*5 );
782 : : std::size_t k = 0;
783 [ + + ]: 24717 : for (const auto& [ed,d] : m_domedgeint) {
784 : 24351 : auto e = m_dsupedge[2].data() + k*2;
785 : 24351 : e[0] = tk::cref_find( lid, ed[0] );
786 : 24351 : e[1] = tk::cref_find( lid, ed[1] );
787 : 24351 : auto i = m_dsupint[2].data() + k*5;
788 : 24351 : i[0] = d[0];
789 : 24351 : i[1] = d[1];
790 : 24351 : i[2] = d[2];
791 : 24351 : i[3] = d[3];
792 : 24351 : i[4] = d[4];
793 : 24351 : ++k;
794 : : }
795 : :
796 [ + + ]: 366 : if (g_cfg.get< tag::fct >()) {
797 : : 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 : 366 : 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 : : 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 : 366 : }
825 : :
826 : : void
827 : : // cppcheck-suppress unusedFunction
828 : 366 : 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 : 366 : bnorm();
835 : :
836 : : // Convert integrals into streamable data structures
837 : 366 : streamable();
838 : :
839 : : // Enforce boundary conditions on initial conditions
840 : 366 : BC( m_u, Disc()->T() );
841 : :
842 : : // Start measuring initial div-free time
843 : 366 : m_timer.emplace_back();
844 : :
845 : : // Compute initial momentum flux
846 [ + - ]: 366 : thisProxy[ thisIndex ].wait4div();
847 [ + - ]: 366 : thisProxy[ thisIndex ].wait4sgrad();
848 : 366 : div( m_u );
849 : 366 : }
850 : :
851 : : void
852 : 11828 : 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 : 11828 : auto d = Disc();
861 : : const auto lid = d->Lid();
862 : :
863 : : // Combine own and communicated contributions
864 [ + + ]: 229624 : for (const auto& [g,r] : gradc) {
865 : 217796 : auto i = tk::cref_find( lid, g );
866 [ + + ]: 1236692 : for (std::size_t c=0; c<r.size(); ++c) grad(i,c) += r[c];
867 : : }
868 : 11828 : tk::destroy(gradc);
869 : :
870 : : // Divide weak result by nodal volume
871 : : const auto& vol = d->Vol();
872 [ + + ]: 1148311 : for (std::size_t p=0; p<grad.nunk(); ++p) {
873 [ + + ]: 6415760 : for (std::size_t c=0; c<grad.nprop(); ++c) {
874 : 5279277 : grad(p,c) /= vol[p];
875 : : }
876 : : }
877 : 11828 : }
878 : :
879 : : void
880 : 4350 : ChoCG::div( const tk::Fields& u )
881 : : // *****************************************************************************
882 : : // Start computing divergence
883 : : // \para[in] u Vector field whose divergence to compute
884 : : // *****************************************************************************
885 : : {
886 : 4350 : auto d = Disc();
887 : : const auto lid = d->Lid();
888 : :
889 : : // Finalize momentum flux communications if needed
890 [ + + ]: 4350 : if (m_np == 1) {
891 [ + - ]: 334 : fingrad( m_flux, m_fluxc );
892 [ + - ]: 334 : physics::symbc( m_flux, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
893 : : }
894 : :
895 : : // Compute divergence
896 : : std::fill( begin(m_div), end(m_div), 0.0 );
897 : 4350 : chorin::div( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel,
898 [ + - ]: 4350 : d->Dt(), m_pr, m_pgrad, u, m_div, m_np>1 );
899 : :
900 : : // Communicate velocity divergence to other chares on chare-boundary
901 [ + + ]: 4350 : if (d->NodeCommMap().empty()) {
902 [ + - ]: 215 : comdiv_complete();
903 : : } else {
904 [ + + ]: 47345 : for (const auto& [c,n] : d->NodeCommMap()) {
905 : : decltype(m_divc) exp;
906 [ + - ][ + + ]: 220868 : for (auto g : n) exp[g] = m_div[ tk::cref_find(lid,g) ];
907 [ + - ][ + - ]: 86420 : thisProxy[c].comdiv( exp );
908 : : }
909 : : }
910 [ + - ]: 4350 : owndiv_complete();
911 : 4350 : }
912 : :
913 : : void
914 : 43210 : 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 [ + + ]: 220868 : 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 [ + + ]: 43210 : if (++m_ndiv == Disc()->NodeCommMap().size()) {
929 : 4135 : m_ndiv = 0;
930 : 4135 : comdiv_complete();
931 : : }
932 : 43210 : }
933 : :
934 : : void
935 : 3494 : ChoCG::velgrad()
936 : : // *****************************************************************************
937 : : // Start computing velocity gradient
938 : : // *****************************************************************************
939 : : {
940 : 3494 : auto d = Disc();
941 : :
942 : : // Compute momentum flux
943 : : m_vgrad.fill( 0.0 );
944 : 3494 : 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 : : const auto& lid = d->Lid();
948 [ + + ]: 3494 : if (d->NodeCommMap().empty()) {
949 : 232 : comvgrad_complete();
950 : : } else {
951 [ + + ]: 41814 : for (const auto& [c,n] : d->NodeCommMap()) {
952 : : decltype(m_vgradc) exp;
953 [ + - ][ + + ]: 185714 : for (auto g : n) exp[g] = m_vgrad[ tk::cref_find(lid,g) ];
954 [ + - ][ + - ]: 77104 : thisProxy[c].comvgrad( exp );
955 : : }
956 : : }
957 : 3494 : ownvgrad_complete();
958 : 3494 : }
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 [ + + ]: 332876 : 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 : 334 : ChoCG::flux()
985 : : // *****************************************************************************
986 : : // Start computing momentum flux
987 : : // *****************************************************************************
988 : : {
989 : 334 : auto d = Disc();
990 : :
991 : : // Finalize computing velocity gradients
992 : 334 : fingrad( m_vgrad, m_vgradc );
993 : :
994 : : // Compute momentum flux
995 : : m_flux.fill( 0.0 );
996 : 334 : chorin::flux( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad,
997 : 334 : m_flux );
998 : :
999 : : // Communicate velocity divergence to other chares on chare-boundary
1000 : : const auto& lid = d->Lid();
1001 [ + + ]: 334 : if (d->NodeCommMap().empty()) {
1002 : 12 : comflux_complete();
1003 : : } else {
1004 [ + + ]: 3874 : for (const auto& [c,n] : d->NodeCommMap()) {
1005 : : decltype(m_fluxc) exp;
1006 [ + - ][ + + ]: 17854 : for (auto g : n) exp[g] = m_flux[ tk::cref_find(lid,g) ];
1007 [ + - ][ + - ]: 7104 : thisProxy[c].comflux( exp );
1008 : : }
1009 : : }
1010 : 334 : ownflux_complete();
1011 : 334 : }
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 [ + + ]: 32156 : 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 : 4350 : ChoCG::pinit()
1038 : : // *****************************************************************************
1039 : : // Initialize Poisson solve
1040 : : // *****************************************************************************
1041 : : {
1042 : 4350 : auto d = Disc();
1043 : : const auto lid = d->Lid();
1044 : : const auto& coord = d->Coord();
1045 : : const auto& x = coord[0];
1046 : : const auto& y = coord[1];
1047 : : const auto& z = coord[2];
1048 : :
1049 : : // Combine own and communicated contributions to velocity divergence
1050 [ + + ]: 86380 : for (const auto& [g,r] : m_divc) m_div[ tk::cref_find(lid,g) ] += r;
1051 : 4350 : tk::destroy(m_divc);
1052 : :
1053 [ + + ][ + + ]: 379500 : 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 : : std::vector< std::pair< int, tk::real > > > dirbc;
1058 [ + + ]: 4350 : if (!g_cfg.get< tag::pre_bc_dir >().empty()) {
1059 [ + - ]: 774 : auto ic = problems::PRESSURE_IC();
1060 : : std::size_t nmask = 1 + 1;
1061 : : Assert( m_dirbcmaskp.size() % nmask == 0, "Size mismatch" );
1062 [ + + ]: 9744 : for (std::size_t i=0; i<m_dirbcmaskp.size()/nmask; ++i) {
1063 [ + + ]: 8970 : auto p = m_dirbcmaskp[i*nmask+0]; // local node id
1064 : 8970 : auto mask = m_dirbcmaskp[i*nmask+1];
1065 [ + + ]: 8970 : if (mask == 1) { // mask == 1: IC value
1066 [ + + ][ - + ]: 5620 : auto val = m_np>1 ? 0.0 : ic( x[p], y[p], z[p] );
1067 [ + - ]: 7380 : 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 [ + - ][ - - ]: 10560 : dirbc[p] = {{ { 1, val } }};
1071 : : }
1072 : : }
1073 : : }
1074 : :
1075 : : // Configure Neumann BCs
1076 : : std::vector< tk::real > neubc;
1077 [ + - ]: 4350 : auto pg = problems::PRESSURE_GRAD();
1078 [ + + ]: 4350 : 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 : : 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 : : }
1103 : :
1104 : : // Set hydrostat
1105 : 4350 : auto h = g_cfg.get< tag::pre_hydrostat >();
1106 [ + + ]: 4350 : if (h != std::numeric_limits< uint64_t >::max()) {
1107 : : 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 [ + + ][ - + ]: 42 : auto val = m_np>1 ? 0.0 : ic( x[p], y[p], z[p] );
1112 : : auto& b = dirbc[p];
1113 [ + - ][ + - ]: 72 : if (b.empty()) b = {{ { 1, val }} };
[ - - ]
1114 : : }
1115 : : }
1116 : :
1117 : : // Configure right hand side
1118 [ + - ]: 4350 : auto pr = problems::PRESSURE_RHS();
1119 [ + + ]: 4350 : if (pr) {
1120 : : const auto& vol = d->Vol();
1121 [ + + ]: 2163 : for (std::size_t i=0; i<x.size(); ++i) {
1122 [ - + ]: 4262 : m_div[i] = pr( x[i], y[i], z[i] ) * vol[i];
1123 : : }
1124 : : }
1125 : :
1126 : : // Initialize Poisson solve
1127 : : const auto& pc = g_cfg.get< tag::pre_pc >();
1128 [ + - ][ + - ]: 17400 : m_cgpre[ thisIndex ].ckLocal()->init( {}, m_div, neubc, dirbc, pc,
[ + + ][ - - ]
1129 [ + - ][ + - ]: 13050 : CkCallback( CkIndex_ChoCG::psolve(), thisProxy[thisIndex] ) );
[ - + ][ - - ]
1130 : 4350 : }
1131 : :
1132 : : void
1133 : 4350 : ChoCG::psolve()
1134 : : // *****************************************************************************
1135 : : // Solve Poisson equation
1136 : : // *****************************************************************************
1137 : : {
1138 : 4350 : auto iter = g_cfg.get< tag::pre_iter >();
1139 : 4350 : auto tol = g_cfg.get< tag::pre_tol >();
1140 : 4350 : auto verbose = g_cfg.get< tag::pre_verbose >();
1141 : :
1142 [ + + ]: 4350 : auto c = m_np != 1 ?
1143 : 4016 : CkCallback( CkIndex_ChoCG::sgrad(), thisProxy[thisIndex] ) :
1144 [ + - ][ + - ]: 9034 : CkCallback( CkIndex_ChoCG::psolved(), thisProxy[thisIndex] );
[ + - ][ + - ]
[ + + ][ - - ]
1145 : :
1146 [ + - ][ + - ]: 13050 : m_cgpre[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose, c );
[ - + ][ - - ]
1147 : 4350 : }
1148 : :
1149 : : void
1150 : 4016 : ChoCG::sgrad()
1151 : : // *****************************************************************************
1152 : : // Compute recent conjugate gradients solution gradient
1153 : : // *****************************************************************************
1154 : : {
1155 : 4016 : auto d = Disc();
1156 : :
1157 [ + - ]: 12048 : auto sol = m_cgpre[ thisIndex ].ckLocal()->solution();
1158 : : m_sgrad.fill( 0.0 );
1159 [ + - ]: 4016 : chorin::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, sol, m_sgrad );
1160 : :
1161 : : // Send gradient contributions to neighbor chares
1162 [ + + ]: 4016 : if (d->NodeCommMap().empty()) {
1163 [ + - ]: 203 : comsgrad_complete();
1164 : : } else {
1165 : : const auto& lid = d->Lid();
1166 [ + + ]: 43471 : for (const auto& [c,n] : d->NodeCommMap()) {
1167 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1168 [ + - ][ + + ]: 203014 : for (auto g : n) exp[g] = m_sgrad[ tk::cref_find(lid,g) ];
1169 [ + - ][ + - ]: 79316 : thisProxy[c].comsgrad( exp );
1170 : : }
1171 : : }
1172 [ + - ]: 4016 : ownsgrad_complete();
1173 : 4016 : }
1174 : :
1175 : : void
1176 : 39658 : 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 [ + + ]: 366370 : 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 [ + + ]: 39658 : if (++m_nsgrad == Disc()->NodeCommMap().size()) {
1193 : 3813 : m_nsgrad = 0;
1194 : 3813 : comsgrad_complete();
1195 : : }
1196 : 39658 : }
1197 : :
1198 : : void
1199 : 4350 : ChoCG::psolved()
1200 : : // *****************************************************************************
1201 : : // Continue setup after Poisson solve and gradient computation
1202 : : // *****************************************************************************
1203 : : {
1204 : 4350 : auto d = Disc();
1205 : :
1206 [ + + ][ + - ]: 5400 : if (thisIndex == 0) d->pit( m_cgpre[ thisIndex ].ckLocal()->it() );
1207 : :
1208 [ + + ]: 4350 : if (m_np != 1) {
1209 : : // Finalize gradient communications
1210 : 4016 : fingrad( m_sgrad, m_sgradc );
1211 : : // Project velocity to divergence-free subspace
1212 [ + + ]: 4016 : auto dt = m_np > 1 ? d->Dt() : 1.0;
1213 [ + + ]: 405435 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1214 : 401419 : m_u(i,0) -= dt * m_sgrad(i,0);
1215 : 401419 : m_u(i,1) -= dt * m_sgrad(i,1);
1216 : 401419 : m_u(i,2) -= dt * m_sgrad(i,2);
1217 : : }
1218 : : // Enforce boundary conditions
1219 : 4016 : BC( m_u, d->T() + d->Dt() );
1220 : : }
1221 : :
1222 [ + + ]: 4350 : if (d->Initial()) {
1223 : :
1224 [ + + ]: 700 : if (g_cfg.get< tag::nstep >() == 1) { // test first Poisson solve only
1225 : :
1226 [ + - ]: 64 : m_pr = m_cgpre[ thisIndex ].ckLocal()->solution();
1227 [ + - ]: 32 : thisProxy[ thisIndex ].wait4step();
1228 [ + - ][ + - ]: 96 : writeFields( CkCallback(CkIndex_ChoCG::diag(), thisProxy[thisIndex]) );
1229 : :
1230 : : } else {
1231 : :
1232 [ + + ]: 668 : if (++m_np < 2) {
1233 : : // Compute momentum flux for initial pressure-Poisson rhs
1234 [ + - ]: 334 : thisProxy[ thisIndex ].wait4vgrad();
1235 [ + - ]: 334 : thisProxy[ thisIndex ].wait4flux();
1236 [ + - ]: 334 : thisProxy[ thisIndex ].wait4div();
1237 : 334 : velgrad();
1238 : : } else {
1239 [ + + ]: 334 : if (thisIndex == 0) {
1240 : : tk::Print() << "Initial div-free time: " << m_timer[0].dsec()
1241 : : << " sec\n";
1242 : : }
1243 : : // Assign initial pressure and compute its gradient
1244 [ + - ]: 668 : m_pr = m_cgpre[ thisIndex ].ckLocal()->solution();
1245 : 334 : pgrad();
1246 : : }
1247 : :
1248 : : }
1249 : :
1250 : : } else {
1251 : :
1252 : : // Update pressure and compute its gradient
1253 : : using tk::operator+=;
1254 [ + - ]: 7300 : m_pr += m_cgpre[ thisIndex ].ckLocal()->solution();
1255 : 3650 : pgrad();
1256 : :
1257 : : }
1258 : 4350 : }
1259 : :
1260 : : void
1261 : 3984 : ChoCG::pgrad()
1262 : : // *****************************************************************************
1263 : : // Compute pressure gradient
1264 : : // *****************************************************************************
1265 : : {
1266 : 3984 : auto d = Disc();
1267 : :
1268 [ + - ]: 7968 : thisProxy[ thisIndex ].wait4pgrad();
1269 : :
1270 : : m_pgrad.fill( 0.0 );
1271 : 3984 : chorin::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_pr, m_pgrad );
1272 : :
1273 : : // Send gradient contributions to neighbor chares
1274 [ + + ]: 3984 : if (d->NodeCommMap().empty()) {
1275 : 202 : compgrad_complete();
1276 : : } else {
1277 : : const auto& lid = d->Lid();
1278 [ + + ]: 43254 : for (const auto& [c,n] : d->NodeCommMap()) {
1279 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1280 [ + - ][ + + ]: 201594 : for (auto g : n) exp[g] = m_pgrad[ tk::cref_find(lid,g) ];
1281 [ + - ][ + - ]: 78944 : thisProxy[c].compgrad( exp );
1282 : : }
1283 : : }
1284 : 3984 : ownpgrad_complete();
1285 : 3984 : }
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 [ + + ]: 363716 : 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 : 3984 : ChoCG::finpgrad()
1312 : : // *****************************************************************************
1313 : : // Compute pressure gradient
1314 : : // *****************************************************************************
1315 : : {
1316 : 3984 : auto d = Disc();
1317 : :
1318 : : // Finalize pressure gradient communications
1319 : 3984 : fingrad( m_pgrad, m_pgradc );
1320 : :
1321 [ + + ]: 3984 : if (d->Initial()) {
1322 [ + - ][ + - ]: 1002 : writeFields( CkCallback(CkIndex_ChoCG::start(), thisProxy[thisIndex]) );
1323 : : } else {
1324 : 3650 : diag();
1325 : : }
1326 : 3984 : }
1327 : :
1328 : : void
1329 : 334 : ChoCG::start()
1330 : : // *****************************************************************************
1331 : : // Start time stepping
1332 : : // *****************************************************************************
1333 : : {
1334 : : // Set flag that indicates that we are now during time stepping
1335 : 334 : Disc()->Initial( 0 );
1336 : : // Start timer measuring time stepping wall clock time
1337 : 334 : Disc()->Timer().zero();
1338 : : // Zero grind-timer
1339 : 334 : Disc()->grindZero();
1340 : : // Continue to first time step
1341 : 334 : dt();
1342 : 334 : }
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 : : const auto ncomp = m_u.nprop();
1353 : : const auto& lid = d->Lid();
1354 : :
1355 : : // Antidiffusive contributions: P+/-
1356 : :
1357 : 2920 : auto ctau = g_cfg.get< tag::fctdif >();
1358 : : 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 : : const auto D = m_dsupint[0].data();
1364 : : 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 [ + + ]: 25570 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1381 : 22650 : const auto N = m_dsupedge[1].data() + e*3;
1382 : : const auto D = m_dsupint[1].data();
1383 : : std::size_t i = 0;
1384 [ + + ]: 90600 : for (const auto& [p,q] : tk::lpoet) {
1385 : 67950 : auto dif = D[(e*3+i)*5+3];
1386 [ + + ]: 271800 : for (std::size_t c=0; c<ncomp; ++c) {
1387 [ - + ]: 203850 : auto aec = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1388 : 203850 : auto a = c*2;
1389 : 203850 : auto b = a+1;
1390 [ - + ]: 203850 : if (aec > 0.0) std::swap(a,b);
1391 : 203850 : m_p(N[p],a) -= aec;
1392 : 203850 : m_p(N[q],b) += aec;
1393 : : }
1394 : 67950 : ++i;
1395 : : }
1396 : : }
1397 : :
1398 : : // edges
1399 [ + + ]: 73670 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1400 : 70750 : const auto N = m_dsupedge[2].data() + e*2;
1401 : 70750 : const auto dif = m_dsupint[2][e*5+3];
1402 [ + + ]: 283000 : for (std::size_t c=0; c<ncomp; ++c) {
1403 [ - + ]: 212250 : auto aec = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1404 : 212250 : auto a = c*2;
1405 : 212250 : auto b = a+1;
1406 [ - + ]: 212250 : if (aec > 0.0) std::swap(a,b);
1407 : 212250 : m_p(N[0],a) -= aec;
1408 : 212250 : 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 : : decltype(m_pc) exp;
1432 [ + - ][ + + ]: 165880 : for (auto g : n) exp[g] = m_p[ tk::cref_find(lid,g) ];
1433 [ + - ][ + - ]: 69760 : thisProxy[c].comaec( exp );
1434 : : }
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 [ + + ]: 296880 : 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 : : const auto npoin = m_u.nunk();
1467 : : const auto ncomp = m_u.nprop();
1468 : : const auto& lid = d->Lid();
1469 : : 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 : : 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 [ - + ][ - + ]: 822240 : 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 : : min(m_rhs(N[q],c), m_u(N[q],c)) );
1521 : : }
1522 [ + + ][ + + ]: 461760 : 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 [ + + ][ + + ]: 540900 : 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 [ + + ]: 25570 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1532 : 22650 : const auto N = m_dsupedge[1].data() + e*3;
1533 [ + + ]: 90600 : for (std::size_t c=0; c<ncomp; ++c) {
1534 : 67950 : auto a = c*2;
1535 : 67950 : auto b = a+1;
1536 [ + + ]: 271800 : for (const auto& [p,q] : tk::lpoet) {
1537 : : tk::real alwp, alwn;
1538 [ - + ]: 203850 : 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 [ - + ][ - + ]: 407700 : alwp = max( max(m_rhs(N[p],c), m_u(N[p],c)),
1543 [ - + ]: 203850 : max(m_rhs(N[q],c), m_u(N[q],c)) );
1544 : 203850 : alwn = min( min(m_rhs(N[p],c), m_u(N[p],c)),
1545 : : min(m_rhs(N[q],c), m_u(N[q],c)) );
1546 : : }
1547 [ + + ][ + + ]: 218520 : m_q(N[p],a) = max(m_q(N[p],a), alwp);
1548 : 203850 : m_q(N[p],b) = min(m_q(N[p],b), alwn);
1549 [ + + ][ + + ]: 232590 : m_q(N[q],a) = max(m_q(N[q],a), alwp);
1550 : 203850 : m_q(N[q],b) = min(m_q(N[q],b), alwn);
1551 : : }
1552 : : }
1553 : : }
1554 : :
1555 : : // edges
1556 [ + + ]: 73670 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1557 : 70750 : const auto N = m_dsupedge[2].data() + e*2;
1558 [ + + ]: 283000 : for (std::size_t c=0; c<ncomp; ++c) {
1559 : 212250 : auto a = c*2;
1560 : 212250 : auto b = a+1;
1561 : : tk::real alwp, alwn;
1562 [ - + ]: 212250 : 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 [ - + ][ - + ]: 424500 : alwp = max( max(m_rhs(N[0],c), m_u(N[0],c)),
1567 [ - + ]: 212250 : max(m_rhs(N[1],c), m_u(N[1],c)) );
1568 : 212250 : alwn = min( min(m_rhs(N[0],c), m_u(N[0],c)),
1569 : : min(m_rhs(N[1],c), m_u(N[1],c)) );
1570 : : }
1571 [ + + ][ + + ]: 216900 : m_q(N[0],a) = max(m_q(N[0],a), alwp);
1572 : 212250 : m_q(N[0],b) = min(m_q(N[0],b), alwn);
1573 [ + + ][ + + ]: 214710 : m_q(N[1],a) = max(m_q(N[1],a), alwp);
1574 : 212250 : 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 : : decltype(m_qc) exp;
1584 [ + - ][ + + ]: 165880 : for (auto g : n) exp[g] = m_q[ tk::cref_find(lid,g) ];
1585 [ + - ][ + - ]: 69760 : thisProxy[c].comalw( exp );
1586 : : }
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 : : 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 : : const auto npoin = m_u.nunk();
1627 : : const auto ncomp = m_u.nprop();
1628 : : 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 : : 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 : : 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 : : const auto D = m_dsupint[0].data();
1688 : : auto C = m_dsuplim[0].data();
1689 : : std::size_t i = 0;
1690 [ + + ]: 159880 : 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 [ - + ][ - + ]: 1233360 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
1699 : : 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 [ + + ]: 25570 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1715 : 22650 : const auto N = m_dsupedge[1].data() + e*3;
1716 : : const auto D = m_dsupint[1].data();
1717 : : auto C = m_dsuplim[0].data();
1718 : : std::size_t i = 0;
1719 [ + + ]: 90600 : for (const auto& [p,q] : tk::lpoet) {
1720 : 67950 : auto dif = D[(e*3+i)*5+3];
1721 : 67950 : auto coef = C + (e*3+i)*ncomp;
1722 : 67950 : tk::real aec[ncomp];
1723 [ + + ]: 271800 : for (std::size_t c=0; c<ncomp; ++c) {
1724 [ - + ]: 203850 : aec[c] = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1725 : 203850 : auto a = c*2;
1726 : 203850 : auto b = a+1;
1727 [ - + ][ - + ]: 611550 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
1728 : : aec[c] > 0.0 ? m_q(N[q],a) : m_q(N[q],b) );
1729 : : }
1730 : 67950 : tk::real cs = 1.0;
1731 [ - - ][ - + ]: 67950 : for (auto c : fctsys) cs = min( cs, coef[c] );
1732 [ - + ]: 67950 : for (auto c : fctsys) coef[c] = cs;
1733 [ + + ]: 271800 : for (std::size_t c=0; c<ncomp; ++c) {
1734 : 203850 : aec[c] *= coef[c];
1735 : 203850 : m_a(N[p],c) -= aec[c];
1736 : 203850 : m_a(N[q],c) += aec[c];
1737 : : }
1738 : 67950 : ++i;
1739 : 67950 : }
1740 : : }
1741 : :
1742 : : // edges
1743 [ + + ]: 73670 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1744 : 70750 : const auto N = m_dsupedge[2].data() + e*2;
1745 : 70750 : const auto dif = m_dsupint[2][e*5+3];
1746 : 70750 : auto coef = m_dsuplim[2].data() + e*ncomp;
1747 : 70750 : tk::real aec[ncomp];
1748 [ + + ]: 283000 : for (std::size_t c=0; c<ncomp; ++c) {
1749 [ - + ]: 212250 : aec[c] = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1750 : 212250 : auto a = c*2;
1751 : 212250 : auto b = a+1;
1752 [ - + ][ - + ]: 636750 : coef[c] = min( aec[c] < 0.0 ? m_q(N[0],a) : m_q(N[0],b),
1753 : : aec[c] > 0.0 ? m_q(N[1],a) : m_q(N[1],b) );
1754 : : }
1755 : 70750 : tk::real cs = 1.0;
1756 [ - - ][ - + ]: 70750 : for (auto c : fctsys) cs = min( cs, coef[c] );
1757 [ - + ]: 70750 : for (auto c : fctsys) coef[c] = cs;
1758 [ + + ]: 283000 : for (std::size_t c=0; c<ncomp; ++c) {
1759 : 212250 : aec[c] *= coef[c];
1760 : 212250 : m_a(N[0],c) -= aec[c];
1761 : 212250 : m_a(N[1],c) += aec[c];
1762 : : }
1763 : 70750 : }
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 : : decltype(m_ac) exp;
1777 [ + - ][ + + ]: 165880 : for (auto g : n) exp[g] = m_a[ tk::cref_find(lid,g) ];
1778 [ + - ][ + - ]: 69760 : thisProxy[c].comlim( exp );
1779 : : }
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 [ + + ]: 296880 : 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 : 8212 : 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 : 8212 : auto d = Disc();
1813 : :
1814 : 8212 : physics::dirbc( u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
1815 : 8212 : physics::symbc( u, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
1816 : 8212 : physics::noslipbc( u, m_noslipbcnodes, /*pos=*/0 );
1817 : 8212 : }
1818 : :
1819 : : void
1820 : 3650 : ChoCG::dt()
1821 : : // *****************************************************************************
1822 : : // Compute time step size
1823 : : // *****************************************************************************
1824 : : {
1825 : 3650 : auto d = Disc();
1826 : : const auto& vol = d->Vol();
1827 : :
1828 : 3650 : tk::real mindt = std::numeric_limits< tk::real >::max();
1829 [ + + ]: 3650 : auto const_dt = g_cfg.get< tag::dt >();
1830 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1831 : :
1832 : : // use constant dt if configured
1833 [ + + ]: 3650 : if (std::abs(const_dt) > eps) {
1834 : :
1835 : : // cppcheck-suppress redundantInitialization
1836 : 2920 : mindt = const_dt;
1837 : :
1838 : : } else {
1839 : :
1840 : 730 : auto cfl = g_cfg.get< tag::cfl >();
1841 : 730 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1842 : : auto large = std::numeric_limits< tk::real >::max();
1843 : :
1844 [ + + ]: 298900 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1845 : 298170 : auto u = m_u(i,0);
1846 : 298170 : auto v = m_u(i,1);
1847 : 298170 : auto w = m_u(i,2);
1848 : 298170 : auto vel = std::sqrt( u*u + v*v + w*w );
1849 : 298170 : auto L = std::cbrt( vol[i] );
1850 [ + + ][ + + ]: 394990 : auto euler_dt = L / std::max( vel, 1.0e-8 );
1851 : 298170 : mindt = std::min( mindt, euler_dt );
1852 [ + + ][ + + ]: 298170 : auto visc_dt = mu > eps ? L * L / mu : large;
1853 : 298170 : mindt = std::min( mindt, visc_dt );
1854 : : }
1855 : 730 : mindt *= cfl;
1856 : :
1857 : : }
1858 : :
1859 : : // Actiavate SDAG waits for next time step stage
1860 [ + - ]: 3650 : thisProxy[ thisIndex ].wait4rhs();
1861 [ + - ]: 3650 : thisProxy[ thisIndex ].wait4aec();
1862 [ + - ]: 3650 : thisProxy[ thisIndex ].wait4alw();
1863 [ + - ]: 3650 : thisProxy[ thisIndex ].wait4sol();
1864 [ + - ]: 3650 : thisProxy[ thisIndex ].wait4div();
1865 [ + - ]: 3650 : thisProxy[ thisIndex ].wait4sgrad();
1866 [ + - ]: 3650 : thisProxy[ thisIndex ].wait4step();
1867 : :
1868 : : // Contribute to minimum dt across all chares and advance to next step
1869 [ + - ]: 3650 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
1870 : 3650 : CkCallback(CkReductionTarget(ChoCG,advance), thisProxy) );
1871 : 3650 : }
1872 : :
1873 : : void
1874 : 3650 : 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 : : // Detect blowup
1881 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1882 [ - + ]: 3650 : if (newdt < eps) m_finished = 1;
1883 : :
1884 : : // Set new time step size
1885 : 3650 : Disc()->setdt( newdt );
1886 : :
1887 : : // Compute lhs and rhs of transport equations
1888 : 3650 : lhs();
1889 : 3650 : rhs();
1890 : 3650 : }
1891 : :
1892 : : void
1893 : 3650 : ChoCG::lhs()
1894 : : // *****************************************************************************
1895 : : // Fill lhs matrix of transport equations
1896 : : // *****************************************************************************
1897 : : {
1898 : 3650 : auto theta = g_cfg.get< tag::theta >();
1899 [ + + ]: 3650 : if (theta < std::numeric_limits< tk::real >::epsilon()) return;
1900 : :
1901 : 40 : auto d = Disc();
1902 : : const auto& inpoel = d->Inpoel();
1903 : : const auto& coord = d->Coord();
1904 : : const auto& X = coord[0];
1905 : : const auto& Y = coord[1];
1906 : : const auto& Z = coord[2];
1907 : : const auto ncomp = m_u.nprop();
1908 : 40 : const auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1909 : :
1910 : : auto dt = d->Dt();
1911 : 40 : auto& A = Lhs();
1912 : 40 : A.zero();
1913 : :
1914 [ + + ]: 60040 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1915 : 60000 : const auto N = inpoel.data() + e*4;
1916 : : const std::array< tk::real, 3 >
1917 : 60000 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
1918 : 60000 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
1919 : 60000 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
1920 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
1921 : : Assert( J > 0, "Element Jacobian non-positive" );
1922 : : std::array< std::array< tk::real, 3 >, 4 > grad;
1923 : 60000 : grad[1] = tk::cross( ca, da );
1924 : 60000 : grad[2] = tk::cross( da, ba );
1925 : 60000 : grad[3] = tk::cross( ba, ca );
1926 [ + + ]: 240000 : for (std::size_t i=0; i<3; ++i)
1927 : 180000 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
1928 [ + + ]: 300000 : for (std::size_t a=0; a<4; ++a) {
1929 [ + + ]: 1200000 : for (std::size_t b=0; b<4; ++b) {
1930 [ + + ]: 960000 : auto v = J/dt/120.0 * ((a == b) ? 2.0 : 1.0);
1931 : 960000 : v += theta * mu * tk::dot(grad[a],grad[b]) / J / 6.0;
1932 [ + + ]: 3840000 : for (std::size_t c=0; c<ncomp; ++c) A(N[a],N[b],c) -= v;
1933 : : }
1934 : : //for (std::size_t c=0; c<ncomp; ++c) A(N[a],N[a],c) -= J/dt/24.0;
1935 : : }
1936 : : }
1937 : : }
1938 : :
1939 : : void
1940 : 3830 : ChoCG::rhs()
1941 : : // *****************************************************************************
1942 : : // Compute right-hand side of transport equations
1943 : : // *****************************************************************************
1944 : : {
1945 : 3830 : auto d = Disc();
1946 : : const auto& lid = d->Lid();
1947 : :
1948 : : // Compute own portion of right-hand side for all equations
1949 : 3830 : auto dt = m_rkcoef[m_stage] * d->Dt();
1950 : 3830 : chorin::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel,
1951 : 3830 : dt, m_pr, m_u, m_vgrad, m_pgrad, m_rhs );
1952 : :
1953 : : // Communicate rhs to other chares on chare-boundary
1954 [ + + ]: 3830 : if (d->NodeCommMap().empty()) {
1955 : 370 : comrhs_complete();
1956 : : } else {
1957 [ + + ]: 39380 : for (const auto& [c,n] : d->NodeCommMap()) {
1958 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1959 [ + - ][ + + ]: 183740 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
1960 [ + - ][ + - ]: 71840 : thisProxy[c].comrhs( exp );
1961 : : }
1962 : : }
1963 : 3830 : ownrhs_complete();
1964 : 3830 : }
1965 : :
1966 : : void
1967 : 35920 : ChoCG::comrhs(
1968 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1969 : : // *****************************************************************************
1970 : : // Receive contributions to right-hand side vector on chare-boundaries
1971 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1972 : : //! global mesh node IDs, value: contributions for all scalar components.
1973 : : //! \details This function receives contributions to m_rhs, which stores the
1974 : : //! right hand side vector at mesh nodes. While m_rhs stores own
1975 : : //! contributions, m_rhsc collects the neighbor chare contributions during
1976 : : //! communication. This way work on m_rhs and m_rhsc is overlapped.
1977 : : // *****************************************************************************
1978 : : {
1979 : : using tk::operator+=;
1980 [ + + ]: 331560 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
1981 : :
1982 : : // When we have heard from all chares we communicate with, this chare is done
1983 [ + + ]: 35920 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
1984 : 3460 : m_nrhs = 0;
1985 : 3460 : comrhs_complete();
1986 : : }
1987 : 35920 : }
1988 : :
1989 : : void
1990 : 3830 : ChoCG::fct()
1991 : : // *****************************************************************************
1992 : : // Continue with flux-corrected transport if enabled
1993 : : // *****************************************************************************
1994 : : {
1995 : : // Combine own and communicated contributions to rhs
1996 : 3830 : const auto lid = Disc()->Lid();
1997 [ + + ]: 72310 : for (const auto& [g,r] : m_rhsc) {
1998 : 68480 : auto i = tk::cref_find( lid, g );
1999 [ + + ]: 273920 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
2000 : : }
2001 : 3830 : tk::destroy(m_rhsc);
2002 : :
2003 : : auto eps = std::numeric_limits< tk::real >::epsilon();
2004 [ + + ][ + + ]: 3830 : if (g_cfg.get< tag::theta >() < eps and g_cfg.get< tag::fct >())
2005 [ + - ]: 2920 : aec();
2006 : : else
2007 [ + - ]: 910 : solve();
2008 : 3830 : }
2009 : :
2010 : : void
2011 : : // cppcheck-suppress unusedFunction
2012 : 3830 : ChoCG::solve()
2013 : : // *****************************************************************************
2014 : : // Advance systems of equations
2015 : : // *****************************************************************************
2016 : : {
2017 : 3830 : auto d = Disc();
2018 : : const auto npoin = m_u.nunk();
2019 : : const auto ncomp = m_u.nprop();
2020 : : const auto& vol = d->Vol();
2021 : :
2022 : : // Combine own and communicated contributions to limited antidiffusive
2023 : : // contributions
2024 [ + + ]: 56910 : for (const auto& [g,a] : m_ac) {
2025 : 53080 : auto i = tk::cref_find( d->Lid(), g );
2026 [ + + ]: 212320 : for (std::size_t c=0; c<a.size(); ++c) m_a(i,c) += a[c];
2027 : : }
2028 : 3830 : tk::destroy(m_ac);
2029 : :
2030 [ + + ]: 3830 : if (m_stage == 0) m_un = m_u;
2031 : :
2032 [ + + ]: 3830 : if (g_cfg.get< tag::fct >()) {
2033 : :
2034 : : // Apply limited antidiffusive contributions to low-order solution
2035 [ + + ]: 79900 : for (std::size_t i=0; i<npoin; ++i) {
2036 [ + + ]: 307920 : for (std::size_t c=0; c<ncomp; ++c) {
2037 : 230940 : m_a(i,c) = m_rhs(i,c) + m_a(i,c)/vol[i];
2038 : : }
2039 : : }
2040 : : // Continue to advective-diffusive prediction
2041 : 2920 : pred();
2042 : :
2043 : : } else {
2044 : :
2045 : : auto eps = std::numeric_limits<tk::real>::epsilon();
2046 [ + + ][ - + ]: 910 : if (g_cfg.get< tag::theta >() < eps || m_stage+1 < m_rkcoef.size()) {
2047 : :
2048 : : // Apply rhs in explicit solve
2049 : 870 : auto dt = m_rkcoef[m_stage] * d->Dt();
2050 [ + + ]: 481020 : for (std::size_t i=0; i<npoin; ++i) {
2051 [ + + ]: 1920600 : for (std::size_t c=0; c<ncomp; ++c) {
2052 : 1440450 : m_a(i,c) = m_un(i,c) - dt*m_rhs(i,c)/vol[i];
2053 : : }
2054 : : }
2055 : : // Continue to advective-diffusive prediction
2056 : 870 : pred();
2057 : :
2058 : : } else {
2059 : :
2060 : : // Configure Dirichlet BCs
2061 : : std::unordered_map< std::size_t,
2062 : : std::vector< std::pair< int, tk::real > > > dirbc;
2063 : 40 : std::size_t nmask = ncomp + 1;
2064 : : Assert( m_dirbcmask.size() % nmask == 0, "Size mismatch" );
2065 [ + + ]: 24520 : for (std::size_t i=0; i<m_dirbcmask.size()/nmask; ++i) {
2066 [ + - ]: 24480 : auto p = m_dirbcmask[i*nmask+0]; // local node id
2067 : : auto& bc = dirbc[p];
2068 [ + - ]: 24480 : bc.resize( ncomp );
2069 [ + + ]: 97920 : for (std::size_t c=0; c<ncomp; ++c) {
2070 : 73440 : bc[c] = { m_dirbcmask[i*nmask+1+c], 0.0 };
2071 : : }
2072 : : }
2073 [ + - ][ + + ]: 8200 : for (auto p : m_noslipbcnodes) {
2074 : : auto& bc = dirbc[p];
2075 [ + - ]: 8160 : bc.resize( ncomp );
2076 [ + + ]: 32640 : for (std::size_t c=0; c<ncomp; ++c) {
2077 : : bc[c] = { 1, 0.0 };
2078 : : }
2079 : : }
2080 : :
2081 : : // Initialize semi-implicit momentum/transport solve
2082 : : const auto& pc = g_cfg.get< tag::mom_pc >();
2083 [ + - ][ + - ]: 160 : m_cgmom[ thisIndex ].ckLocal()->init( {}, m_rhs.vec(), {}, dirbc, pc,
[ - + ][ - - ]
2084 [ + - ][ + - ]: 120 : CkCallback( CkIndex_ChoCG::msolve(), thisProxy[thisIndex] ) );
[ - + ][ - - ]
2085 : :
2086 : : }
2087 : :
2088 : : }
2089 : 3830 : }
2090 : :
2091 : : void
2092 : 40 : ChoCG::msolve()
2093 : : // *****************************************************************************
2094 : : // Solve for momentum/transport system of equations
2095 : : // *****************************************************************************
2096 : : {
2097 : 40 : auto iter = g_cfg.get< tag::mom_iter >();
2098 : 40 : auto tol = g_cfg.get< tag::mom_tol >();
2099 : 40 : auto verbose = g_cfg.get< tag::mom_verbose >();
2100 : :
2101 [ + - ]: 80 : m_cgmom[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose,
2102 [ + - ][ + - ]: 120 : CkCallback( CkIndex_ChoCG::msolved(), thisProxy[thisIndex] ) );
[ - + ][ - - ]
2103 : 40 : }
2104 : :
2105 : : void
2106 : 40 : ChoCG::msolved()
2107 : : // *****************************************************************************
2108 : : // Continue after momentum/transport solve in semi-implcit solve
2109 : : // *****************************************************************************
2110 : : {
2111 : 40 : auto d = Disc();
2112 : : const auto npoin = m_u.nunk();
2113 : : const auto ncomp = m_u.nprop();
2114 : :
2115 [ + + ][ + - ]: 80 : if (thisIndex == 0) d->mit( m_cgmom[ thisIndex ].ckLocal()->it() );
2116 : :
2117 : : // Update momentum/transport solution in semi-implicit solve
2118 : 40 : auto& du = m_cgmom[ thisIndex ].ckLocal()->solution();
2119 [ + + ]: 24520 : for (std::size_t i=0; i<npoin; ++i) {
2120 [ + + ]: 97920 : for (std::size_t c=0; c<ncomp; ++c) {
2121 : 73440 : m_a(i,c) = m_un(i,c) + du[i*ncomp+c];
2122 : : }
2123 : : }
2124 : :
2125 : : // Continue to advective-diffusive prediction
2126 : 40 : pred();
2127 : 40 : }
2128 : :
2129 : : void
2130 : 3830 : ChoCG::pred()
2131 : : // *****************************************************************************
2132 : : // Compute advective-diffusive prediction of momentum/transport
2133 : : // *****************************************************************************
2134 : : {
2135 : 3830 : auto d = Disc();
2136 : :
2137 : : // Configure and apply scalar source to solution (if defined)
2138 : 3830 : auto src = problems::PHYS_SRC();
2139 [ - + ][ - - ]: 3830 : if (src) src( d->Coord(), d->T(), m_a );
2140 : :
2141 : : // Enforce boundary conditions
2142 [ + - ]: 3830 : BC( m_a, d->T() + m_rkcoef[m_stage] * d->Dt() );
2143 : :
2144 : : // Update momentum/transport solution
2145 : : m_u = m_a;
2146 : : m_a.fill( 0.0 );
2147 : :
2148 : : // Compute velocity gradients if needed
2149 [ + + ]: 3830 : if (g_cfg.get< tag::flux >() == "damp4") {
2150 [ + - ][ + - ]: 3160 : thisProxy[ thisIndex ].wait4vgrad();
[ - - ]
2151 [ + - ]: 3160 : velgrad();
2152 : : } else {
2153 [ + - ]: 670 : corr();
2154 : : }
2155 : 3830 : }
2156 : :
2157 : : void
2158 : 3830 : ChoCG::corr()
2159 : : // *****************************************************************************
2160 : : // Compute pressure correction
2161 : : // *****************************************************************************
2162 : : {
2163 : : // Finalize computing velocity gradients
2164 [ + + ]: 3830 : if (g_cfg.get< tag::flux >() == "damp4") fingrad( m_vgrad, m_vgradc );
2165 : :
2166 [ + + ]: 3830 : if (++m_stage < m_rkcoef.size()) {
2167 : :
2168 : : // Activate SDAG wait for next time step stage
2169 [ + - ]: 180 : thisProxy[ thisIndex ].wait4rhs();
2170 : : // Continue to next time stage of momentum/transport prediction
2171 : 180 : rhs();
2172 : :
2173 : : } else {
2174 : :
2175 : : // Reset Runge-Kutta stage counter
2176 : 3650 : m_stage = 0;
2177 : : // Continue to pressure correction and projection
2178 : 3650 : div( m_u );
2179 : :
2180 : : }
2181 : 3830 : }
2182 : :
2183 : : void
2184 : 3682 : ChoCG::diag()
2185 : : // *****************************************************************************
2186 : : // Compute diagnostics
2187 : : // *****************************************************************************
2188 : : {
2189 : 3682 : auto d = Disc();
2190 : :
2191 : : // Increase number of iterations and physical time
2192 : 3682 : d->next();
2193 : :
2194 : : // Compute diagnostics, e.g., residuals
2195 : 3682 : auto diag_iter = g_cfg.get< tag::diag_iter >();
2196 : 3682 : const auto& dp = m_cgpre[ thisIndex ].ckLocal()->solution();
2197 : 3682 : auto diagnostics = m_diag.precompute( *d, m_u, m_un, m_pr, dp, diag_iter );
2198 : :
2199 : : // Evaluate residuals
2200 [ - + ][ - - ]: 3682 : if (!diagnostics) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
2201 : 3682 : }
2202 : :
2203 : : void
2204 : 3682 : ChoCG::evalres( const std::vector< tk::real >& )
2205 : : // *****************************************************************************
2206 : : // Evaluate residuals
2207 : : // *****************************************************************************
2208 : : {
2209 : 3682 : refine();
2210 : 3682 : }
2211 : :
2212 : : void
2213 : 3682 : ChoCG::refine()
2214 : : // *****************************************************************************
2215 : : // Optionally refine/derefine mesh
2216 : : // *****************************************************************************
2217 : : {
2218 : 3682 : auto d = Disc();
2219 : :
2220 : : // See if this is the last time step
2221 [ + + ]: 3682 : if (d->finished()) m_finished = 1;
2222 : :
2223 : 3682 : auto dtref = g_cfg.get< tag::href_dt >();
2224 : 3682 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
2225 : :
2226 : : // if t>0 refinement enabled and we hit the frequency
2227 [ - + ][ - - ]: 3682 : if (dtref && !(d->It() % dtfreq)) { // refine
2228 : :
2229 : 0 : d->refined() = 1;
2230 : 0 : d->startvol();
2231 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
2232 : :
2233 : : // Activate SDAG waits for re-computing the integrals
2234 [ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
2235 : :
2236 : : } else { // do not refine
2237 : :
2238 : 3682 : d->refined() = 0;
2239 : 3682 : feop_complete();
2240 : 3682 : resize_complete();
2241 : :
2242 : : }
2243 : 3682 : }
2244 : :
2245 : : void
2246 : 0 : ChoCG::resizePostAMR(
2247 : : const std::vector< std::size_t >& /*ginpoel*/,
2248 : : const tk::UnsMesh::Chunk& chunk,
2249 : : const tk::UnsMesh::Coords& coord,
2250 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
2251 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
2252 : : const std::set< std::size_t >& removedNodes,
2253 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
2254 : : nodeCommMap,
2255 : : const std::map< int, std::vector< std::size_t > >& bface,
2256 : : const std::map< int, std::vector< std::size_t > >& bnode,
2257 : : const std::vector< std::size_t >& triinpoel )
2258 : : // *****************************************************************************
2259 : : // Receive new mesh from Refiner
2260 : : //! \param[in] ginpoel Mesh connectivity with global node ids
2261 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
2262 : : //! \param[in] coord New mesh node coordinates
2263 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
2264 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
2265 : : //! \param[in] removedNodes Newly removed mesh node local ids
2266 : : //! \param[in] nodeCommMap New node communication map
2267 : : //! \param[in] bface Boundary-faces mapped to side set ids
2268 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
2269 : : //! \param[in] triinpoel Boundary-face connectivity
2270 : : // *****************************************************************************
2271 : : {
2272 : 0 : auto d = Disc();
2273 : :
2274 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
2275 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
2276 : :
2277 : : // Resize mesh data structures after mesh refinement
2278 : 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
2279 : :
2280 : : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
2281 : : "Incorrect vector length post-AMR: expected length after resizing = " +
2282 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
2283 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
2284 : :
2285 : : // Remove newly removed nodes from solution vectors
2286 : 0 : m_u.rm( removedNodes );
2287 : : //m_pr.rm( removedNodes );
2288 : 0 : m_rhs.rm( removedNodes );
2289 : :
2290 : : // Resize auxiliary solution vectors
2291 : : auto npoin = coord[0].size();
2292 : : m_u.resize( npoin );
2293 : 0 : m_pr.resize( npoin );
2294 : : m_rhs.resize( npoin );
2295 : :
2296 : : // Update solution on new mesh
2297 [ - - ]: 0 : for (const auto& n : addedNodes)
2298 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
2299 : : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
2300 : : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
2301 : : "Indices of parent-edge nodes out of bounds post-AMR");
2302 : 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
2303 : : }
2304 : :
2305 : : // Update physical-boundary node-, face-, and element lists
2306 : : m_bnode = bnode;
2307 : : m_bface = bface;
2308 : 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
2309 : :
2310 : 0 : auto meshid = d->MeshId();
2311 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2312 : 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
2313 : 0 : }
2314 : :
2315 : : void
2316 : 877 : ChoCG::writeFields( CkCallback cb )
2317 : : // *****************************************************************************
2318 : : // Output mesh-based fields to file
2319 : : //! \param[in] cb Function to continue with after the write
2320 : : // *****************************************************************************
2321 : : {
2322 [ + + ]: 877 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
2323 : :
2324 : 293 : auto d = Disc();
2325 : :
2326 : : // Field output
2327 : :
2328 : : std::vector< std::string > nodefieldnames{
2329 [ - + ][ + + ]: 1758 : "velocityx", "velocityy", "velocityz", "divergence", "pressure" };
[ - + ][ - - ]
[ - - ]
2330 : :
2331 : : std::vector< std::vector< tk::real > > nodefields;
2332 : :
2333 [ + - ]: 293 : nodefields.push_back( m_u.extract(0) );
2334 [ + - ]: 293 : nodefields.push_back( m_u.extract(1) );
2335 [ + - ]: 293 : nodefields.push_back( m_u.extract(2) );
2336 : :
2337 : : // Divide weak result by nodal volume
2338 : : const auto& vol = d->Vol();
2339 [ + + ]: 83768 : for (std::size_t i=0; i<m_div.size(); ++i) m_div[i] /= vol[i];
2340 [ + - ]: 293 : nodefields.push_back( m_div );
2341 : :
2342 [ + - ]: 293 : nodefields.push_back( m_pr ) ;
2343 : :
2344 : : //nodefieldnames.push_back( "dp/dx" );
2345 : : //nodefieldnames.push_back( "dp/dy" );
2346 : : //nodefieldnames.push_back( "dp/dz" );
2347 : : //nodefields.push_back( m_pgrad.extract(0) );
2348 : : //nodefields.push_back( m_pgrad.extract(1) );
2349 : : //nodefields.push_back( m_pgrad.extract(2) );
2350 : :
2351 : : //nodefieldnames.push_back( "fx" );
2352 : : //nodefieldnames.push_back( "fy" );
2353 : : //nodefieldnames.push_back( "fz" );
2354 : : //nodefields.push_back( m_flux.extract(0) );
2355 : : //nodefields.push_back( m_flux.extract(1) );
2356 : : //nodefields.push_back( m_flux.extract(2) );
2357 : :
2358 : : //nodefieldnames.push_back( "du/dx" );
2359 : : //nodefieldnames.push_back( "du/dy" );
2360 : : //nodefieldnames.push_back( "du/dz" );
2361 : : //nodefieldnames.push_back( "dv/dx" );
2362 : : //nodefieldnames.push_back( "dv/dy" );
2363 : : //nodefieldnames.push_back( "dv/dz" );
2364 : : //nodefieldnames.push_back( "dw/dx" );
2365 : : //nodefieldnames.push_back( "dw/dy" );
2366 : : //nodefieldnames.push_back( "dw/dz" );
2367 : : //nodefields.push_back( m_vgrad.extract(0) );
2368 : : //nodefields.push_back( m_vgrad.extract(1) );
2369 : : //nodefields.push_back( m_vgrad.extract(2) );
2370 : : //nodefields.push_back( m_vgrad.extract(3) );
2371 : : //nodefields.push_back( m_vgrad.extract(4) );
2372 : : //nodefields.push_back( m_vgrad.extract(5) );
2373 : : //nodefields.push_back( m_vgrad.extract(6) );
2374 : : //nodefields.push_back( m_vgrad.extract(7) );
2375 : : //nodefields.push_back( m_vgrad.extract(8) );
2376 : :
2377 : : // also output analytic pressure (if defined)
2378 [ + - ]: 293 : auto pressure_sol = problems::PRESSURE_SOL();
2379 [ + + ]: 293 : if (pressure_sol) {
2380 : : const auto& coord = d->Coord();
2381 : : const auto& x = coord[0];
2382 : : const auto& y = coord[1];
2383 : : const auto& z = coord[2];
2384 [ + - ]: 64 : auto ap = m_pr;
2385 [ + + ]: 4326 : for (std::size_t i=0; i<ap.size(); ++i) {
2386 [ - + ]: 8524 : ap[i] = pressure_sol( x[i], y[i], z[i] );
2387 : : }
2388 [ + - ][ - - ]: 64 : nodefieldnames.push_back( "analytic" );
2389 [ + - ]: 64 : nodefields.push_back( ap );
2390 : : }
2391 : :
2392 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
2393 : :
2394 : : // Surface output
2395 : :
2396 : : std::vector< std::string > nodesurfnames;
2397 : : std::vector< std::vector< tk::real > > nodesurfs;
2398 : :
2399 : : const auto& f = g_cfg.get< tag::fieldout >();
2400 : :
2401 [ + + ]: 293 : if (!f.empty()) {
2402 : : std::size_t ncomp = 5;
2403 [ + - ]: 3 : nodesurfnames.push_back( "velocityx" );
2404 [ + - ]: 3 : nodesurfnames.push_back( "velocityy" );
2405 [ + - ]: 3 : nodesurfnames.push_back( "velocityz" );
2406 [ + - ]: 3 : nodesurfnames.push_back( "divergence" );
2407 [ + - ]: 3 : nodesurfnames.push_back( "pressure" );
2408 : :
2409 [ + - ]: 3 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
2410 [ + - ]: 3 : std::set< int > outsets( begin(f), end(f) );
2411 [ + + ]: 6 : for (auto sideset : outsets) {
2412 : : auto b = bnode.find(sideset);
2413 [ - + ]: 3 : if (b == end(bnode)) continue;
2414 : : const auto& nodes = b->second;
2415 : : auto i = nodesurfs.size();
2416 : : nodesurfs.insert( end(nodesurfs), ncomp,
2417 [ + - ]: 6 : std::vector< tk::real >( nodes.size() ) );
2418 : : std::size_t j = 0;
2419 [ + + ]: 1278 : for (auto n : nodes) {
2420 [ + - ]: 1275 : const auto s = m_u[n];
2421 : : std::size_t p = 0;
2422 : 1275 : nodesurfs[i+(p++)][j] = s[0];
2423 : 1275 : nodesurfs[i+(p++)][j] = s[1];
2424 : 1275 : nodesurfs[i+(p++)][j] = s[2];
2425 : 1275 : nodesurfs[i+(p++)][j] = m_div[n];
2426 : 1275 : nodesurfs[i+(p++)][j] = m_pr[n];
2427 : 1275 : ++j;
2428 : : }
2429 : : }
2430 : : }
2431 : :
2432 : : // Send mesh and fields data (solution dump) for output to file
2433 [ + - ][ + - ]: 586 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
2434 [ + - ]: 293 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
2435 : : {}, nodefields, {}, nodesurfs, cb );
2436 [ + - ][ + - ]: 1172 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
2437 : :
2438 : : void
2439 : 3682 : ChoCG::out()
2440 : : // *****************************************************************************
2441 : : // Output mesh field data
2442 : : // *****************************************************************************
2443 : : {
2444 : 3682 : auto d = Disc();
2445 : :
2446 : : // Time history
2447 [ + - ][ + - ]: 3682 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
2448 : : auto ncomp = m_u.nprop();
2449 : : const auto& inpoel = d->Inpoel();
2450 : 0 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
2451 : : std::size_t j = 0;
2452 [ - - ]: 0 : for (const auto& p : d->Hist()) {
2453 [ - - ]: 0 : auto e = p.get< tag::elem >(); // host element id
2454 : : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
2455 [ - - ]: 0 : hist[j].resize( ncomp+1, 0.0 );
2456 [ - - ]: 0 : for (std::size_t i=0; i<4; ++i) {
2457 [ - - ]: 0 : const auto u = m_u[ inpoel[e*4+i] ];
2458 : 0 : hist[j][0] += n[i] * u[0];
2459 : 0 : hist[j][1] += n[i] * u[1]/u[0];
2460 : 0 : hist[j][2] += n[i] * u[2]/u[0];
2461 : 0 : hist[j][3] += n[i] * u[3]/u[0];
2462 : 0 : hist[j][4] += n[i] * u[4]/u[0];
2463 : 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];
2464 : 0 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
2465 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
2466 : : }
2467 : 0 : ++j;
2468 : : }
2469 [ - - ]: 0 : d->history( std::move(hist) );
2470 : 0 : }
2471 : :
2472 : : // Field data
2473 [ + + ][ + - ]: 3682 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
2474 [ + - ][ + - ]: 1533 : writeFields( CkCallback(CkIndex_ChoCG::integrals(), thisProxy[thisIndex]) );
2475 : : } else {
2476 : 3171 : integrals();
2477 : : }
2478 : 3682 : }
2479 : :
2480 : : void
2481 : 3682 : ChoCG::integrals()
2482 : : // *****************************************************************************
2483 : : // Compute integral quantities for output
2484 : : // *****************************************************************************
2485 : : {
2486 : 3682 : auto d = Disc();
2487 : :
2488 [ + + ][ + - ]: 3682 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ]
2489 : :
2490 : : using namespace integrals;
2491 [ + - ]: 20 : std::vector< std::map< int, tk::real > > ints( NUMINT );
2492 : : // Prepend integral vector with metadata on the current time step:
2493 : : // current iteration count, current physical time, time step size
2494 [ + - ][ + - ]: 20 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
2495 [ + - ][ + - ]: 20 : ints[ TIME ][ 0 ] = d->T();
2496 [ + - ][ + - ]: 20 : ints[ DT ][ 0 ] = d->Dt();
2497 : :
2498 : : // Compute integrals requested for surfaces requested
2499 : : const auto& reqv = g_cfg.get< tag::integout_integrals >();
2500 : 20 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
2501 [ + - ][ - + ]: 40 : if (req.count("mass_flow_rate")) {
2502 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
2503 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
2504 : : const auto& nodes = sint.first;
2505 : : const auto& ndA = sint.second;
2506 : : auto n = ndA.data();
2507 [ - - ]: 0 : for (auto p : nodes) {
2508 : 0 : mfr += n[0]*m_u(p,0) + n[1]*m_u(p,1) + n[2]*m_u(p,2);
2509 : 0 : n += 3;
2510 : : }
2511 : : }
2512 : : }
2513 [ + - ][ + - ]: 40 : if (req.count("force")) {
2514 : 20 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
2515 [ + + ]: 40 : for (const auto& [s,sint] : m_surfint) {
2516 [ + - ]: 20 : auto& fx = ints[ FORCE_X ][ s ];
2517 [ + - ]: 20 : auto& fy = ints[ FORCE_Y ][ s ];
2518 [ + - ]: 20 : auto& fz = ints[ FORCE_Z ][ s ];
2519 : : const auto& nodes = sint.first;
2520 : : const auto& ndA = sint.second;
2521 : : auto n = ndA.data();
2522 [ + + ]: 8520 : for (auto p : nodes) {
2523 : : // pressure force
2524 : 8500 : fx -= n[0]*m_pr[p];
2525 : 8500 : fy -= n[1]*m_pr[p];
2526 : 8500 : fz -= n[2]*m_pr[p];
2527 : : // viscous force
2528 : 8500 : fx += mu*(m_vgrad(p,0)*n[0] + m_vgrad(p,1)*n[1] + m_vgrad(p,2)*n[2]);
2529 : 8500 : fy += mu*(m_vgrad(p,3)*n[0] + m_vgrad(p,4)*n[1] + m_vgrad(p,5)*n[2]);
2530 : 8500 : fz += mu*(m_vgrad(p,6)*n[0] + m_vgrad(p,7)*n[1] + m_vgrad(p,8)*n[2]);
2531 : 8500 : n += 3;
2532 : : }
2533 : : }
2534 : : }
2535 : :
2536 [ + - ]: 20 : auto stream = serialize( d->MeshId(), ints );
2537 [ + - ][ + - ]: 40 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
2538 [ + - ][ - - ]: 20 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
2539 : :
2540 : 20 : } else {
2541 : :
2542 : 3662 : step();
2543 : :
2544 : : }
2545 : 3682 : }
2546 : :
2547 : : void
2548 : 3316 : ChoCG::evalLB( int nrestart )
2549 : : // *****************************************************************************
2550 : : // Evaluate whether to do load balancing
2551 : : //! \param[in] nrestart Number of times restarted
2552 : : // *****************************************************************************
2553 : : {
2554 : 3316 : auto d = Disc();
2555 : :
2556 : : // Detect if just returned from a checkpoint and if so, zero timers and
2557 : : // finished flag
2558 [ + + ]: 3316 : if (d->restarted( nrestart )) m_finished = 0;
2559 : :
2560 : 3316 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
2561 [ + + ]: 3316 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
2562 : :
2563 : : // Load balancing if user frequency is reached or after the second time-step
2564 [ + + ][ + + ]: 3316 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
2565 : :
2566 : 3176 : AtSync();
2567 [ - + ]: 3176 : if (nonblocking) dt();
2568 : :
2569 : : } else {
2570 : :
2571 : 140 : dt();
2572 : :
2573 : : }
2574 : 3316 : }
2575 : :
2576 : : void
2577 : 3311 : ChoCG::evalRestart()
2578 : : // *****************************************************************************
2579 : : // Evaluate whether to save checkpoint/restart
2580 : : // *****************************************************************************
2581 : : {
2582 : 3311 : auto d = Disc();
2583 : :
2584 : 3311 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
2585 : 3311 : const auto benchmark = g_cfg.get< tag::benchmark >();
2586 : :
2587 [ + + ][ - + ]: 3311 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
2588 : :
2589 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
2590 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
2591 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
2592 : :
2593 : : } else {
2594 : :
2595 : 3311 : evalLB( /* nrestart = */ -1 );
2596 : :
2597 : : }
2598 : 3311 : }
2599 : :
2600 : : void
2601 : 3682 : ChoCG::step()
2602 : : // *****************************************************************************
2603 : : // Evaluate whether to continue with next time step
2604 : : // *****************************************************************************
2605 : : {
2606 : 3682 : auto d = Disc();
2607 : :
2608 : : // Output one-liner status report to screen
2609 [ + + ]: 3682 : if(thisIndex == 0) d->status();
2610 : :
2611 [ + + ]: 3682 : if (not m_finished) {
2612 : :
2613 : 3311 : evalRestart();
2614 : :
2615 : : } else {
2616 : :
2617 : 371 : auto meshid = d->MeshId();
2618 [ + - ]: 742 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2619 : 371 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
2620 : :
2621 : : }
2622 : 3682 : }
2623 : :
2624 : : #include "NoWarning/chocg.def.h"
|
