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