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