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