Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ChoCG.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2025 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief ChoCG: Projection-based solver for incompressible flow
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "XystBuildConfig.hpp"
14 : : #include "ChoCG.hpp"
15 : : #include "Vector.hpp"
16 : : #include "Reader.hpp"
17 : : #include "ContainerUtil.hpp"
18 : : #include "UnsMesh.hpp"
19 : : #include "ExodusIIMeshWriter.hpp"
20 : : #include "InciterConfig.hpp"
21 : : #include "DerivedData.hpp"
22 : : #include "Discretization.hpp"
23 : : #include "DiagReducer.hpp"
24 : : #include "IntegralReducer.hpp"
25 : : #include "Integrals.hpp"
26 : : #include "Refiner.hpp"
27 : : #include "Reorder.hpp"
28 : : #include "Around.hpp"
29 : : #include "Chorin.hpp"
30 : : #include "Problems.hpp"
31 : : #include "EOS.hpp"
32 : : #include "BC.hpp"
33 : : #include "Print.hpp"
34 : :
35 : : namespace inciter {
36 : :
37 : : extern ctr::Config g_cfg;
38 : :
39 : : static CkReduction::reducerType IntegralsMerger;
40 : :
41 : : //! Runge-Kutta coefficients
42 : : //! Runge-Kutta coefficients
43 : : static const std::array< std::vector< tk::real >, 4 > rkcoef{{
44 : : { 1.0 },
45 : : { 1.0/2.0, 1.0 },
46 : : { 1.0/3.0, 1.0/2.0, 1.0 },
47 : : { 1.0/4.0, 1.0/3.0, 1.0/2.0, 1.0 }
48 : : }};
49 : :
50 : : } // inciter::
51 : :
52 : : using inciter::g_cfg;
53 : : using inciter::ChoCG;
54 : :
55 : 416 : 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 : 416 : const std::vector< std::size_t >& triinpoel ) :
61 : : m_disc( disc ),
62 [ + - ]: 416 : m_cgpre( cgpre ),
63 [ + - ]: 416 : m_cgmom( cgmom ),
64 : 416 : m_nrhs( 0 ),
65 : 416 : m_nnorm( 0 ),
66 : 416 : m_nsgrad( 0 ),
67 : 416 : m_npgrad( 0 ),
68 : 416 : m_nvgrad( 0 ),
69 : 416 : m_nflux( 0 ),
70 : 416 : m_ndiv( 0 ),
71 : 416 : m_nbpint( 0 ),
72 [ + - ]: 416 : m_np( 0 ),
73 : : m_bnode( bnode ),
74 : : m_bface( bface ),
75 [ + - ][ + - ]: 416 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
76 [ + - ]: 416 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
77 : : m_un( m_u.nunk(), m_u.nprop() ),
78 [ + - ][ + - ]: 416 : m_pr( m_u.nunk(), 0.0 ),
[ - - ]
79 : : m_rhs( m_u.nunk(), m_u.nprop() ),
80 : : m_sgrad( m_u.nunk(), 3UL ),
81 : : m_pgrad( m_u.nunk(), 3UL ),
82 [ + - ]: 416 : m_vgrad( m_u.nunk(), m_u.nprop()*3 ),
83 : : m_flux( m_u.nunk(), 3UL ),
84 [ + - ]: 416 : m_div( m_u.nunk() ),
85 : 416 : m_stage( 0 ),
86 : 416 : m_finished( 0 ),
87 : 416 : m_freezeflow( 1.0 ),
88 [ + - ][ + - ]: 832 : m_rkcoef( rkcoef[ g_cfg.get< tag::rk >() - 1 ] )
89 : : // *****************************************************************************
90 : : // Constructor
91 : : //! \param[in] disc Discretization proxy
92 : : //! \param[in] cgpre ConjugateGradients Charm++ proxy for pressure solve
93 : : //! \param[in] cgmom ConjugateGradients Charm++ proxy for momentum solve
94 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
95 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
96 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
97 : : // *****************************************************************************
98 : : {
99 : 416 : usesAtSync = true; // enable migration at AtSync
100 : :
101 [ + - ]: 416 : auto d = Disc();
102 : :
103 : : // Create new local ids based on mesh locality
104 : : std::unordered_map< std::size_t, std::size_t > map;
105 : : std::size_t n = 0;
106 : :
107 [ + - ][ + - ]: 416 : auto psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
108 [ + + ]: 32426 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
109 [ + - ]: 1328 : if (!map.count(p)) map[p] = n++;
110 [ + + ][ + + ]: 311634 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
111 [ + - ]: 30682 : if (!map.count(q)) map[q] = n++;
112 : : }
113 : : }
114 : :
115 : : Assert( map.size() == d->Gid().size(),
116 : : "Mesh-locality reorder map size mismatch" );
117 : :
118 : : // Remap data in bound Discretization object
119 [ + - ]: 416 : d->remap( map );
120 : : // Remap boundary triangle face connectivity
121 [ + - ]: 416 : tk::remap( m_triinpoel, map );
122 : : // Recompute points surrounding points
123 [ + - ][ + - ]: 416 : psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
124 : :
125 : : // Compute total box IC volume
126 [ + - ]: 416 : d->boxvol();
127 : :
128 : : // Setup LHS matrix for pressure solve
129 [ + - ][ + - ]: 832 : m_cgpre[ thisIndex ].insert( prelhs( psup ),
[ + - ][ - + ]
[ - - ]
130 : : d->Gid(),
131 : : d->Lid(),
132 : : d->NodeCommMap() );
133 : :
134 : : // Setup empty LHS matrix for momentum solve if needed
135 [ + + ]: 416 : if (g_cfg.get< tag::theta >() > std::numeric_limits< tk::real >::epsilon()) {
136 [ + - ][ + - ]: 4 : m_cgmom[ thisIndex ].insert( momlhs( psup ),
[ + - ][ - + ]
[ - - ]
137 : : d->Gid(),
138 : : d->Lid(),
139 : : d->NodeCommMap() );
140 : : }
141 : :
142 : : // Activate SDAG waits for setup
143 [ + - ][ + - ]: 416 : thisProxy[ thisIndex ].wait4int();
144 : 832 : }
145 : :
146 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
147 : 416 : ChoCG::prelhs( const std::pair< std::vector< std::size_t >,
148 : : std::vector< std::size_t > >& psup )
149 : : // *****************************************************************************
150 : : // Setup lhs matrix for pressure solve
151 : : //! \param[in] psup Points surrounding points
152 : : //! \return { A, x, b } in linear system A * x = b to solve for pressure
153 : : // *****************************************************************************
154 : : {
155 : 416 : auto d = Disc();
156 : : const auto& inpoel = d->Inpoel();
157 : : const auto& coord = d->Coord();
158 : : const auto& X = coord[0];
159 : : const auto& Y = coord[1];
160 : : const auto& Z = coord[2];
161 : :
162 : : // Matrix with compressed sparse row storage
163 : 416 : tk::CSR A( /*DOF=*/ 1, psup );
164 : :
165 : : // Fill matrix with Laplacian
166 [ + + ]: 78797 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
167 : 78381 : const auto N = inpoel.data() + e*4;
168 : : const std::array< tk::real, 3 >
169 : 78381 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
170 : 78381 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
171 : 78381 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
172 : 78381 : const auto J = tk::triple( ba, ca, da ) * 6.0;
173 : : std::array< std::array< tk::real, 3 >, 4 > grad;
174 : 78381 : grad[1] = tk::cross( ca, da );
175 : 78381 : grad[2] = tk::cross( da, ba );
176 : 78381 : grad[3] = tk::cross( ba, ca );
177 [ + + ]: 313524 : for (std::size_t i=0; i<3; ++i)
178 : 235143 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
179 [ + + ]: 391905 : for (std::size_t a=0; a<4; ++a)
180 [ + + ]: 1567620 : for (std::size_t b=0; b<4; ++b)
181 [ + - ]: 1254096 : A(N[a],N[b]) -= tk::dot( grad[a], grad[b] ) / J;
182 : : }
183 : :
184 : : auto nunk = X.size();
185 [ + - ][ + - ]: 416 : std::vector< tk::real > x( nunk, 0.0 ), b( nunk, 0.0 );
[ - - ]
186 : :
187 : 416 : return { std::move(A), std::move(x), std::move(b) };
188 : 416 : }
189 : :
190 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
191 : 2 : ChoCG::momlhs( const std::pair< std::vector< std::size_t >,
192 : : std::vector< std::size_t > >& psup )
193 : : // *****************************************************************************
194 : : // Setup empty lhs matrix for momentum solve
195 : : //! \param[in] psup Points surrounding points
196 : : //! \return { A, x, b } in linear system A * x = b to solve for momentum
197 : : // *****************************************************************************
198 : : {
199 : : auto ncomp = m_u.nprop();
200 : :
201 : : // Matrix with compressed sparse row storage
202 : 2 : tk::CSR A( /*DOF=*/ ncomp, psup );
203 : :
204 : 2 : auto nunk = (psup.second.size() - 1) * ncomp;
205 [ + - ][ + - ]: 2 : std::vector< tk::real > x( nunk, 0.0 ), b( nunk, 0.0 );
[ - - ]
206 : :
207 : 2 : return { std::move(A), std::move(x), std::move(b) };
208 : 2 : }
209 : :
210 : : void
211 : 832 : ChoCG::setupDirBC( const std::vector< std::vector< int > >& cfgmask,
212 : : const std::vector< std::vector< double > >& cfgval,
213 : : std::size_t ncomp,
214 : : std::vector< std::size_t >& mask,
215 : : std::vector< double >& val )
216 : : // *****************************************************************************
217 : : // Prepare Dirichlet boundary condition data structures
218 : : //! \param[in] cfgmask Boundary condition mask config data to use
219 : : //! \param[in] cfgval Boundary condition values config data to use
220 : : //! \param[in] ncomp Number of scalar component BCs expected per mesh node
221 : : //! \param[in,out] mask Mesh nodes and their Dirichlet BC masks
222 : : //! \param[in,out] val Mesh nodes and their Dirichlet BC values
223 : : // *****************************************************************************
224 : : {
225 : : // Query Dirichlet BC nodes associated to side sets
226 : : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
227 [ + + ]: 1255 : for (const auto& s : cfgmask) {
228 : : auto k = m_bface.find(s[0]);
229 [ + + ]: 423 : if (k != end(m_bface)) {
230 [ + - ]: 259 : auto& n = dir[ k->first ];
231 [ + + ]: 34035 : for (auto f : k->second) {
232 [ + - ]: 33776 : n.insert( m_triinpoel[f*3+0] );
233 [ + - ]: 33776 : n.insert( m_triinpoel[f*3+1] );
234 [ + - ]: 33776 : n.insert( m_triinpoel[f*3+2] );
235 : : }
236 : : }
237 : : }
238 : :
239 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
240 [ + - ]: 832 : const auto& lid = Disc()->Lid();
241 [ + + ]: 1255 : for (const auto& s : cfgmask) {
242 : : auto k = m_bnode.find(s[0]);
243 [ + + ]: 423 : if (k != end(m_bnode)) {
244 [ + - ]: 260 : auto& n = dir[ k->first ];
245 [ + - ][ + + ]: 31491 : for (auto g : k->second) {
246 : : n.insert( tk::cref_find(lid,g) );
247 : : }
248 : : }
249 : : }
250 : :
251 : : // Associate sidesets to Dirichlet BC values if configured by user
252 : : std::unordered_map< int, std::vector< double > > dirval;
253 [ + + ]: 896 : for (const auto& s : cfgval) {
254 [ + + ]: 64 : auto k = dir.find( static_cast<int>(s[0]) );
255 [ + + ]: 64 : if (k != end(dir)) {
256 [ + - ]: 22 : auto& v = dirval[ k->first ];
257 [ + - ]: 22 : v.resize( s.size()-1 );
258 [ + + ]: 52 : for (std::size_t i=1; i<s.size(); ++i) v[i-1] = s[i];
259 : : }
260 : : }
261 : :
262 : : // Collect unique set of nodes + Dirichlet BC components mask and value
263 : 832 : auto nmask = ncomp + 1;
264 : : std::unordered_map< std::size_t,
265 : : std::pair< std::vector< int >,
266 : : std::vector< double > > > dirbcset;
267 [ + + ]: 1255 : for (const auto& vec : cfgmask) {
268 [ - + ][ - - ]: 423 : ErrChk( vec.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
269 : : auto n = dir.find( vec[0] );
270 [ + + ]: 423 : if (n != end(dir)) {
271 [ + - ][ + + ]: 260 : std::vector< double > v( ncomp, 0.0 );
272 : : auto m = dirval.find( vec[0] );
273 [ + + ]: 260 : if (m != end(dirval)) {
274 [ - + ][ - - ]: 22 : ErrChk( m->second.size() == ncomp, "Incorrect Dirichlet BC val ncomp" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
275 [ + - ]: 22 : v = m->second;
276 : : }
277 [ + - ][ + + ]: 22213 : for (auto p : n->second) {
278 : : auto& mv = dirbcset[p]; // mask & value
279 [ + - ]: 21953 : mv.second = v;
280 [ + + ]: 21953 : auto& mval = mv.first;
281 [ + + ][ + - ]: 21953 : if (mval.empty()) mval.resize( ncomp, 0 );
[ - - ]
282 [ + + ]: 90721 : for (std::size_t c=0; c<ncomp; ++c)
283 [ + + ]: 68768 : if (!mval[c]) mval[c] = vec[c+1]; // overwrite mask if 0 -> 1
284 : : }
285 : : }
286 : : }
287 : :
288 : : // Compile streamable list of nodes + Dirichlet BC components mask and values
289 : : tk::destroy( mask );
290 [ + + ]: 21353 : for (const auto& [p,mv] : dirbcset) {
291 [ + - ]: 20521 : mask.push_back( p );
292 [ + - ]: 20521 : mask.insert( end(mask), begin(mv.first), end(mv.first) );
293 [ + - ][ + - ]: 20521 : val.push_back( static_cast< double >( p ) );
294 [ + - ]: 20521 : val.insert( end(val), begin(mv.second), end(mv.second) );
295 : : }
296 : :
297 [ - + ][ - - ]: 832 : ErrChk( mask.size() % nmask == 0, "Dirichlet BC mask incomplete" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
298 [ - + ][ - - ]: 832 : ErrChk( val.size() % nmask == 0, "Dirichlet BC val incomplete" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
299 [ - + ][ - - ]: 832 : ErrChk( mask.size() == val.size(), "Dirichlet BC mask & val size mismatch" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
300 : 832 : }
301 : :
302 : : void
303 : 416 : ChoCG::feop()
304 : : // *****************************************************************************
305 : : // Start (re-)computing finite element domain and boundary operators
306 : : // *****************************************************************************
307 : : {
308 : 416 : auto d = Disc();
309 : :
310 : : // Prepare Dirichlet boundary conditions data structures
311 : 416 : setupDirBC( g_cfg.get< tag::bc_dir >(), g_cfg.get< tag::bc_dirval >(),
312 : 416 : m_u.nprop(), m_dirbcmask, m_dirbcval );
313 : 416 : setupDirBC( g_cfg.get< tag::pre_bc_dir >(), g_cfg.get< tag::pre_bc_dirval >(),
314 : 416 : 1, m_dirbcmaskp, m_dirbcvalp );
315 : :
316 : : // Compute local contributions to boundary normals and integrals
317 : 416 : bndint();
318 : : // Compute local contributions to domain edge integrals
319 : 416 : domint();
320 : :
321 : : // Send boundary point normal contributions to neighbor chares
322 [ + + ]: 416 : if (d->NodeCommMap().empty()) {
323 : 15 : comnorm_complete();
324 : : } else {
325 [ + + ]: 4385 : for (const auto& [c,nodes] : d->NodeCommMap()) {
326 : : decltype(m_bnorm) exp;
327 [ + + ]: 21740 : for (auto i : nodes) {
328 [ + + ]: 53996 : for (const auto& [s,b] : m_bnorm) {
329 : : auto k = b.find(i);
330 [ + + ]: 46467 : if (k != end(b)) exp[s][i] = k->second;
331 : : }
332 : : }
333 [ + - ][ + - ]: 7968 : thisProxy[c].comnorm( exp );
334 : : }
335 : : }
336 : 416 : ownnorm_complete();
337 : 416 : }
338 : :
339 : : void
340 : 416 : ChoCG::bndint()
341 : : // *****************************************************************************
342 : : // Compute local contributions to boundary normals and integrals
343 : : // *****************************************************************************
344 : : {
345 : 416 : auto d = Disc();
346 : : const auto& coord = d->Coord();
347 : : const auto& gid = d->Gid();
348 : : const auto& x = coord[0];
349 : : const auto& y = coord[1];
350 : : const auto& z = coord[2];
351 : :
352 : : // Lambda to compute the inverse distance squared between boundary face
353 : : // centroid and boundary point. Here p is the global node id and c is the
354 : : // the boundary face centroid.
355 : 146760 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
356 : 146760 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
357 : 146760 : (c[1] - y[p]) * (c[1] - y[p]) +
358 : 146760 : (c[2] - z[p]) * (c[2] - z[p]) );
359 : 416 : };
360 : :
361 : 416 : tk::destroy( m_bnorm );
362 : 416 : tk::destroy( m_bndpoinint );
363 : :
364 [ + + ]: 1371 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
365 [ + + ]: 49875 : for (auto f : faceids) { // for all side set triangles
366 : 48920 : const auto N = m_triinpoel.data() + f*3;
367 : : const std::array< tk::real, 3 >
368 : 48920 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
369 : 48920 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
370 : : auto n = tk::cross( ba, ca );
371 : : auto A2 = tk::length( n );
372 : 48920 : n[0] /= A2;
373 : 48920 : n[1] /= A2;
374 : 48920 : n[2] /= A2;
375 : : const tk::real centroid[3] = {
376 : 48920 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
377 : 48920 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
378 : 48920 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
379 [ + + ]: 195680 : for (const auto& [i,j] : tk::lpoet) {
380 : 146760 : auto p = N[i];
381 : 146760 : tk::real r = invdistsq( centroid, p );
382 : : auto& v = m_bnorm[setid]; // associate side set id
383 : : auto& bpn = v[gid[p]]; // associate global node id of bnd pnt
384 : 146760 : bpn[0] += r * n[0]; // inv.dist.sq-weighted normal
385 : 146760 : bpn[1] += r * n[1];
386 : 146760 : bpn[2] += r * n[2];
387 : 146760 : bpn[3] += r; // inv.dist.sq of node from centroid
388 : : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
389 : 146760 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
390 : 146760 : b[1] += n[1] * A2 / 6.0;
391 : 146760 : b[2] += n[2] * A2 / 6.0;
392 : : }
393 : : }
394 : : }
395 : 416 : }
396 : :
397 : : void
398 : 416 : ChoCG::domint()
399 : : // *****************************************************************************
400 : : //! Compute local contributions to domain edge integrals
401 : : // *****************************************************************************
402 : : {
403 : 416 : auto d = Disc();
404 : :
405 : : const auto& gid = d->Gid();
406 : : const auto& inpoel = d->Inpoel();
407 : :
408 : : const auto& coord = d->Coord();
409 : : const auto& x = coord[0];
410 : : const auto& y = coord[1];
411 : : const auto& z = coord[2];
412 : :
413 : 416 : tk::destroy( m_domedgeint );
414 : :
415 [ + + ]: 78797 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
416 : 78381 : const auto N = inpoel.data() + e*4;
417 : : const std::array< tk::real, 3 >
418 : 78381 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
419 : 78381 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
420 : 78381 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
421 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
422 : : Assert( J > 0, "Element Jacobian non-positive" );
423 : : std::array< std::array< tk::real, 3 >, 4 > grad;
424 : 78381 : grad[1] = tk::cross( ca, da );
425 : 78381 : grad[2] = tk::cross( da, ba );
426 : 78381 : grad[3] = tk::cross( ba, ca );
427 [ + + ]: 313524 : for (std::size_t i=0; i<3; ++i)
428 : 235143 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
429 [ + + ]: 548667 : for (const auto& [p,q] : tk::lpoed) {
430 [ + + ]: 470286 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
431 : : tk::real sig = 1.0;
432 [ + + ]: 470286 : if (ed[0] > ed[1]) {
433 : : std::swap( ed[0], ed[1] );
434 : : sig = -1.0;
435 : : }
436 : : auto& n = m_domedgeint[ ed ];
437 : 470286 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
438 : 470286 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
439 : 470286 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
440 : 470286 : n[3] += J / 120.0;
441 : 470286 : n[4] += tk::dot( grad[p], grad[q] ) / J / 6.0;
442 : : }
443 : : }
444 : 416 : }
445 : :
446 : : void
447 : 3984 : ChoCG::comnorm( const decltype(m_bnorm)& inbnd )
448 : : // *****************************************************************************
449 : : // Receive contributions to boundary point normals on chare-boundaries
450 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
451 : : // *****************************************************************************
452 : : {
453 : : // Buffer up incoming boundary point normal vector contributions
454 [ + + ]: 7168 : for (const auto& [s,b] : inbnd) {
455 : : auto& bndnorm = m_bnormc[s];
456 [ + + ]: 13411 : for (const auto& [p,n] : b) {
457 : : auto& norm = bndnorm[p];
458 : 10227 : norm[0] += n[0];
459 : 10227 : norm[1] += n[1];
460 : 10227 : norm[2] += n[2];
461 : 10227 : norm[3] += n[3];
462 : : }
463 : : }
464 : :
465 [ + + ]: 3984 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
466 : 401 : m_nnorm = 0;
467 : 401 : comnorm_complete();
468 : : }
469 : 3984 : }
470 : :
471 : : void
472 : 803 : ChoCG::registerReducers()
473 : : // *****************************************************************************
474 : : // Configure Charm++ reduction types initiated from this chare array
475 : : //! \details Since this is a [initnode] routine, the runtime system executes the
476 : : //! routine exactly once on every logical node early on in the Charm++ init
477 : : //! sequence. Must be static as it is called without an object. See also:
478 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
479 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
480 : : // *****************************************************************************
481 : : {
482 : 803 : NodeDiagnostics::registerReducers();
483 : 803 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
484 : 803 : }
485 : :
486 : : void
487 : : // cppcheck-suppress unusedFunction
488 : 3446 : ChoCG::ResumeFromSync()
489 : : // *****************************************************************************
490 : : // Return from migration
491 : : //! \details This is called when load balancing (LB) completes. The presence of
492 : : //! this function does not affect whether or not we block on LB.
493 : : // *****************************************************************************
494 : : {
495 [ - + ][ - - ]: 3446 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
496 : :
497 [ + - ]: 3446 : if (!g_cfg.get< tag::nonblocking >()) dt();
498 : 3446 : }
499 : :
500 : : void
501 : 416 : ChoCG::setup( tk::real v )
502 : : // *****************************************************************************
503 : : // Start setup for solution
504 : : //! \param[in] v Total volume within user-specified box
505 : : // *****************************************************************************
506 : : {
507 : 416 : auto d = Disc();
508 : :
509 : : // Store user-defined box IC volume
510 : 416 : Disc()->Boxvol() = v;
511 : :
512 : : // Set initial conditions
513 : 416 : problems::initialize( d->Coord(), m_u, d->T(), d->BoxNodes() );
514 : :
515 : : // Query time history field output labels from all PDEs integrated
516 [ - + ]: 416 : if (!g_cfg.get< tag::histout >().empty()) {
517 : : std::vector< std::string > var
518 [ - - ][ - - ]: 0 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ][ - - ]
[ - - ]
519 : : auto ncomp = m_u.nprop();
520 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c)
521 [ - - ]: 0 : var.push_back( "c" + std::to_string(c-5) );
522 [ - - ]: 0 : d->histheader( std::move(var) );
523 : 0 : }
524 : :
525 : : // Compute finite element operators
526 : 416 : feop();
527 [ - - ][ - - ]: 416 : }
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
528 : :
529 : : void
530 : 416 : ChoCG::bnorm()
531 : : // *****************************************************************************
532 : : // Combine own and communicated portions of the boundary point normals
533 : : // *****************************************************************************
534 : : {
535 : 416 : const auto& lid = Disc()->Lid();
536 : :
537 : : // Combine own and communicated contributions to boundary point normals
538 [ + + ]: 1333 : for (const auto& [s,b] : m_bnormc) {
539 : : auto& bndnorm = m_bnorm[s];
540 [ + + ]: 8504 : for (const auto& [g,n] : b) {
541 : : auto& norm = bndnorm[g];
542 : 7587 : norm[0] += n[0];
543 : 7587 : norm[1] += n[1];
544 : 7587 : norm[2] += n[2];
545 : 7587 : norm[3] += n[3];
546 : : }
547 : : }
548 : 416 : tk::destroy( m_bnormc );
549 : :
550 : : // Divide summed point normals by the sum of the inverse distance squared
551 [ + + ]: 1426 : for (auto& [s,b] : m_bnorm) {
552 [ + + ]: 36384 : for (auto& [g,n] : b) {
553 : 35374 : n[0] /= n[3];
554 : 35374 : n[1] /= n[3];
555 : 35374 : n[2] /= n[3];
556 : : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
557 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
558 : : "Non-unit normal" );
559 : : }
560 : : }
561 : :
562 : : // Replace global->local ids associated to boundary point normals
563 : : decltype(m_bnorm) loc;
564 [ + + ]: 1426 : for (auto& [s,b] : m_bnorm) {
565 : : auto& bnd = loc[s];
566 [ + + ]: 36384 : for (auto&& [g,n] : b) {
567 : 35374 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
568 : : }
569 : : }
570 : : m_bnorm = std::move(loc);
571 : 416 : }
572 : :
573 : : void
574 [ + - ]: 416 : ChoCG::streamable()
575 : : // *****************************************************************************
576 : : // Convert integrals into streamable data structures
577 : : // *****************************************************************************
578 : : {
579 : : // Query surface integral output nodes
580 : : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
581 : : const auto& is = g_cfg.get< tag::integout >();
582 [ + - ]: 416 : std::set< int > outsets( begin(is), end(is) );
583 [ + + ]: 417 : for (auto s : outsets) {
584 : : auto m = m_bface.find(s);
585 [ + - ]: 1 : if (m != end(m_bface)) {
586 [ + - ]: 1 : auto& n = surfintnodes[ m->first ]; // associate set id
587 [ + - ][ + + ]: 847 : for (auto f : m->second) { // face ids on side set
588 : 846 : auto t = m_triinpoel.data() + f*3;
589 [ + - ]: 846 : n.push_back( t[0] ); // nodes on side set
590 [ + - ]: 846 : n.push_back( t[1] );
591 [ + - ]: 846 : n.push_back( t[2] );
592 : : }
593 : : }
594 : : }
595 [ + - ][ + + ]: 417 : for (auto& [s,n] : surfintnodes) tk::unique( n );
596 : : // Prepare surface integral data
597 : 416 : tk::destroy( m_surfint );
598 [ + - ]: 416 : const auto& gid = Disc()->Gid();
599 [ + + ]: 417 : for (auto&& [s,n] : surfintnodes) {
600 [ + - ]: 1 : auto& sint = m_surfint[s]; // associate set id
601 : 1 : auto& nodes = sint.first;
602 : 1 : auto& ndA = sint.second;
603 : : nodes = std::move(n);
604 [ + - ]: 1 : ndA.resize( nodes.size()*3 );
605 : : auto a = ndA.data();
606 [ + + ]: 426 : for (auto p : nodes) {
607 : : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
608 : 425 : a[0] = b[0]; // store ni * dA
609 : 425 : a[1] = b[1];
610 : 425 : a[2] = b[2];
611 : 425 : a += 3;
612 : : }
613 : : }
614 : 416 : tk::destroy( m_bndpoinint );
615 : :
616 : : // Generate domain superedges
617 [ + - ]: 416 : domsuped();
618 : :
619 : : // Prepare symmetry boundary condition data structures
620 : :
621 : : // Query symmetry BC nodes associated to side sets
622 : : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
623 [ + + ]: 718 : for (auto s : g_cfg.get< tag::bc_sym >()) {
624 : : auto k = m_bface.find(s);
625 [ + + ]: 302 : if (k != end(m_bface)) {
626 [ + - ]: 117 : auto& n = sym[ k->first ];
627 [ + - ][ + + ]: 3561 : for (auto f : k->second) {
628 [ + - ]: 3444 : const auto& t = m_triinpoel.data() + f*3;
629 : : n.insert( t[0] );
630 [ + - ]: 3444 : n.insert( t[1] );
631 [ + - ]: 3444 : n.insert( t[2] );
632 : : }
633 : : }
634 : : }
635 : :
636 : : // Generate unique set of symmetry BC nodes of all symmetryc BC side sets
637 : : std::set< std::size_t > symbcnodeset;
638 [ + + ]: 533 : for (const auto& [s,n] : sym) symbcnodeset.insert( begin(n), end(n) );
639 : :
640 : : // Generate symmetry BC data as streamable data structures
641 [ - + ]: 416 : tk::destroy( m_symbcnodes );
642 [ - + ]: 416 : tk::destroy( m_symbcnorms );
643 [ + + ]: 2853 : for (auto p : symbcnodeset) {
644 [ + + ]: 6093 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
645 : : auto m = m_bnorm.find( s );
646 [ + + ]: 3656 : if (m != end(m_bnorm)) {
647 : : auto r = m->second.find( p );
648 [ + + ]: 3138 : if (r != end(m->second)) {
649 [ + - ]: 2437 : m_symbcnodes.push_back( p );
650 [ + - ]: 2437 : m_symbcnorms.push_back( r->second[0] );
651 [ + - ]: 2437 : m_symbcnorms.push_back( r->second[1] );
652 [ + - ]: 2437 : m_symbcnorms.push_back( r->second[2] );
653 : : }
654 : : }
655 : : }
656 : : }
657 : 416 : tk::destroy( m_bnorm );
658 : :
659 : : // Prepare noslip boundary condition data structures
660 : :
661 : : // Query noslip BC nodes associated to side sets
662 : : std::unordered_map< int, std::unordered_set< std::size_t > > noslip;
663 [ + + ]: 505 : for (auto s : g_cfg.get< tag::bc_noslip >()) {
664 : : auto k = m_bface.find(s);
665 [ + + ]: 89 : if (k != end(m_bface)) {
666 [ + - ]: 83 : auto& n = noslip[ k->first ];
667 [ + - ][ + + ]: 5549 : for (auto f : k->second) {
668 [ + - ]: 5466 : const auto& t = m_triinpoel.data() + f*3;
669 : : n.insert( t[0] );
670 [ + - ]: 5466 : n.insert( t[1] );
671 [ + - ]: 5466 : n.insert( t[2] );
672 : : }
673 : : }
674 : : }
675 : :
676 : : // Generate unique set of noslip BC nodes of all noslip BC side sets
677 : : std::set< std::size_t > noslipbcnodeset;
678 [ + + ]: 499 : for (const auto& [s,n] : noslip) noslipbcnodeset.insert( begin(n), end(n) );
679 : :
680 : : // Generate noslip BC data as streamable data structures
681 [ - + ]: 416 : tk::destroy( m_noslipbcnodes );
682 [ + - ]: 416 : m_noslipbcnodes.insert( m_noslipbcnodes.end(),
683 : : begin(noslipbcnodeset), end(noslipbcnodeset) );
684 : 416 : }
685 : :
686 : : void
687 : 416 : ChoCG::domsuped()
688 : : // *****************************************************************************
689 : : // Generate superedge-groups for domain-edge loops
690 : : //! \see See Lohner, Sec. 15.1.6.2, An Introduction to Applied CFD Techniques,
691 : : //! Wiley, 2008.
692 : : // *****************************************************************************
693 : : {
694 : : Assert( !m_domedgeint.empty(), "No domain edges to group" );
695 : :
696 : : #ifndef NDEBUG
697 : : auto nedge = m_domedgeint.size();
698 : : #endif
699 : :
700 : 416 : const auto& inpoel = Disc()->Inpoel();
701 : 416 : const auto& lid = Disc()->Lid();
702 : 416 : const auto& gid = Disc()->Gid();
703 : :
704 : : tk::destroy( m_dsupedge[0] );
705 : : tk::destroy( m_dsupedge[1] );
706 : : tk::destroy( m_dsupedge[2] );
707 : :
708 : : tk::destroy( m_dsupint[0] );
709 : : tk::destroy( m_dsupint[1] );
710 : : tk::destroy( m_dsupint[2] );
711 : :
712 : : tk::UnsMesh::FaceSet untri;
713 [ + + ]: 78797 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
714 : : std::size_t N[4] = {
715 : 78381 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
716 [ + - ][ + + ]: 391905 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
717 : : }
718 : :
719 [ + + ]: 78797 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
720 : : std::size_t N[4] = {
721 : 78381 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
722 : : int f = 0;
723 : : tk::real sig[6];
724 [ + + ]: 548667 : decltype(m_domedgeint)::const_iterator d[6];
725 [ + + ]: 231185 : for (const auto& [p,q] : tk::lpoed) {
726 [ + + ]: 218992 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
727 [ + + ]: 298709 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
728 : 218992 : d[f] = m_domedgeint.find( ed );
729 [ + + ]: 218992 : if (d[f] == end(m_domedgeint)) break; else ++f;
730 : : }
731 [ + + ]: 78381 : if (f == 6) {
732 [ + - ]: 12193 : m_dsupedge[0].push_back( N[0] );
733 [ + - ]: 12193 : m_dsupedge[0].push_back( N[1] );
734 [ + - ]: 12193 : m_dsupedge[0].push_back( N[2] );
735 [ + - ]: 12193 : m_dsupedge[0].push_back( N[3] );
736 [ + + ]: 60965 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
737 [ + + ]: 85351 : for (int ed=0; ed<6; ++ed) {
738 : : const auto& ded = d[ed]->second;
739 [ + - ]: 73158 : m_dsupint[0].push_back( sig[ed] * ded[0] );
740 [ + - ]: 73158 : m_dsupint[0].push_back( sig[ed] * ded[1] );
741 [ + - ][ + - ]: 73158 : m_dsupint[0].push_back( sig[ed] * ded[2] );
742 [ + - ]: 73158 : m_dsupint[0].push_back( ded[3] );
743 [ + - ]: 73158 : m_dsupint[0].push_back( ded[4] );
744 : 73158 : m_domedgeint.erase( d[ed] );
745 : : }
746 : : }
747 : : }
748 : :
749 [ + + ]: 138246 : for (const auto& N : untri) {
750 : : int f = 0;
751 : : tk::real sig[3];
752 [ + + ]: 551320 : decltype(m_domedgeint)::const_iterator d[3];
753 [ + + ]: 226975 : for (const auto& [p,q] : tk::lpoet) {
754 [ + + ]: 215392 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
755 [ + + ]: 328676 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
756 : 215392 : d[f] = m_domedgeint.find( ed );
757 [ + + ]: 215392 : if (d[f] == end(m_domedgeint)) break; else ++f;
758 : : }
759 [ + + ]: 137830 : if (f == 3) {
760 [ + - ]: 11583 : m_dsupedge[1].push_back( N[0] );
761 [ + - ]: 11583 : m_dsupedge[1].push_back( N[1] );
762 [ + - ]: 11583 : m_dsupedge[1].push_back( N[2] );
763 [ + + ]: 46332 : for (int ed=0; ed<3; ++ed) {
764 : : const auto& ded = d[ed]->second;
765 [ + - ]: 34749 : m_dsupint[1].push_back( sig[ed] * ded[0] );
766 [ + - ]: 34749 : m_dsupint[1].push_back( sig[ed] * ded[1] );
767 [ + - ][ + - ]: 34749 : m_dsupint[1].push_back( sig[ed] * ded[2] );
768 [ + - ]: 34749 : m_dsupint[1].push_back( ded[3] );
769 [ + - ]: 34749 : m_dsupint[1].push_back( ded[4] );
770 : 34749 : m_domedgeint.erase( d[ed] );
771 : : }
772 : : }
773 : : }
774 : :
775 [ + - ]: 416 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
776 [ + - ]: 416 : m_dsupint[2].resize( m_domedgeint.size()*5 );
777 : : std::size_t k = 0;
778 [ + + ]: 32321 : for (const auto& [ed,d] : m_domedgeint) {
779 : 31905 : auto e = m_dsupedge[2].data() + k*2;
780 : 31905 : e[0] = tk::cref_find( lid, ed[0] );
781 : 31905 : e[1] = tk::cref_find( lid, ed[1] );
782 : 31905 : auto i = m_dsupint[2].data() + k*5;
783 : 31905 : i[0] = d[0];
784 : 31905 : i[1] = d[1];
785 : 31905 : i[2] = d[2];
786 : 31905 : i[3] = d[3];
787 : 31905 : i[4] = d[4];
788 : 31905 : ++k;
789 : : }
790 : :
791 : 416 : tk::destroy( m_domedgeint );
792 : :
793 : : //std::cout << std::setprecision(2)
794 : : // << "superedges: ntet:" << m_dsupedge[0].size()/4 << "(nedge:"
795 : : // << m_dsupedge[0].size()/4*6 << ","
796 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[0].size()/4*6 ) /
797 : : // static_cast< tk::real >( nedge )
798 : : // << "%) + ntri:" << m_dsupedge[1].size()/3
799 : : // << "(nedge:" << m_dsupedge[1].size() << ","
800 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[1].size() ) /
801 : : // static_cast< tk::real >( nedge )
802 : : // << "%) + nedge:"
803 : : // << m_dsupedge[2].size()/2 << "("
804 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[2].size()/2 ) /
805 : : // static_cast< tk::real >( nedge )
806 : : // << "%) = " << m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
807 : : // m_dsupedge[2].size()/2 << " of "<< nedge << " total edges\n";
808 : :
809 : : Assert( m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
810 : : m_dsupedge[2].size()/2 == nedge,
811 : : "Not all edges accounted for in superedge groups" );
812 : 416 : }
813 : :
814 : : void
815 : : // cppcheck-suppress unusedFunction
816 : 416 : ChoCG::merge()
817 : : // *****************************************************************************
818 : : // Combine own and communicated portions of the integrals
819 : : // *****************************************************************************
820 : : {
821 : : // Combine own and communicated contributions to boundary point normals
822 : 416 : bnorm();
823 : :
824 : : // Convert integrals into streamable data structures
825 : 416 : streamable();
826 : :
827 : : // Enforce boundary conditions on initial conditions
828 : 416 : BC( m_u, Disc()->T() );
829 : :
830 : : // Start measuring initial div-free time
831 : 416 : m_timer.emplace_back();
832 : :
833 : : // Compute initial momentum flux
834 [ + - ]: 416 : thisProxy[ thisIndex ].wait4div();
835 [ + - ]: 416 : thisProxy[ thisIndex ].wait4sgrad();
836 : 416 : div( m_u );
837 : 416 : }
838 : :
839 : : void
840 : 14188 : ChoCG::fingrad( tk::Fields& grad,
841 : : std::unordered_map< std::size_t, std::vector< tk::real > >& gradc )
842 : : // *****************************************************************************
843 : : // Finalize computing gradient
844 : : //! \param[in,out] grad Gradient to finalize
845 : : //! \param[in,out] gradc Gradient communication buffer to finalize
846 : : // *****************************************************************************
847 : : {
848 : 14188 : auto d = Disc();
849 : : const auto lid = d->Lid();
850 : :
851 : : // Combine own and communicated contributions
852 [ + + ]: 309292 : for (const auto& [g,r] : gradc) {
853 : 295104 : auto i = tk::cref_find( lid, g );
854 [ + + ]: 1561737 : for (std::size_t c=0; c<r.size(); ++c) grad(i,c) += r[c];
855 : : }
856 : 14188 : tk::destroy(gradc);
857 : :
858 : : // Divide weak result by nodal volume
859 : : const auto& vol = d->Vol();
860 [ + + ]: 1595435 : for (std::size_t p=0; p<grad.nunk(); ++p) {
861 [ + + ]: 9979561 : for (std::size_t c=0; c<grad.nprop(); ++c) {
862 : 8398314 : grad(p,c) /= vol[p];
863 : : }
864 : : }
865 : 14188 : }
866 : :
867 : : void
868 : 5450 : ChoCG::div( const tk::Fields& u )
869 : : // *****************************************************************************
870 : : // Start computing divergence
871 : : // \para[in] u Vector field whose divergence to compute
872 : : // *****************************************************************************
873 : : {
874 : 5450 : auto d = Disc();
875 : : const auto lid = d->Lid();
876 : :
877 : : // Finalize momentum flux communications if needed
878 [ + + ]: 5450 : if (m_np == 1) {
879 [ + - ]: 384 : fingrad( m_flux, m_fluxc );
880 [ + - ]: 384 : physics::symbc( m_flux, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
881 : : }
882 : :
883 : : // Compute divergence
884 : : std::fill( begin(m_div), end(m_div), 0.0 );
885 : 5450 : chorin::div( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel,
886 [ + - ]: 5450 : d->Dt(), m_pr, m_pgrad, u, m_div, m_np>1 );
887 : :
888 : : // Communicate velocity divergence to other chares on chare-boundary
889 [ + + ]: 5450 : if (d->NodeCommMap().empty()) {
890 [ + - ]: 259 : comdiv_complete();
891 : : } else {
892 [ + + ]: 53813 : for (const auto& [c,n] : d->NodeCommMap()) {
893 : : decltype(m_divc) exp;
894 [ + - ][ + + ]: 275120 : for (auto g : n) exp[g] = m_div[ tk::cref_find(lid,g) ];
895 [ + - ][ + - ]: 97244 : thisProxy[c].comdiv( exp );
896 : : }
897 : : }
898 [ + - ]: 5450 : owndiv_complete();
899 : 5450 : }
900 : :
901 : : void
902 : 48622 : ChoCG::comdiv( const std::unordered_map< std::size_t, tk::real >& indiv )
903 : : // *****************************************************************************
904 : : // Receive contributions to velocity divergence on chare-boundaries
905 : : //! \param[in] indiv Partial contributions of velocity divergence to
906 : : //! chare-boundary nodes. Key: global mesh node IDs, value: contribution.
907 : : //! \details This function receives contributions to m_div, which stores the
908 : : //! velocity divergence at mesh nodes. While m_div stores own contributions,
909 : : //! m_divc collects the neighbor chare contributions during communication.
910 : : //! This way work on m_div and m_divc is overlapped.
911 : : // *****************************************************************************
912 : : {
913 [ + + ]: 275120 : for (const auto& [g,r] : indiv) m_divc[g] += r;
914 : :
915 : : // When we have heard from all chares we communicate with, this chare is done
916 [ + + ]: 48622 : if (++m_ndiv == Disc()->NodeCommMap().size()) {
917 : 5191 : m_ndiv = 0;
918 : 5191 : comdiv_complete();
919 : : }
920 : 48622 : }
921 : :
922 : : void
923 : 3704 : ChoCG::velgrad()
924 : : // *****************************************************************************
925 : : // Start computing velocity gradient
926 : : // *****************************************************************************
927 : : {
928 : 3704 : auto d = Disc();
929 : :
930 : : // Compute momentum flux
931 : : m_vgrad.fill( 0.0 );
932 : 3704 : chorin::vgrad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad );
933 : :
934 : : // Communicate velocity divergence to other chares on chare-boundary
935 : : const auto& lid = d->Lid();
936 [ + + ]: 3704 : if (d->NodeCommMap().empty()) {
937 : 394 : comvgrad_complete();
938 : : } else {
939 [ + + ]: 42108 : for (const auto& [c,n] : d->NodeCommMap()) {
940 : : decltype(m_vgradc) exp;
941 [ + - ][ + + ]: 188180 : for (auto g : n) exp[g] = m_vgrad[ tk::cref_find(lid,g) ];
942 [ + - ][ + - ]: 77596 : thisProxy[c].comvgrad( exp );
943 : : }
944 : : }
945 : 3704 : ownvgrad_complete();
946 : 3704 : }
947 : :
948 : : void
949 : 38798 : ChoCG::comvgrad(
950 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
951 : : // *****************************************************************************
952 : : // Receive contributions to velocity gradients on chare-boundaries
953 : : //! \param[in] ingrad Partial contributions of momentum flux to
954 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
955 : : //! \details This function receives contributions to m_vgrad, which stores the
956 : : //! velocity gradients at mesh nodes. While m_vgrad stores own contributions,
957 : : //! m_vgradc collects the neighbor chare contributions during communication.
958 : : //! This way work on m_vgrad and m_vgradc is overlapped.
959 : : // *****************************************************************************
960 : : {
961 : : using tk::operator+=;
962 [ + + ]: 337562 : for (const auto& [g,r] : ingrad) m_vgradc[g] += r;
963 : :
964 : : // When we have heard from all chares we communicate with, this chare is done
965 [ + + ]: 38798 : if (++m_nvgrad == Disc()->NodeCommMap().size()) {
966 : 3310 : m_nvgrad = 0;
967 : 3310 : comvgrad_complete();
968 : : }
969 : 38798 : }
970 : :
971 : : void
972 : 384 : ChoCG::flux()
973 : : // *****************************************************************************
974 : : // Start computing momentum flux
975 : : // *****************************************************************************
976 : : {
977 : 384 : auto d = Disc();
978 : :
979 : : // Finalize computing velocity gradients
980 : 384 : fingrad( m_vgrad, m_vgradc );
981 : :
982 : : // Compute momentum flux
983 : : m_flux.fill( 0.0 );
984 : 384 : chorin::flux( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad,
985 : 384 : m_flux );
986 : :
987 : : // Communicate velocity divergence to other chares on chare-boundary
988 : : const auto& lid = d->Lid();
989 [ + + ]: 384 : if (d->NodeCommMap().empty()) {
990 : 14 : comflux_complete();
991 : : } else {
992 [ + + ]: 4168 : for (const auto& [c,n] : d->NodeCommMap()) {
993 : : decltype(m_fluxc) exp;
994 [ + - ][ + + ]: 20320 : for (auto g : n) exp[g] = m_flux[ tk::cref_find(lid,g) ];
995 [ + - ][ + - ]: 7596 : thisProxy[c].comflux( exp );
996 : : }
997 : : }
998 : 384 : ownflux_complete();
999 : 384 : }
1000 : :
1001 : : void
1002 : 3798 : ChoCG::comflux(
1003 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& influx )
1004 : : // *****************************************************************************
1005 : : // Receive contributions to momentum flux on chare-boundaries
1006 : : //! \param[in] influx Partial contributions of momentum flux to
1007 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
1008 : : //! \details This function receives contributions to m_flux, which stores the
1009 : : //! momentum flux at mesh nodes. While m_flux stores own contributions,
1010 : : //! m_fluxc collects the neighbor chare contributions during communication.
1011 : : //! This way work on m_flux and m_fluxc is overlapped.
1012 : : // *****************************************************************************
1013 : : {
1014 : : using tk::operator+=;
1015 [ + + ]: 36842 : for (const auto& [g,r] : influx) m_fluxc[g] += r;
1016 : :
1017 : : // When we have heard from all chares we communicate with, this chare is done
1018 [ + + ]: 3798 : if (++m_nflux == Disc()->NodeCommMap().size()) {
1019 : 370 : m_nflux = 0;
1020 : 370 : comflux_complete();
1021 : : }
1022 : 3798 : }
1023 : :
1024 : : void
1025 : 5450 : ChoCG::pinit()
1026 : : // *****************************************************************************
1027 : : // Initialize Poisson solve
1028 : : // *****************************************************************************
1029 : : {
1030 : 5450 : auto d = Disc();
1031 : : const auto lid = d->Lid();
1032 : : const auto& coord = d->Coord();
1033 : : const auto& x = coord[0];
1034 : : const auto& y = coord[1];
1035 : : const auto& z = coord[2];
1036 : :
1037 : : // Combine own and communicated contributions to velocity divergence
1038 [ + + ]: 126134 : for (const auto& [g,r] : m_divc) m_div[ tk::cref_find(lid,g) ] += r;
1039 : 5450 : tk::destroy(m_divc);
1040 : :
1041 : : // Divide Poisson rhs by dt if solving for time-increment
1042 [ + + ][ + + ]: 495420 : if (m_np > 1) for (auto& div : m_div) div /= d->Dt();
1043 : :
1044 : : // Configure Poisson BCs
1045 : : std::unordered_map< std::size_t,
1046 : : std::vector< std::pair< int, tk::real > > > dirbc;
1047 : : std::vector< tk::real > neubc;
1048 : :
1049 [ + - ]: 5450 : if (m_np < 3) { // only during startup and first time step
1050 : : // Configure Dirichlet BCs
1051 [ + + ]: 5450 : if (!g_cfg.get< tag::pre_bc_dir >().empty()) {
1052 [ + - ]: 752 : auto ic = problems::PRESSURE_IC();
1053 : : std::size_t nmask = 1 + 1;
1054 : : Assert( m_dirbcmaskp.size() % nmask == 0, "Size mismatch" );
1055 [ + + ]: 9194 : for (std::size_t i=0; i<m_dirbcmaskp.size()/nmask; ++i) {
1056 [ + + ]: 8442 : auto p = m_dirbcmaskp[i*nmask+0]; // local node id
1057 : 8442 : auto mask = m_dirbcmaskp[i*nmask+1];
1058 [ + + ]: 8442 : if (mask == 1) { // mask == 1: IC value
1059 [ + + ][ - + ]: 5620 : auto val = m_np>1 ? 0.0 : ic( x[p], y[p], z[p] );
1060 [ + - ]: 7380 : dirbc[p] = {{ { 1, val } }};
1061 [ + - ][ + - ]: 4752 : } else if (mask == 2 && !m_dirbcvalp.empty()) { // mask == 2: BC value
1062 [ + + ]: 4752 : auto val = m_np>1 ? 0.0 : m_dirbcvalp[i*nmask+1];
1063 [ + - ][ - - ]: 9504 : dirbc[p] = {{ { 1, val } }};
1064 : : }
1065 : : }
1066 : : }
1067 : :
1068 : : // Configure Neumann BCs
1069 [ + - ]: 5450 : auto pg = problems::PRESSURE_GRAD();
1070 [ + + ]: 5450 : if (pg) {
1071 : : // Collect Neumann BC elements
1072 [ + - ][ - - ]: 2 : std::vector< std::uint8_t > besym( m_triinpoel.size(), 0 );
1073 [ + + ]: 8 : for (auto s : g_cfg.get< tag::pre_bc_sym >()) {
1074 : : auto k = m_bface.find(s);
1075 [ + + ][ + + ]: 210 : if (k != end(m_bface)) for (auto f : k->second) besym[f] = 1;
1076 : : }
1077 : : // Setup Neumann BCs
1078 [ + - ][ - - ]: 2 : neubc.resize( x.size(), 0.0 );
1079 [ + + ]: 410 : for (std::size_t e=0; e<m_triinpoel.size()/3; ++e) {
1080 [ + + ]: 408 : if (besym[e]) {
1081 : 204 : const auto N = m_triinpoel.data() + e*3;
1082 : : tk::real n[3];
1083 [ - + ]: 204 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
1084 [ - + ]: 204 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
1085 : : n[0], n[1], n[2] );
1086 : 204 : auto g = pg( x[N[0]], y[N[0]], z[N[0]] );
1087 [ - + ]: 204 : neubc[ N[0] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1088 [ - + ]: 204 : g = pg( x[N[1]], y[N[1]], z[N[1]] );
1089 [ - + ]: 204 : neubc[ N[1] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1090 [ - + ]: 204 : g = pg( x[N[2]], y[N[2]], z[N[2]] );
1091 : 204 : neubc[ N[2] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1092 : : }
1093 : : }
1094 : : }
1095 : :
1096 : : // Set hydrostat
1097 : 5450 : auto h = g_cfg.get< tag::pre_hydrostat >();
1098 [ + + ]: 5450 : if (h != std::numeric_limits< uint64_t >::max()) {
1099 : : auto pi = lid.find( h );
1100 [ + + ]: 1194 : if (pi != end(lid)) {
1101 : 146 : auto p = pi->second;
1102 [ + - ]: 146 : auto ic = problems::PRESSURE_IC();
1103 [ + + ][ - + ]: 162 : auto val = m_np>1 ? 0.0 : ic( x[p], y[p], z[p] );
1104 : : auto& b = dirbc[p];
1105 [ + - ][ + - ]: 292 : if (b.empty()) b = {{ { 1, val }} };
[ - - ]
1106 : : }
1107 : : }
1108 : :
1109 : : // Configure right hand side
1110 [ + - ]: 5450 : auto pr = problems::PRESSURE_RHS();
1111 [ + + ]: 5450 : if (pr) {
1112 : : const auto& vol = d->Vol();
1113 [ + + ]: 2163 : for (std::size_t i=0; i<x.size(); ++i) {
1114 [ - + ]: 4262 : m_div[i] = pr( x[i], y[i], z[i] ) * vol[i];
1115 : : }
1116 : : }
1117 : : }
1118 : :
1119 : : // Initialize Poisson solve setting ICs and BCs
1120 : 5450 : m_cgpre[ thisIndex ].ckLocal()->
1121 [ + - ][ + - ]: 21800 : init( {}, m_div, neubc, dirbc, m_np<3, g_cfg.get< tag::pre_pc >(),
[ + + ][ - - ]
1122 [ + - ][ + - ]: 16350 : CkCallback( CkIndex_ChoCG::psolve(), thisProxy[thisIndex] ) );
[ - + ][ - - ]
1123 : 5450 : }
1124 : :
1125 : : void
1126 : 5450 : ChoCG::psolve()
1127 : : // *****************************************************************************
1128 : : // Solve Poisson equation
1129 : : // *****************************************************************************
1130 : : {
1131 : 5450 : auto iter = g_cfg.get< tag::pre_iter >();
1132 : 5450 : auto tol = g_cfg.get< tag::pre_tol >();
1133 : 5450 : auto verbose = g_cfg.get< tag::pre_verbose >();
1134 : :
1135 [ + + ]: 5450 : auto c = m_np != 1 ?
1136 : 5066 : CkCallback( CkIndex_ChoCG::sgrad(), thisProxy[thisIndex] ) :
1137 [ + - ][ + - ]: 11284 : CkCallback( CkIndex_ChoCG::psolved(), thisProxy[thisIndex] );
[ + - ][ + - ]
[ + + ][ - - ]
1138 : :
1139 [ + - ][ + - ]: 16350 : m_cgpre[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose, c );
[ - + ][ - - ]
1140 : 5450 : }
1141 : :
1142 : : void
1143 : 5066 : ChoCG::sgrad()
1144 : : // *****************************************************************************
1145 : : // Compute recent conjugate gradients solution gradient
1146 : : // *****************************************************************************
1147 : : {
1148 : 5066 : auto d = Disc();
1149 : :
1150 [ + - ]: 15198 : auto sol = m_cgpre[ thisIndex ].ckLocal()->solution();
1151 : : m_sgrad.fill( 0.0 );
1152 [ + - ]: 5066 : chorin::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, sol, m_sgrad );
1153 : :
1154 : : // Send gradient contributions to neighbor chares
1155 [ + + ]: 5066 : if (d->NodeCommMap().empty()) {
1156 [ + - ]: 245 : comsgrad_complete();
1157 : : } else {
1158 : : const auto& lid = d->Lid();
1159 [ + + ]: 49645 : for (const auto& [c,n] : d->NodeCommMap()) {
1160 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1161 [ + - ][ + + ]: 254800 : for (auto g : n) exp[g] = m_sgrad[ tk::cref_find(lid,g) ];
1162 [ + - ][ + - ]: 89648 : thisProxy[c].comsgrad( exp );
1163 : : }
1164 : : }
1165 [ + - ]: 5066 : ownsgrad_complete();
1166 : 5066 : }
1167 : :
1168 : : void
1169 : 44824 : 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 [ + + ]: 464776 : 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 [ + + ]: 44824 : if (++m_nsgrad == Disc()->NodeCommMap().size()) {
1186 : 4821 : m_nsgrad = 0;
1187 : 4821 : comsgrad_complete();
1188 : : }
1189 : 44824 : }
1190 : :
1191 : : void
1192 : 5450 : ChoCG::psolved()
1193 : : // *****************************************************************************
1194 : : // Continue setup after Poisson solve and gradient computation
1195 : : // *****************************************************************************
1196 : : {
1197 : 5450 : auto d = Disc();
1198 : :
1199 [ + + ][ + - ]: 6676 : if (thisIndex == 0) d->pit( m_cgpre[ thisIndex ].ckLocal()->it() );
1200 : :
1201 [ + + ]: 5450 : if (m_np != 1) {
1202 : : // Finalize gradient communications
1203 : 5066 : fingrad( m_sgrad, m_sgradc );
1204 : : // Project velocity to divergence-free subspace
1205 [ + + ]: 5066 : auto dt = m_np > 1 ? d->Dt() : 1.0;
1206 [ + + ]: 527046 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1207 : 521980 : m_u(i,0) -= dt * m_sgrad(i,0);
1208 : 521980 : m_u(i,1) -= dt * m_sgrad(i,1);
1209 : 521980 : m_u(i,2) -= dt * m_sgrad(i,2);
1210 : : }
1211 : : // Enforce boundary conditions
1212 : 5066 : BC( m_u, d->T() + d->Dt() );
1213 : : }
1214 : :
1215 [ + + ]: 5450 : if (d->Initial()) {
1216 : :
1217 [ + + ]: 800 : if (g_cfg.get< tag::nstep >() == 1) { // test first Poisson solve only
1218 : :
1219 [ + - ]: 64 : m_pr = m_cgpre[ thisIndex ].ckLocal()->solution();
1220 [ + - ]: 32 : thisProxy[ thisIndex ].wait4step();
1221 [ + - ][ + - ]: 96 : writeFields( CkCallback(CkIndex_ChoCG::diag(), thisProxy[thisIndex]) );
1222 : :
1223 : : } else {
1224 : :
1225 [ + + ]: 768 : if (++m_np < 2) {
1226 : : // Compute momentum flux for initial pressure-Poisson rhs
1227 [ + - ]: 384 : thisProxy[ thisIndex ].wait4vgrad();
1228 [ + - ]: 384 : thisProxy[ thisIndex ].wait4flux();
1229 [ + - ]: 384 : thisProxy[ thisIndex ].wait4div();
1230 : 384 : velgrad();
1231 : : } else {
1232 [ + + ]: 384 : if (thisIndex == 0) {
1233 : : tk::Print() << "Initial div-free time: " << m_timer[0].dsec()
1234 : : << " sec\n";
1235 : : }
1236 : : // Assign initial pressure and compute its gradient
1237 [ + - ]: 768 : m_pr = m_cgpre[ thisIndex ].ckLocal()->solution();
1238 : 384 : pgrad();
1239 : : }
1240 : :
1241 : : }
1242 : :
1243 : : } else {
1244 : :
1245 : : // Update pressure and compute its gradient
1246 : : using tk::operator+=;
1247 [ + - ]: 9300 : m_pr += m_cgpre[ thisIndex ].ckLocal()->solution();
1248 : 4650 : pgrad();
1249 : :
1250 : : }
1251 : 5450 : }
1252 : :
1253 : : void
1254 : 5034 : ChoCG::pgrad()
1255 : : // *****************************************************************************
1256 : : // Compute pressure gradient
1257 : : // *****************************************************************************
1258 : : {
1259 : 5034 : auto d = Disc();
1260 : :
1261 [ + - ]: 10068 : thisProxy[ thisIndex ].wait4pgrad();
1262 : :
1263 : : m_pgrad.fill( 0.0 );
1264 : 5034 : chorin::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_pr, m_pgrad );
1265 : :
1266 : : // Send gradient contributions to neighbor chares
1267 [ + + ]: 5034 : if (d->NodeCommMap().empty()) {
1268 : 244 : compgrad_complete();
1269 : : } else {
1270 : : const auto& lid = d->Lid();
1271 [ + + ]: 49428 : for (const auto& [c,n] : d->NodeCommMap()) {
1272 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1273 [ + - ][ + + ]: 253380 : for (auto g : n) exp[g] = m_pgrad[ tk::cref_find(lid,g) ];
1274 [ + - ][ + - ]: 89276 : thisProxy[c].compgrad( exp );
1275 : : }
1276 : : }
1277 : 5034 : ownpgrad_complete();
1278 : 5034 : }
1279 : :
1280 : : void
1281 : 44638 : 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 [ + + ]: 462122 : 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 [ + + ]: 44638 : if (++m_npgrad == Disc()->NodeCommMap().size()) {
1298 : 4790 : m_npgrad = 0;
1299 : 4790 : compgrad_complete();
1300 : : }
1301 : 44638 : }
1302 : :
1303 : : void
1304 : 5034 : ChoCG::finpgrad()
1305 : : // *****************************************************************************
1306 : : // Compute pressure gradient
1307 : : // *****************************************************************************
1308 : : {
1309 : 5034 : auto d = Disc();
1310 : :
1311 : : // Finalize pressure gradient communications
1312 : 5034 : fingrad( m_pgrad, m_pgradc );
1313 : :
1314 [ + + ]: 5034 : if (d->Initial()) {
1315 [ + - ][ + - ]: 1152 : writeFields( CkCallback(CkIndex_ChoCG::start(), thisProxy[thisIndex]) );
1316 : : } else {
1317 : 4650 : diag();
1318 : : }
1319 : 5034 : }
1320 : :
1321 : : void
1322 : 384 : ChoCG::start()
1323 : : // *****************************************************************************
1324 : : // Start time stepping
1325 : : // *****************************************************************************
1326 : : {
1327 : : // Set flag that indicates that we are now during time stepping
1328 : 384 : Disc()->Initial( 0 );
1329 : : // Start timer measuring time stepping wall clock time
1330 : 384 : Disc()->Timer().zero();
1331 : : // Zero grind-timer
1332 : 384 : Disc()->grindZero();
1333 : : // Continue to first time step
1334 : 384 : dt();
1335 : 384 : }
1336 : :
1337 : : void
1338 : 12392 : ChoCG::BC( tk::Fields& u, tk::real t )
1339 : : // *****************************************************************************
1340 : : // Apply boundary conditions
1341 : : //! \param[in,out] u Solution to apply BCs to
1342 : : //! \param[in] t Physical time
1343 : : // *****************************************************************************
1344 : : {
1345 : 12392 : auto d = Disc();
1346 : :
1347 : 12392 : physics::dirbc( u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
1348 : 12392 : physics::symbc( u, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
1349 : 12392 : physics::noslipbc( u, m_noslipbcnodes, /*pos=*/0 );
1350 : 12392 : }
1351 : :
1352 : : void
1353 : 4650 : ChoCG::dt()
1354 : : // *****************************************************************************
1355 : : // Compute time step size
1356 : : // *****************************************************************************
1357 : : {
1358 : 4650 : auto d = Disc();
1359 : : const auto& vol = d->Vol();
1360 : :
1361 : 4650 : tk::real mindt = std::numeric_limits< tk::real >::max();
1362 [ + + ]: 4650 : auto const_dt = g_cfg.get< tag::dt >();
1363 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1364 : :
1365 : : // use constant dt if configured
1366 [ + + ]: 4650 : if (std::abs(const_dt) > eps) {
1367 : :
1368 : : // cppcheck-suppress redundantInitialization
1369 : 2920 : mindt = const_dt;
1370 : :
1371 : : } else {
1372 : :
1373 : 1730 : auto cfl = g_cfg.get< tag::cfl >();
1374 : : auto large = std::numeric_limits< tk::real >::max();
1375 : 1730 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1376 [ - + ]: 1730 : auto dif = g_cfg.get< tag::mat_dyn_diffusivity >();
1377 : 1730 : dif = std::max( mu, dif );
1378 : :
1379 [ + + ]: 414720 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1380 : 412990 : auto u = m_u(i,0);
1381 : 412990 : auto v = m_u(i,1);
1382 : 412990 : auto w = m_u(i,2);
1383 : 412990 : auto vel = std::sqrt( u*u + v*v + w*w );
1384 : 412990 : auto L = std::cbrt( vol[i] );
1385 [ + + ][ + + ]: 501950 : auto euler_dt = L / std::max( vel, 1.0e-8 );
1386 : 412990 : mindt = std::min( mindt, euler_dt );
1387 [ + + ][ + + ]: 412990 : auto dif_dt = dif > eps ? L * L / dif : large;
1388 : 412990 : mindt = std::min( mindt, dif_dt );
1389 : : }
1390 [ + + ]: 1730 : mindt *= cfl;
1391 : :
1392 [ + + ]: 1730 : if (d->T() > g_cfg.get< tag::freezetime >()) {
1393 : 1634 : m_freezeflow = g_cfg.get< tag::freezeflow >();
1394 : : }
1395 : 1730 : mindt *= m_freezeflow;
1396 : :
1397 : : }
1398 : :
1399 : : // Actiavate SDAG waits for next time step stage
1400 [ + - ]: 4650 : thisProxy[ thisIndex ].wait4rhs();
1401 [ + - ]: 4650 : thisProxy[ thisIndex ].wait4div();
1402 [ + - ]: 4650 : thisProxy[ thisIndex ].wait4sgrad();
1403 [ + - ]: 4650 : thisProxy[ thisIndex ].wait4step();
1404 : :
1405 : : // Contribute to minimum dt across all chares and advance to next step
1406 [ + - ]: 4650 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
1407 : 4650 : CkCallback(CkReductionTarget(ChoCG,advance), thisProxy) );
1408 : 4650 : }
1409 : :
1410 : : void
1411 : 4650 : ChoCG::advance( tk::real newdt )
1412 : : // *****************************************************************************
1413 : : // Advance equations to next time step
1414 : : //! \param[in] newdt The smallest dt across the whole problem
1415 : : // *****************************************************************************
1416 : : {
1417 : : // Detect blowup
1418 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1419 [ - + ]: 4650 : if (newdt < eps) m_finished = 1;
1420 : :
1421 : : // Set new time step size
1422 : 4650 : Disc()->setdt( newdt );
1423 : :
1424 : : // Compute lhs and rhs of transport equations
1425 : 4650 : lhs();
1426 : 4650 : rhs();
1427 : 4650 : }
1428 : :
1429 : : void
1430 : 4650 : ChoCG::lhs()
1431 : : // *****************************************************************************
1432 : : // Fill lhs matrix of transport equations
1433 : : // *****************************************************************************
1434 : : {
1435 : 4650 : auto theta = g_cfg.get< tag::theta >();
1436 [ + + ]: 4650 : if (theta < std::numeric_limits< tk::real >::epsilon()) return;
1437 : :
1438 : 40 : auto d = Disc();
1439 : : const auto& inpoel = d->Inpoel();
1440 : : const auto& coord = d->Coord();
1441 : : const auto& X = coord[0];
1442 : : const auto& Y = coord[1];
1443 : : const auto& Z = coord[2];
1444 : : const auto ncomp = m_u.nprop();
1445 : 40 : const auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1446 : :
1447 : : auto dt = d->Dt();
1448 : 40 : auto& A = Lhs();
1449 : 40 : A.zero();
1450 : :
1451 [ + + ]: 60040 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1452 : 60000 : const auto N = inpoel.data() + e*4;
1453 : : const std::array< tk::real, 3 >
1454 : 60000 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
1455 : 60000 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
1456 : 60000 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
1457 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
1458 : : Assert( J > 0, "Element Jacobian non-positive" );
1459 : : std::array< std::array< tk::real, 3 >, 4 > grad;
1460 : 60000 : grad[1] = tk::cross( ca, da );
1461 : 60000 : grad[2] = tk::cross( da, ba );
1462 : 60000 : grad[3] = tk::cross( ba, ca );
1463 [ + + ]: 240000 : for (std::size_t i=0; i<3; ++i)
1464 : 180000 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
1465 [ + + ]: 300000 : for (std::size_t a=0; a<4; ++a) {
1466 [ + + ]: 1200000 : for (std::size_t b=0; b<4; ++b) {
1467 [ + + ]: 960000 : auto v = J/dt/120.0 * ((a == b) ? 2.0 : 1.0);
1468 : 960000 : v += theta * mu * tk::dot(grad[a],grad[b]) / J / 6.0;
1469 [ + + ]: 3840000 : for (std::size_t c=0; c<ncomp; ++c) A(N[a],N[b],c) -= v;
1470 : : }
1471 : : }
1472 : : }
1473 : : }
1474 : :
1475 : : void
1476 : 6910 : ChoCG::rhs()
1477 : : // *****************************************************************************
1478 : : // Compute right-hand side of transport equations
1479 : : // *****************************************************************************
1480 : : {
1481 : 6910 : auto d = Disc();
1482 : : const auto& lid = d->Lid();
1483 : :
1484 : : // Compute own portion of right-hand side for all equations
1485 : 6910 : chorin::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, d->V(), d->T(),
1486 : 6910 : m_pr, m_u, m_vgrad, m_rhs );
1487 : :
1488 : : // Communicate rhs to other chares on chare-boundary
1489 [ + + ]: 6910 : if (d->NodeCommMap().empty()) {
1490 : 570 : comrhs_complete();
1491 : : } else {
1492 [ + + ]: 57020 : for (const auto& [c,n] : d->NodeCommMap()) {
1493 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1494 [ + - ][ + + ]: 331700 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
1495 [ + - ][ + - ]: 101360 : thisProxy[c].comrhs( exp );
1496 : : }
1497 : : }
1498 : 6910 : ownrhs_complete();
1499 : 6910 : }
1500 : :
1501 : : void
1502 : 50680 : ChoCG::comrhs(
1503 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1504 : : // *****************************************************************************
1505 : : // Receive contributions to right-hand side vector on chare-boundaries
1506 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1507 : : //! global mesh node IDs, value: contributions for all scalar components.
1508 : : //! \details This function receives contributions to m_rhs, which stores the
1509 : : //! right hand side vector at mesh nodes. While m_rhs stores own
1510 : : //! contributions, m_rhsc collects the neighbor chare contributions during
1511 : : //! communication. This way work on m_rhs and m_rhsc is overlapped.
1512 : : // *****************************************************************************
1513 : : {
1514 : : using tk::operator+=;
1515 [ + + ]: 612720 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
1516 : :
1517 : : // When we have heard from all chares we communicate with, this chare is done
1518 [ + + ]: 50680 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
1519 : 6340 : m_nrhs = 0;
1520 : 6340 : comrhs_complete();
1521 : : }
1522 : 50680 : }
1523 : :
1524 : : void
1525 : : // cppcheck-suppress unusedFunction
1526 : 6910 : ChoCG::solve()
1527 : : // *****************************************************************************
1528 : : // Advance systems of equations
1529 : : // *****************************************************************************
1530 : : {
1531 : 6910 : auto d = Disc();
1532 : : const auto npoin = m_u.nunk();
1533 : : const auto ncomp = m_u.nprop();
1534 : : const auto& vol = d->Vol();
1535 : : const auto& lid = d->Lid();
1536 : :
1537 : : // Combine own and communicated contributions to rhs
1538 [ + + ]: 180810 : for (const auto& [g,r] : m_rhsc) {
1539 : 173900 : auto i = tk::cref_find( lid, g );
1540 [ + + ]: 801020 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1541 : : }
1542 : 6910 : tk::destroy(m_rhsc);
1543 : :
1544 [ + + ]: 6910 : if (m_stage == 0) m_un = m_u;
1545 : :
1546 : : tk::Fields u;
1547 [ + + ]: 6910 : std::size_t cstart = m_freezeflow > 1.0 ? 3 : 0;
1548 : : if (cstart) u = m_u;
1549 : :
1550 : : auto eps = std::numeric_limits<tk::real>::epsilon();
1551 [ + + ][ + - ]: 6910 : if (g_cfg.get< tag::theta >() < eps || m_stage+1 < m_rkcoef.size()) {
1552 : :
1553 : : // Apply rhs in explicit solve
1554 : 6870 : auto dt = m_rkcoef[m_stage] * d->Dt();
1555 [ + + ]: 1004980 : for (std::size_t i=0; i<npoin; ++i)
1556 [ + + ]: 4457660 : for (std::size_t c=0; c<ncomp; ++c)
1557 : 3459550 : m_u(i,c) = m_un(i,c) - dt*m_rhs(i,c)/vol[i];
1558 : :
1559 : : // Continue to advective-diffusive prediction
1560 [ + - ]: 6870 : pred();
1561 : :
1562 [ + + ]: 6870 : if (cstart) {
1563 [ + + ]: 72120 : for (std::size_t i=0; i<npoin; ++i) {
1564 [ + + ]: 288240 : for (std::size_t c=0; c<cstart; ++c) {
1565 : 216180 : m_u(i,c) = u(i,c);
1566 : : }
1567 : : }
1568 : : }
1569 : :
1570 : : } else {
1571 : :
1572 : : // Configure momentum BCs
1573 : : std::unordered_map< std::size_t,
1574 : : std::vector< std::pair< int, tk::real > > > dirbc;
1575 : :
1576 [ + - ]: 40 : if (m_np < 3) { // only during startup and first time step
1577 : : // Configure Dirichlet BCs
1578 : 40 : std::size_t nmask = ncomp + 1;
1579 : : Assert( m_dirbcmask.size() % nmask == 0, "Size mismatch" );
1580 [ + + ]: 24520 : for (std::size_t i=0; i<m_dirbcmask.size()/nmask; ++i) {
1581 [ + - ]: 24480 : auto p = m_dirbcmask[i*nmask+0]; // local node id
1582 : : auto& bc = dirbc[p];
1583 [ + - ]: 24480 : bc.resize( ncomp );
1584 [ + + ]: 97920 : for (std::size_t c=0; c<ncomp; ++c) {
1585 : 73440 : bc[c] = { m_dirbcmask[i*nmask+1+c], 0.0 };
1586 : : }
1587 : : }
1588 [ + - ][ + + ]: 8200 : for (auto p : m_noslipbcnodes) {
1589 : : auto& bc = dirbc[p];
1590 [ + - ]: 8160 : bc.resize( ncomp );
1591 [ + + ]: 32640 : for (std::size_t c=0; c<ncomp; ++c) {
1592 : : bc[c] = { 1, 0.0 };
1593 : : }
1594 : : }
1595 : : }
1596 : :
1597 : : // Initialize semi-implicit momentum/transport solve
1598 : 40 : m_cgmom[ thisIndex ].ckLocal()->
1599 [ + - ][ + - ]: 160 : init( {}, m_rhs.vec(), {}, dirbc, m_np<3, g_cfg.get< tag::mom_pc >(),
[ - + ][ - - ]
1600 [ + - ][ + - ]: 120 : CkCallback( CkIndex_ChoCG::msolve(), thisProxy[thisIndex] ) );
[ - + ][ - - ]
1601 : :
1602 : : }
1603 : 6910 : }
1604 : :
1605 : : void
1606 : 40 : ChoCG::msolve()
1607 : : // *****************************************************************************
1608 : : // Solve for momentum/transport system of equations
1609 : : // *****************************************************************************
1610 : : {
1611 : 40 : auto iter = g_cfg.get< tag::mom_iter >();
1612 : 40 : auto tol = g_cfg.get< tag::mom_tol >();
1613 : 40 : auto verbose = g_cfg.get< tag::mom_verbose >();
1614 : :
1615 [ + - ]: 80 : m_cgmom[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose,
1616 [ + - ][ + - ]: 120 : CkCallback( CkIndex_ChoCG::msolved(), thisProxy[thisIndex] ) );
[ - + ][ - - ]
1617 : 40 : }
1618 : :
1619 : : void
1620 : 40 : ChoCG::msolved()
1621 : : // *****************************************************************************
1622 : : // Continue after momentum/transport solve in semi-implcit solve
1623 : : // *****************************************************************************
1624 : : {
1625 : 40 : auto d = Disc();
1626 : : const auto npoin = m_u.nunk();
1627 : : const auto ncomp = m_u.nprop();
1628 : :
1629 [ + + ][ + - ]: 80 : if (thisIndex == 0) d->mit( m_cgmom[ thisIndex ].ckLocal()->it() );
1630 : :
1631 : : // Update momentum/transport solution in semi-implicit solve
1632 : 40 : auto& du = m_cgmom[ thisIndex ].ckLocal()->solution();
1633 [ + + ]: 24520 : for (std::size_t i=0; i<npoin; ++i) {
1634 [ + + ]: 97920 : for (std::size_t c=0; c<ncomp; ++c) {
1635 : 73440 : m_u(i,c) = m_un(i,c) + du[i*ncomp+c];
1636 : : }
1637 : : }
1638 : :
1639 : : // Continue to advective-diffusive prediction
1640 : 40 : pred();
1641 : 40 : }
1642 : :
1643 : : void
1644 : 6910 : ChoCG::pred()
1645 : : // *****************************************************************************
1646 : : // Compute advective-diffusive prediction of momentum/transport
1647 : : // *****************************************************************************
1648 : : {
1649 : 6910 : auto d = Disc();
1650 : :
1651 : : // Configure and apply scalar source to solution (if defined)
1652 : 6910 : auto src = problems::PHYS_SRC();
1653 [ - + ][ - - ]: 6910 : if (src) src( d->Coord(), d->T(), m_u );
1654 : :
1655 : : // Enforce boundary conditions
1656 [ + - ]: 6910 : BC( m_u, d->T() + m_rkcoef[m_stage] * d->Dt() );
1657 : :
1658 : : // Compute velocity gradients if needed
1659 [ + + ]: 6910 : if (g_cfg.get< tag::flux >() == "damp4") {
1660 [ + - ][ + - ]: 3320 : thisProxy[ thisIndex ].wait4vgrad();
[ - - ]
1661 [ + - ]: 3320 : velgrad();
1662 : : } else {
1663 [ + - ]: 3590 : corr();
1664 : : }
1665 : 6910 : }
1666 : :
1667 : : void
1668 : 6910 : ChoCG::corr()
1669 : : // *****************************************************************************
1670 : : // Compute pressure correction
1671 : : // *****************************************************************************
1672 : : {
1673 : : // Finalize computing velocity gradients
1674 [ + + ]: 6910 : if (g_cfg.get< tag::flux >() == "damp4") fingrad( m_vgrad, m_vgradc );
1675 : :
1676 [ + + ]: 6910 : if (++m_stage < m_rkcoef.size()) {
1677 : :
1678 : : // Activate SDAG wait for next time step stage
1679 [ + - ]: 2260 : thisProxy[ thisIndex ].wait4rhs();
1680 : : // Continue to next time stage of momentum/transport prediction
1681 [ + - ]: 4520 : contribute( CkCallback(CkReductionTarget(ChoCG,rhs), thisProxy) );
1682 : :
1683 : : } else {
1684 : :
1685 : : // Reset Runge-Kutta stage counter
1686 : 4650 : m_stage = 0;
1687 : : // Continue to pressure correction and projection
1688 : 4650 : div( m_u );
1689 : :
1690 : : }
1691 : 6910 : }
1692 : :
1693 : : void
1694 : 4682 : ChoCG::diag()
1695 : : // *****************************************************************************
1696 : : // Compute diagnostics
1697 : : // *****************************************************************************
1698 : : {
1699 : 4682 : auto d = Disc();
1700 : :
1701 : : // Increase number of iterations and physical time
1702 : 4682 : d->next();
1703 : :
1704 : : // Compute diagnostics, e.g., residuals
1705 : 4682 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1706 : 4682 : const auto& dp = m_cgpre[ thisIndex ].ckLocal()->solution();
1707 : 4682 : auto diagnostics = m_diag.precompute( *d, m_u, m_un, m_pr, dp, diag_iter );
1708 : :
1709 : : // Evaluate residuals
1710 [ - + ][ - - ]: 4682 : if (!diagnostics) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
1711 : 4682 : }
1712 : :
1713 : : void
1714 : 4682 : ChoCG::evalres( const std::vector< tk::real >& )
1715 : : // *****************************************************************************
1716 : : // Evaluate residuals
1717 : : // *****************************************************************************
1718 : : {
1719 : 4682 : refine();
1720 : 4682 : }
1721 : :
1722 : : void
1723 : 4682 : ChoCG::refine()
1724 : : // *****************************************************************************
1725 : : // Optionally refine/derefine mesh
1726 : : // *****************************************************************************
1727 : : {
1728 : 4682 : auto d = Disc();
1729 : :
1730 : : // See if this is the last time step
1731 [ + + ]: 4682 : if (d->finished()) m_finished = 1;
1732 : :
1733 : 4682 : auto dtref = g_cfg.get< tag::href_dt >();
1734 : 4682 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
1735 : :
1736 : : // if t>0 refinement enabled and we hit the frequency
1737 [ - + ][ - - ]: 4682 : if (dtref && !(d->It() % dtfreq)) { // refine
1738 : :
1739 : 0 : d->refined() = 1;
1740 : 0 : d->startvol();
1741 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1742 : :
1743 : : // Activate SDAG waits for re-computing the integrals
1744 [ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1745 : :
1746 : : } else { // do not refine
1747 : :
1748 : 4682 : d->refined() = 0;
1749 : 4682 : feop_complete();
1750 : 4682 : resize_complete();
1751 : :
1752 : : }
1753 : 4682 : }
1754 : :
1755 : : void
1756 : 0 : ChoCG::resizePostAMR(
1757 : : const std::vector< std::size_t >& /*ginpoel*/,
1758 : : const tk::UnsMesh::Chunk& chunk,
1759 : : const tk::UnsMesh::Coords& coord,
1760 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1761 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1762 : : const std::set< std::size_t >& removedNodes,
1763 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1764 : : nodeCommMap,
1765 : : const std::map< int, std::vector< std::size_t > >& bface,
1766 : : const std::map< int, std::vector< std::size_t > >& bnode,
1767 : : const std::vector< std::size_t >& triinpoel )
1768 : : // *****************************************************************************
1769 : : // Receive new mesh from Refiner
1770 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1771 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1772 : : //! \param[in] coord New mesh node coordinates
1773 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1774 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1775 : : //! \param[in] removedNodes Newly removed mesh node local ids
1776 : : //! \param[in] nodeCommMap New node communication map
1777 : : //! \param[in] bface Boundary-faces mapped to side set ids
1778 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1779 : : //! \param[in] triinpoel Boundary-face connectivity
1780 : : // *****************************************************************************
1781 : : {
1782 : 0 : auto d = Disc();
1783 : :
1784 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1785 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1786 : :
1787 : : // Resize mesh data structures after mesh refinement
1788 : 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1789 : :
1790 : : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
1791 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1792 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1793 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1794 : :
1795 : : // Remove newly removed nodes from solution vectors
1796 : 0 : m_u.rm( removedNodes );
1797 : : //m_pr.rm( removedNodes );
1798 : 0 : m_rhs.rm( removedNodes );
1799 : :
1800 : : // Resize auxiliary solution vectors
1801 : : auto npoin = coord[0].size();
1802 : : m_u.resize( npoin );
1803 : 0 : m_pr.resize( npoin );
1804 : : m_rhs.resize( npoin );
1805 : :
1806 : : // Update solution on new mesh
1807 [ - - ]: 0 : for (const auto& n : addedNodes)
1808 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1809 : : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
1810 : : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
1811 : : "Indices of parent-edge nodes out of bounds post-AMR");
1812 : 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
1813 : : }
1814 : :
1815 : : // Update physical-boundary node-, face-, and element lists
1816 : : m_bnode = bnode;
1817 : : m_bface = bface;
1818 : 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1819 : :
1820 : 0 : auto meshid = d->MeshId();
1821 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1822 : 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1823 : 0 : }
1824 : :
1825 : : void
1826 : 974 : ChoCG::writeFields( CkCallback cb )
1827 : : // *****************************************************************************
1828 : : // Output mesh-based fields to file
1829 : : //! \param[in] cb Function to continue with after the write
1830 : : // *****************************************************************************
1831 : : {
1832 [ + + ]: 974 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1833 : :
1834 : 390 : auto d = Disc();
1835 : :
1836 : : // Field output
1837 : :
1838 : : std::vector< std::string > nodefieldnames{
1839 [ - + ][ + + ]: 2340 : "velocityx", "velocityy", "velocityz", "divergence", "pressure" };
[ - + ][ - - ]
[ - - ]
1840 : :
1841 : : std::vector< std::vector< tk::real > > nodefields;
1842 : :
1843 [ + - ]: 390 : nodefields.push_back( m_u.extract(0) );
1844 [ + - ]: 390 : nodefields.push_back( m_u.extract(1) );
1845 [ + - ]: 390 : nodefields.push_back( m_u.extract(2) );
1846 : :
1847 : : // Divide weak result by nodal volume
1848 : : const auto& vol = d->Vol();
1849 [ + + ]: 91711 : for (std::size_t i=0; i<m_div.size(); ++i) m_div[i] /= vol[i];
1850 [ + - ]: 390 : nodefields.push_back( m_div );
1851 : :
1852 [ + - ]: 390 : nodefields.push_back( m_pr ) ;
1853 : :
1854 : : //nodefieldnames.push_back( "dp/dx" );
1855 : : //nodefieldnames.push_back( "dp/dy" );
1856 : : //nodefieldnames.push_back( "dp/dz" );
1857 : : //nodefields.push_back( m_pgrad.extract(0) );
1858 : : //nodefields.push_back( m_pgrad.extract(1) );
1859 : : //nodefields.push_back( m_pgrad.extract(2) );
1860 : :
1861 : : //nodefieldnames.push_back( "fx" );
1862 : : //nodefieldnames.push_back( "fy" );
1863 : : //nodefieldnames.push_back( "fz" );
1864 : : //nodefields.push_back( m_flux.extract(0) );
1865 : : //nodefields.push_back( m_flux.extract(1) );
1866 : : //nodefields.push_back( m_flux.extract(2) );
1867 : :
1868 : : //nodefieldnames.push_back( "du/dx" );
1869 : : //nodefieldnames.push_back( "du/dy" );
1870 : : //nodefieldnames.push_back( "du/dz" );
1871 : : //nodefieldnames.push_back( "dv/dx" );
1872 : : //nodefieldnames.push_back( "dv/dy" );
1873 : : //nodefieldnames.push_back( "dv/dz" );
1874 : : //nodefieldnames.push_back( "dw/dx" );
1875 : : //nodefieldnames.push_back( "dw/dy" );
1876 : : //nodefieldnames.push_back( "dw/dz" );
1877 : : //nodefields.push_back( m_vgrad.extract(0) );
1878 : : //nodefields.push_back( m_vgrad.extract(1) );
1879 : : //nodefields.push_back( m_vgrad.extract(2) );
1880 : : //nodefields.push_back( m_vgrad.extract(3) );
1881 : : //nodefields.push_back( m_vgrad.extract(4) );
1882 : : //nodefields.push_back( m_vgrad.extract(5) );
1883 : : //nodefields.push_back( m_vgrad.extract(6) );
1884 : : //nodefields.push_back( m_vgrad.extract(7) );
1885 : : //nodefields.push_back( m_vgrad.extract(8) );
1886 : :
1887 : : auto ncomp = m_u.nprop();
1888 [ + + ]: 492 : for (std::size_t c=0; c<ncomp-3; ++c) {
1889 [ + - ]: 102 : nodefieldnames.push_back( "c" + std::to_string(c) );
1890 [ + - ]: 204 : nodefields.push_back( m_u.extract(3+c) );
1891 : : }
1892 : :
1893 : : // query function to evaluate analytic solution (if defined)
1894 [ + - ]: 390 : auto sol = problems::SOL();
1895 : :
1896 [ + + ]: 390 : if (sol) {
1897 : : const auto& coord = d->Coord();
1898 : : const auto& x = coord[0];
1899 : : const auto& y = coord[1];
1900 : : const auto& z = coord[2];
1901 : : auto an = m_u;
1902 [ + + ]: 18334 : for (std::size_t i=0; i<an.nunk(); ++i) {
1903 [ - + ]: 18168 : auto s = sol( x[i], y[i], z[i], d->T() );
1904 [ + + ]: 86578 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
1905 : : }
1906 [ + - ]: 166 : nodefieldnames.push_back( "velocity_analyticx" );
1907 [ + - ]: 166 : nodefields.push_back( an.extract(0) );
1908 [ + - ]: 166 : nodefieldnames.push_back( "velocity_analyticy" );
1909 [ + - ]: 166 : nodefields.push_back( an.extract(1) );
1910 [ + - ]: 166 : nodefieldnames.push_back( "velocity_analyticz" );
1911 [ + - ]: 166 : nodefields.push_back( an.extract(2) );
1912 [ + + ]: 268 : for (std::size_t c=0; c<ncomp-3; ++c) {
1913 [ + - ]: 102 : nodefieldnames.push_back( nodefieldnames[5+c] + "_analytic" );
1914 [ + - ][ - - ]: 204 : nodefields.push_back( an.extract(3+c) );
1915 : : }
1916 : : }
1917 : :
1918 : : // also output analytic pressure (if defined)
1919 [ + - ]: 390 : auto pressure_sol = problems::PRESSURE_SOL();
1920 [ + + ]: 390 : if (pressure_sol) {
1921 : : const auto& coord = d->Coord();
1922 : : const auto& x = coord[0];
1923 : : const auto& y = coord[1];
1924 : : const auto& z = coord[2];
1925 [ + - ]: 64 : auto ap = m_pr;
1926 [ + + ]: 4326 : for (std::size_t i=0; i<ap.size(); ++i) {
1927 [ - + ]: 8524 : ap[i] = pressure_sol( x[i], y[i], z[i] );
1928 : : }
1929 [ + - ][ - - ]: 64 : nodefieldnames.push_back( "analytic" );
1930 [ + - ]: 64 : nodefields.push_back( ap );
1931 : : }
1932 : :
1933 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
1934 : :
1935 : : // Surface output
1936 : :
1937 : : std::vector< std::string > nodesurfnames;
1938 : : std::vector< std::vector< tk::real > > nodesurfs;
1939 : :
1940 : : const auto& f = g_cfg.get< tag::fieldout >();
1941 : :
1942 [ + + ]: 390 : if (!f.empty()) {
1943 : : std::size_t nc = 5;
1944 [ + - ]: 3 : nodesurfnames.push_back( "velocityx" );
1945 [ + - ]: 3 : nodesurfnames.push_back( "velocityy" );
1946 [ + - ]: 3 : nodesurfnames.push_back( "velocityz" );
1947 [ + - ]: 3 : nodesurfnames.push_back( "divergence" );
1948 [ + - ]: 3 : nodesurfnames.push_back( "pressure" );
1949 : :
1950 [ + - ]: 3 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
1951 [ + - ]: 3 : std::set< int > outsets( begin(f), end(f) );
1952 [ + + ]: 6 : for (auto sideset : outsets) {
1953 : : auto b = bnode.find(sideset);
1954 [ - + ]: 3 : if (b == end(bnode)) continue;
1955 : : const auto& nodes = b->second;
1956 : : auto i = nodesurfs.size();
1957 : : nodesurfs.insert( end(nodesurfs), nc,
1958 [ + - ]: 6 : std::vector< tk::real >( nodes.size() ) );
1959 : : std::size_t j = 0;
1960 [ + + ]: 1278 : for (auto n : nodes) {
1961 [ + - ]: 1275 : const auto s = m_u[n];
1962 : : std::size_t p = 0;
1963 : 1275 : nodesurfs[i+(p++)][j] = s[0];
1964 : 1275 : nodesurfs[i+(p++)][j] = s[1];
1965 : 1275 : nodesurfs[i+(p++)][j] = s[2];
1966 : 1275 : nodesurfs[i+(p++)][j] = m_div[n];
1967 : 1275 : nodesurfs[i+(p++)][j] = m_pr[n];
1968 : 1275 : ++j;
1969 : : }
1970 : : }
1971 : : }
1972 : :
1973 : : // Send mesh and fields data (solution dump) for output to file
1974 [ + - ][ + - ]: 780 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
1975 [ + - ]: 390 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
1976 : : {}, nodefields, {}, nodesurfs, cb );
1977 [ + - ][ + - ]: 1560 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
1978 : :
1979 : : void
1980 : 4682 : ChoCG::out()
1981 : : // *****************************************************************************
1982 : : // Output mesh field data
1983 : : // *****************************************************************************
1984 : : {
1985 : 4682 : auto d = Disc();
1986 : :
1987 : : // Time history
1988 [ + - ][ + - ]: 4682 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
1989 : : auto ncomp = m_u.nprop();
1990 : : const auto& inpoel = d->Inpoel();
1991 : 0 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
1992 : : std::size_t j = 0;
1993 [ - - ]: 0 : for (const auto& p : d->Hist()) {
1994 [ - - ]: 0 : auto e = p.get< tag::elem >(); // host element id
1995 : : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
1996 [ - - ]: 0 : hist[j].resize( ncomp+1, 0.0 );
1997 [ - - ]: 0 : for (std::size_t i=0; i<4; ++i) {
1998 [ - - ]: 0 : const auto u = m_u[ inpoel[e*4+i] ];
1999 : 0 : hist[j][0] += n[i] * u[0];
2000 : 0 : hist[j][1] += n[i] * u[1]/u[0];
2001 : 0 : hist[j][2] += n[i] * u[2]/u[0];
2002 : 0 : hist[j][3] += n[i] * u[3]/u[0];
2003 : 0 : hist[j][4] += n[i] * u[4]/u[0];
2004 : 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];
2005 : 0 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
2006 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
2007 : : }
2008 : 0 : ++j;
2009 : : }
2010 [ - - ]: 0 : d->history( std::move(hist) );
2011 : 0 : }
2012 : :
2013 : : // Field data
2014 [ + + ][ + - ]: 4682 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
2015 [ + - ][ + - ]: 1674 : writeFields( CkCallback(CkIndex_ChoCG::integrals(), thisProxy[thisIndex]) );
2016 : : } else {
2017 : 4124 : integrals();
2018 : : }
2019 : 4682 : }
2020 : :
2021 : : void
2022 : 4682 : ChoCG::integrals()
2023 : : // *****************************************************************************
2024 : : // Compute integral quantities for output
2025 : : // *****************************************************************************
2026 : : {
2027 : 4682 : auto d = Disc();
2028 : :
2029 [ + + ][ + - ]: 4682 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ]
2030 : :
2031 : : using namespace integrals;
2032 [ + - ]: 20 : std::vector< std::map< int, tk::real > > ints( NUMINT );
2033 : : // Prepend integral vector with metadata on the current time step:
2034 : : // current iteration count, current physical time, time step size
2035 [ + - ][ + - ]: 20 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
2036 [ + - ][ + - ]: 20 : ints[ TIME ][ 0 ] = d->T();
2037 [ + - ][ + - ]: 20 : ints[ DT ][ 0 ] = d->Dt();
2038 : :
2039 : : // Compute integrals requested for surfaces requested
2040 : : const auto& reqv = g_cfg.get< tag::integout_integrals >();
2041 : 20 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
2042 [ + - ][ - + ]: 40 : if (req.count("mass_flow_rate")) {
2043 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
2044 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
2045 : : const auto& nodes = sint.first;
2046 : : const auto& ndA = sint.second;
2047 : : auto n = ndA.data();
2048 [ - - ]: 0 : for (auto p : nodes) {
2049 : 0 : mfr += n[0]*m_u(p,0) + n[1]*m_u(p,1) + n[2]*m_u(p,2);
2050 : 0 : n += 3;
2051 : : }
2052 : : }
2053 : : }
2054 [ + - ][ + - ]: 40 : if (req.count("force")) {
2055 : 20 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
2056 [ + + ]: 40 : for (const auto& [s,sint] : m_surfint) {
2057 [ + - ]: 20 : auto& fx = ints[ FORCE_X ][ s ];
2058 [ + - ]: 20 : auto& fy = ints[ FORCE_Y ][ s ];
2059 [ + - ]: 20 : auto& fz = ints[ FORCE_Z ][ s ];
2060 : : const auto& nodes = sint.first;
2061 : : const auto& ndA = sint.second;
2062 : : auto n = ndA.data();
2063 [ + + ]: 8520 : for (auto p : nodes) {
2064 : : // pressure force
2065 : 8500 : fx -= n[0]*m_pr[p];
2066 : 8500 : fy -= n[1]*m_pr[p];
2067 : 8500 : fz -= n[2]*m_pr[p];
2068 : : // viscous force
2069 : 8500 : fx += mu*(m_vgrad(p,0)*n[0] + m_vgrad(p,1)*n[1] + m_vgrad(p,2)*n[2]);
2070 : 8500 : fy += mu*(m_vgrad(p,3)*n[0] + m_vgrad(p,4)*n[1] + m_vgrad(p,5)*n[2]);
2071 : 8500 : fz += mu*(m_vgrad(p,6)*n[0] + m_vgrad(p,7)*n[1] + m_vgrad(p,8)*n[2]);
2072 : 8500 : n += 3;
2073 : : }
2074 : : }
2075 : : }
2076 : :
2077 [ + - ]: 20 : auto stream = serialize( d->MeshId(), ints );
2078 [ + - ][ + - ]: 40 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
2079 [ + - ][ - - ]: 20 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
2080 : :
2081 : 20 : } else {
2082 : :
2083 : 4662 : step();
2084 : :
2085 : : }
2086 : 4682 : }
2087 : :
2088 : : void
2089 : 4266 : ChoCG::evalLB( int nrestart )
2090 : : // *****************************************************************************
2091 : : // Evaluate whether to do load balancing
2092 : : //! \param[in] nrestart Number of times restarted
2093 : : // *****************************************************************************
2094 : : {
2095 : 4266 : auto d = Disc();
2096 : :
2097 : : // Detect if just returned from a checkpoint and if so, zero timers and
2098 : : // finished flag
2099 [ + + ]: 4266 : if (d->restarted( nrestart )) m_finished = 0;
2100 : :
2101 : 4266 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
2102 [ + + ]: 4266 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
2103 : :
2104 : : // Load balancing if user frequency is reached or after the second time-step
2105 [ + + ][ + + ]: 4266 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
2106 : :
2107 : 3446 : AtSync();
2108 [ - + ]: 3446 : if (nonblocking) dt();
2109 : :
2110 : : } else {
2111 : :
2112 : 820 : dt();
2113 : :
2114 : : }
2115 : 4266 : }
2116 : :
2117 : : void
2118 : 4261 : ChoCG::evalRestart()
2119 : : // *****************************************************************************
2120 : : // Evaluate whether to save checkpoint/restart
2121 : : // *****************************************************************************
2122 : : {
2123 : 4261 : auto d = Disc();
2124 : :
2125 : 4261 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
2126 : 4261 : const auto benchmark = g_cfg.get< tag::benchmark >();
2127 : :
2128 [ + + ][ - + ]: 4261 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
2129 : :
2130 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
2131 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
2132 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
2133 : :
2134 : : } else {
2135 : :
2136 : 4261 : evalLB( /* nrestart = */ -1 );
2137 : :
2138 : : }
2139 : 4261 : }
2140 : :
2141 : : void
2142 : 4682 : ChoCG::step()
2143 : : // *****************************************************************************
2144 : : // Evaluate whether to continue with next time step
2145 : : // *****************************************************************************
2146 : : {
2147 : 4682 : auto d = Disc();
2148 : :
2149 : : // Output one-liner status report to screen
2150 [ + + ]: 4682 : if(thisIndex == 0) d->status();
2151 : :
2152 [ + + ]: 4682 : if (not m_finished) {
2153 : :
2154 : 4261 : evalRestart();
2155 : :
2156 : : } else {
2157 : :
2158 : 421 : auto meshid = d->MeshId();
2159 [ + - ]: 842 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2160 : 421 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
2161 : :
2162 : : }
2163 : 4682 : }
2164 : :
2165 : : #include "NoWarning/chocg.def.h"
|
