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