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 [ + + ]: 23682 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
104 [ + - ][ + + ]: 23297 : if (!map.count(p)) map[p] = n++;
[ + - ]
105 [ + + ]: 220897 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
106 [ + - ][ + + ]: 197600 : 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_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 [ + - ]: 135 : auto& n = dir[ k->first ];
223 [ + + ]: 21023 : 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 [ + - ]: 137 : auto& n = dir[ k->first ];
237 [ + + ]: 22062 : for (auto g : k->second) {
238 [ + - ][ + - ]: 21925 : 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 [ + - ]: 137 : std::vector< double > v( ncomp, 0.0 );
264 [ + - ]: 137 : auto m = dirval.find( vec[0] );
265 [ + + ]: 137 : if (m != end(dirval)) {
266 [ - + ][ - - ]: 15 : ErrChk( m->second.size() == ncomp, "Incorrect Dirichlet BC val ncomp" );
[ - - ][ - - ]
267 [ + - ]: 15 : v = m->second;
268 : : }
269 [ + + ]: 13812 : for (auto p : n->second) {
270 [ + - ]: 13675 : auto& mv = dirbcset[p]; // mask & value
271 [ + + ][ + - ]: 13675 : if (mv.second.empty()) mv.second = v;
272 : 13675 : auto& mval = mv.first;
273 [ + + ][ + - ]: 13675 : if (mval.empty()) mval.resize( ncomp, 0 );
274 [ + + ]: 74186 : for (std::size_t c=0; c<ncomp; ++c) {
275 [ + + ]: 60511 : 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 : 137 : }
281 : : }
282 : :
283 : : // Compile streamable list of nodes + Dirichlet BC components mask and values
284 : 770 : tk::destroy( mask );
285 [ + + ]: 13657 : for (const auto& [p,mv] : dirbcset) {
286 [ + - ]: 12887 : mask.push_back( p );
287 [ + - ]: 12887 : mask.insert( end(mask), begin(mv.first), end(mv.first) );
288 [ + - ]: 12887 : val.push_back( static_cast< double >( p ) );
289 [ + - ]: 12887 : 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 : : // Prepare Dirichlet boundary conditions data structures
306 : 385 : setupDirBC( g_cfg.get< tag::bc_dir >(), g_cfg.get< tag::bc_dirval >(),
307 : 385 : m_u.nprop(), m_dirbcmask, m_dirbcval );
308 : 385 : setupDirBC( g_cfg.get< tag::pre_bc_dir >(), g_cfg.get< tag::pre_bc_dirval >(),
309 : 385 : 1, m_dirbcmaskp, m_dirbcvalp );
310 : :
311 : : // Compute local contributions to boundary normals and integrals
312 : 385 : bndint();
313 : : // Compute local contributions to domain edge integrals
314 : 385 : domint();
315 : :
316 : : // Send boundary point normal contributions to neighbor chares
317 [ + + ]: 385 : if (d->NodeCommMap().empty()) {
318 : 9 : comnorm_complete();
319 : : } else {
320 [ + + ]: 4230 : for (const auto& [c,nodes] : d->NodeCommMap()) {
321 : 3854 : decltype(m_bnorm) exp;
322 [ + + ]: 21076 : for (auto i : nodes) {
323 [ + + ]: 50198 : for (const auto& [s,b] : m_bnorm) {
324 [ + - ]: 32976 : auto k = b.find(i);
325 [ + + ][ + - ]: 32976 : if (k != end(b)) exp[s][i] = k->second;
[ + - ]
326 : : }
327 : : }
328 [ + - ][ + - ]: 3854 : thisProxy[c].comnorm( exp );
329 : 3854 : }
330 : : }
331 : 385 : ownnorm_complete();
332 : 385 : }
333 : :
334 : : void
335 : 385 : LohCG::bndint()
336 : : // *****************************************************************************
337 : : // Compute local contributions to boundary normals and integrals
338 : : // *****************************************************************************
339 : : {
340 [ + - ]: 385 : auto d = Disc();
341 : 385 : const auto& coord = d->Coord();
342 : 385 : const auto& gid = d->Gid();
343 : 385 : const auto& x = coord[0];
344 : 385 : const auto& y = coord[1];
345 : 385 : const auto& z = coord[2];
346 : :
347 : : // Lambda to compute the inverse distance squared between boundary face
348 : : // centroid and boundary point. Here p is the global node id and c is the
349 : : // the boundary face centroid.
350 : 99456 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
351 : 99456 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
352 : 99456 : (c[1] - y[p]) * (c[1] - y[p]) +
353 : 99456 : (c[2] - z[p]) * (c[2] - z[p]) );
354 : 385 : };
355 : :
356 : 385 : tk::destroy( m_bnorm );
357 : 385 : tk::destroy( m_bndpoinint );
358 : :
359 [ + + ]: 1216 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
360 [ + + ]: 33983 : for (auto f : faceids) { // for all side set triangles
361 : 33152 : const auto N = m_triinpoel.data() + f*3;
362 : : const std::array< tk::real, 3 >
363 : 33152 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
364 : 33152 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
365 : 33152 : auto n = tk::cross( ba, ca );
366 : 33152 : auto A2 = tk::length( n );
367 : 33152 : n[0] /= A2;
368 : 33152 : n[1] /= A2;
369 : 33152 : n[2] /= A2;
370 : : const tk::real centroid[3] = {
371 : 33152 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
372 : 33152 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
373 : 66304 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
374 [ + + ]: 132608 : for (const auto& [i,j] : tk::lpoet) {
375 : 99456 : auto p = N[i];
376 : 99456 : tk::real r = invdistsq( centroid, p );
377 [ + - ]: 99456 : auto& v = m_bnorm[setid]; // associate side set id
378 [ + - ]: 99456 : auto& bpn = v[gid[p]]; // associate global node id of bnd pnt
379 : 99456 : bpn[0] += r * n[0]; // inv.dist.sq-weighted normal
380 : 99456 : bpn[1] += r * n[1];
381 : 99456 : bpn[2] += r * n[2];
382 : 99456 : bpn[3] += r; // inv.dist.sq of node from centroid
383 [ + - ]: 99456 : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
384 : 99456 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
385 : 99456 : b[1] += n[1] * A2 / 6.0;
386 : 99456 : b[2] += n[2] * A2 / 6.0;
387 : : }
388 : : }
389 : : }
390 : 385 : }
391 : :
392 : : void
393 : 385 : LohCG::domint()
394 : : // *****************************************************************************
395 : : //! Compute local contributions to domain edge integrals
396 : : // *****************************************************************************
397 : : {
398 : 385 : auto d = Disc();
399 : :
400 : 385 : const auto& gid = d->Gid();
401 : 385 : const auto& inpoel = d->Inpoel();
402 : :
403 : 385 : const auto& coord = d->Coord();
404 : 385 : const auto& x = coord[0];
405 : 385 : const auto& y = coord[1];
406 : 385 : const auto& z = coord[2];
407 : :
408 : 385 : tk::destroy( m_domedgeint );
409 : :
410 [ + + ]: 54501 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
411 : 54116 : const auto N = inpoel.data() + e*4;
412 : : const std::array< tk::real, 3 >
413 : 54116 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
414 : 54116 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
415 : 54116 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
416 : 54116 : const auto J = tk::triple( ba, ca, da ); // J = 6V
417 [ - + ][ - - ]: 54116 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
418 : : std::array< std::array< tk::real, 3 >, 4 > grad;
419 : 54116 : grad[1] = tk::cross( ca, da );
420 : 54116 : grad[2] = tk::cross( da, ba );
421 : 54116 : grad[3] = tk::cross( ba, ca );
422 [ + + ]: 216464 : for (std::size_t i=0; i<3; ++i)
423 : 162348 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
424 [ + + ]: 378812 : for (const auto& [p,q] : tk::lpoed) {
425 : 324696 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
426 : 324696 : tk::real sig = 1.0;
427 [ + + ]: 324696 : if (ed[0] > ed[1]) {
428 : 122178 : std::swap( ed[0], ed[1] );
429 : 122178 : sig = -1.0;
430 : : }
431 [ + - ]: 324696 : auto& n = m_domedgeint[ ed ];
432 : 324696 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
433 : 324696 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
434 : 324696 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
435 : 324696 : n[3] += tk::dot( grad[p], grad[q] ) / J / 6.0;
436 : : }
437 : : }
438 : 385 : }
439 : :
440 : : void
441 : 3854 : LohCG::comnorm( const decltype(m_bnorm)& inbnd )
442 : : // *****************************************************************************
443 : : // Receive contributions to boundary point normals on chare-boundaries
444 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
445 : : // *****************************************************************************
446 : : {
447 : : // Buffer up incoming boundary point normal vector contributions
448 [ + + ]: 6817 : for (const auto& [s,b] : inbnd) {
449 [ + - ]: 2963 : auto& bndnorm = m_bnormc[s];
450 [ + + ]: 12570 : for (const auto& [p,n] : b) {
451 [ + - ]: 9607 : auto& norm = bndnorm[p];
452 : 9607 : norm[0] += n[0];
453 : 9607 : norm[1] += n[1];
454 : 9607 : norm[2] += n[2];
455 : 9607 : norm[3] += n[3];
456 : : }
457 : : }
458 : :
459 [ + + ]: 3854 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
460 : 376 : m_nnorm = 0;
461 : 376 : comnorm_complete();
462 : : }
463 : 3854 : }
464 : :
465 : : void
466 : 802 : LohCG::registerReducers()
467 : : // *****************************************************************************
468 : : // Configure Charm++ reduction types initiated from this chare array
469 : : //! \details Since this is a [initnode] routine, the runtime system executes the
470 : : //! routine exactly once on every logical node early on in the Charm++ init
471 : : //! sequence. Must be static as it is called without an object. See also:
472 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
473 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
474 : : // *****************************************************************************
475 : : {
476 : 802 : NodeDiagnostics::registerReducers();
477 : 802 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
478 : 802 : }
479 : :
480 : : void
481 : : // cppcheck-suppress unusedFunction
482 : 3525 : LohCG::ResumeFromSync()
483 : : // *****************************************************************************
484 : : // Return from migration
485 : : //! \details This is called when load balancing (LB) completes. The presence of
486 : : //! this function does not affect whether or not we block on LB.
487 : : // *****************************************************************************
488 : : {
489 [ - + ][ - - ]: 3525 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
490 : :
491 [ + - ]: 3525 : if (!g_cfg.get< tag::nonblocking >()) dt();
492 : 3525 : }
493 : :
494 : : void
495 : 385 : LohCG::setup( tk::real v )
496 : : // *****************************************************************************
497 : : // Start setup for solution
498 : : //! \param[in] v Total volume within user-specified box
499 : : // *****************************************************************************
500 : : {
501 : 385 : auto d = Disc();
502 : :
503 : : // Store user-defined box IC volume
504 : 385 : Disc()->Boxvol() = v;
505 : :
506 : : // Set initial conditions
507 : 385 : problems::initialize( d->Coord(), m_u, d->T(), d->BoxNodes() );
508 : :
509 : : // Query time history field output labels from all PDEs integrated
510 [ - + ]: 385 : if (!g_cfg.get< tag::histout >().empty()) {
511 : : std::vector< std::string > var
512 [ - - ][ - - ]: 0 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ]
513 : 0 : auto ncomp = m_u.nprop();
514 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c)
515 [ - - ][ - - ]: 0 : var.push_back( "c" + std::to_string(c-5) );
[ - - ]
516 [ - - ]: 0 : d->histheader( std::move(var) );
517 : 0 : }
518 : :
519 : : // Compute finite element operators
520 : 385 : feop();
521 [ - - ][ - - ]: 385 : }
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
522 : :
523 : : void
524 : 385 : LohCG::bnorm()
525 : : // *****************************************************************************
526 : : // Combine own and communicated portions of the boundary point normals
527 : : // *****************************************************************************
528 : : {
529 [ + - ]: 385 : const auto& lid = Disc()->Lid();
530 : :
531 : : // Combine own and communicated contributions to boundary point normals
532 [ + + ]: 1215 : for (const auto& [s,b] : m_bnormc) {
533 [ + - ]: 830 : auto& bndnorm = m_bnorm[s];
534 [ + + ]: 7937 : for (const auto& [g,n] : b) {
535 [ + - ]: 7107 : auto& norm = bndnorm[g];
536 : 7107 : norm[0] += n[0];
537 : 7107 : norm[1] += n[1];
538 : 7107 : norm[2] += n[2];
539 : 7107 : norm[3] += n[3];
540 : : }
541 : : }
542 : 385 : tk::destroy( m_bnormc );
543 : :
544 : : // Divide summed point normals by the sum of the inverse distance squared
545 [ + + ]: 1272 : for (auto& [s,b] : m_bnorm) {
546 [ + + ]: 25795 : for (auto& [g,n] : b) {
547 : 24908 : n[0] /= n[3];
548 : 24908 : n[1] /= n[3];
549 : 24908 : n[2] /= n[3];
550 [ - + ][ - - ]: 24908 : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
[ - - ][ - - ]
551 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
552 : : "Non-unit normal" );
553 : : }
554 : : }
555 : :
556 : : // Replace global->local ids associated to boundary point normals
557 : 385 : decltype(m_bnorm) loc;
558 [ + + ]: 1272 : for (auto& [s,b] : m_bnorm) {
559 [ + - ]: 887 : auto& bnd = loc[s];
560 [ + + ]: 25795 : for (auto&& [g,n] : b) {
561 [ + - ][ + - ]: 24908 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
562 : : }
563 : : }
564 : 385 : m_bnorm = std::move(loc);
565 : 385 : }
566 : :
567 : : void
568 : 385 : LohCG::streamable()
569 : : // *****************************************************************************
570 : : // Convert integrals into streamable data structures
571 : : // *****************************************************************************
572 : : {
573 : : // Query surface integral output nodes
574 : 385 : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
575 : 385 : const auto& is = g_cfg.get< tag::integout >();
576 [ + - ]: 385 : std::set< int > outsets( begin(is), end(is) );
577 [ + + ]: 396 : for (auto s : outsets) {
578 [ + - ]: 11 : auto m = m_bface.find(s);
579 [ + + ]: 11 : if (m != end(m_bface)) {
580 [ + - ]: 8 : auto& n = surfintnodes[ m->first ]; // associate set id
581 [ + + ]: 1700 : for (auto f : m->second) { // face ids on side set
582 : 1692 : auto t = m_triinpoel.data() + f*3;
583 [ + - ]: 1692 : n.push_back( t[0] ); // nodes on side set
584 [ + - ]: 1692 : n.push_back( t[1] );
585 [ + - ]: 1692 : n.push_back( t[2] );
586 : : }
587 : : }
588 : : }
589 [ + - ][ + + ]: 393 : for (auto& [s,n] : surfintnodes) tk::unique( n );
590 : : // Prepare surface integral data
591 : 385 : tk::destroy( m_surfint );
592 [ + - ]: 385 : const auto& gid = Disc()->Gid();
593 [ + + ]: 393 : for (auto&& [s,n] : surfintnodes) {
594 [ + - ]: 8 : auto& sint = m_surfint[s]; // associate set id
595 : 8 : auto& nodes = sint.first;
596 : 8 : auto& ndA = sint.second;
597 : 8 : nodes = std::move(n);
598 [ + - ]: 8 : ndA.resize( nodes.size()*3 );
599 : 8 : auto a = ndA.data();
600 [ + + ]: 966 : for (auto p : nodes) {
601 [ + - ]: 958 : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
602 : 958 : a[0] = b[0]; // store ni * dA
603 : 958 : a[1] = b[1];
604 : 958 : a[2] = b[2];
605 : 958 : a += 3;
606 : : }
607 : : }
608 : 385 : tk::destroy( m_bndpoinint );
609 : :
610 : : // Generate domain superedges
611 [ + - ]: 385 : domsuped();
612 : :
613 : : // Prepare symmetry boundary condition data structures
614 : :
615 : : // Query symmetry BC nodes associated to side sets
616 : 385 : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
617 [ + + ]: 677 : for (auto s : g_cfg.get< tag::bc_sym >()) {
618 [ + - ]: 292 : auto k = m_bface.find(s);
619 [ + + ]: 292 : if (k != end(m_bface)) {
620 [ + - ]: 110 : auto& n = sym[ k->first ];
621 [ + + ]: 1334 : for (auto f : k->second) {
622 : 1224 : const auto& t = m_triinpoel.data() + f*3;
623 [ + - ]: 1224 : n.insert( t[0] );
624 [ + - ]: 1224 : n.insert( t[1] );
625 [ + - ]: 1224 : n.insert( t[2] );
626 : : }
627 : : }
628 : : }
629 : :
630 : : // Generate unique set of symmetry BC nodes of all symmetryc BC side sets
631 : 385 : std::set< std::size_t > symbcnodeset;
632 [ + - ][ + + ]: 495 : for (const auto& [s,n] : sym) symbcnodeset.insert( begin(n), end(n) );
633 : :
634 : : // Generate symmetry BC data as streamable data structures
635 : 385 : tk::destroy( m_symbcnodes );
636 : 385 : tk::destroy( m_symbcnorms );
637 [ + + ]: 1603 : for (auto p : symbcnodeset) {
638 [ + + ]: 2436 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
639 [ + - ]: 1218 : auto m = m_bnorm.find( s );
640 [ + - ]: 1218 : if (m != end(m_bnorm)) {
641 [ + - ]: 1218 : auto r = m->second.find( p );
642 [ + - ]: 1218 : if (r != end(m->second)) {
643 [ + - ]: 1218 : m_symbcnodes.push_back( p );
644 [ + - ]: 1218 : m_symbcnorms.push_back( r->second[0] );
645 [ + - ]: 1218 : m_symbcnorms.push_back( r->second[1] );
646 [ + - ]: 1218 : m_symbcnorms.push_back( r->second[2] );
647 : : }
648 : : }
649 : : }
650 : : }
651 : 385 : tk::destroy( m_bnorm );
652 : :
653 : : // Prepare noslip boundary condition data structures
654 : :
655 : : // Query noslip BC nodes associated to side sets
656 : 385 : std::unordered_map< int, std::unordered_set< std::size_t > > noslip;
657 [ + + ]: 484 : for (auto s : g_cfg.get< tag::bc_noslip >()) {
658 [ + - ]: 99 : auto k = m_bface.find(s);
659 [ + + ]: 99 : if (k != end(m_bface)) {
660 [ + - ]: 89 : auto& n = noslip[ k->first ];
661 [ + + ]: 5141 : for (auto f : k->second) {
662 : 5052 : const auto& t = m_triinpoel.data() + f*3;
663 [ + - ]: 5052 : n.insert( t[0] );
664 [ + - ]: 5052 : n.insert( t[1] );
665 [ + - ]: 5052 : n.insert( t[2] );
666 : : }
667 : : }
668 : : }
669 : :
670 : : // Generate unique set of noslip BC nodes of all noslip BC side sets
671 : 385 : std::set< std::size_t > noslipbcnodeset;
672 [ + - ][ + + ]: 474 : for (const auto& [s,n] : noslip) noslipbcnodeset.insert( begin(n), end(n) );
673 : :
674 : : // Generate noslip BC data as streamable data structures
675 : 385 : tk::destroy( m_noslipbcnodes );
676 [ + - ]: 385 : m_noslipbcnodes.insert( m_noslipbcnodes.end(),
677 : : begin(noslipbcnodeset), end(noslipbcnodeset) );
678 : 385 : }
679 : :
680 : : void
681 : 385 : LohCG::domsuped()
682 : : // *****************************************************************************
683 : : // Generate superedge-groups for domain-edge loops
684 : : //! \see See Lohner, Sec. 15.1.6.2, An Introduction to Applied CFD Techniques,
685 : : //! Wiley, 2008.
686 : : // *****************************************************************************
687 : : {
688 [ - + ][ - - ]: 385 : Assert( !m_domedgeint.empty(), "No domain edges to group" );
[ - - ][ - - ]
689 : :
690 : : #ifndef NDEBUG
691 : 385 : auto nedge = m_domedgeint.size();
692 : : #endif
693 : :
694 [ + - ]: 385 : const auto& inpoel = Disc()->Inpoel();
695 [ + - ]: 385 : const auto& lid = Disc()->Lid();
696 [ + - ]: 385 : const auto& gid = Disc()->Gid();
697 : :
698 : 385 : tk::destroy( m_dsupedge[0] );
699 : 385 : tk::destroy( m_dsupedge[1] );
700 : 385 : tk::destroy( m_dsupedge[2] );
701 : :
702 : 385 : tk::destroy( m_dsupint[0] );
703 : 385 : tk::destroy( m_dsupint[1] );
704 : 385 : tk::destroy( m_dsupint[2] );
705 : :
706 : 385 : tk::UnsMesh::FaceSet untri;
707 [ + + ]: 54501 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
708 : : std::size_t N[4] = {
709 : 54116 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
710 [ + - ][ + + ]: 270580 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
711 : : }
712 : :
713 [ + + ]: 54501 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
714 : : std::size_t N[4] = {
715 : 54116 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
716 : 54116 : int f = 0;
717 : : tk::real sig[6];
718 [ + - ][ + + ]: 378812 : decltype(m_domedgeint)::const_iterator d[6];
719 [ + + ]: 163139 : for (const auto& [p,q] : tk::lpoed) {
720 : 154638 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
721 [ + + ]: 154638 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
722 [ + - ]: 154638 : d[f] = m_domedgeint.find( ed );
723 [ + + ]: 154638 : if (d[f] == end(m_domedgeint)) break; else ++f;
724 : : }
725 [ + + ]: 54116 : if (f == 6) {
726 [ + - ]: 8501 : m_dsupedge[0].push_back( N[0] );
727 [ + - ]: 8501 : m_dsupedge[0].push_back( N[1] );
728 [ + - ]: 8501 : m_dsupedge[0].push_back( N[2] );
729 [ + - ]: 8501 : m_dsupedge[0].push_back( N[3] );
730 [ + - ][ + + ]: 42505 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
731 [ + + ]: 59507 : for (int ed=0; ed<6; ++ed) {
732 : 51006 : const auto& ded = d[ed]->second;
733 [ + - ]: 51006 : m_dsupint[0].push_back( sig[ed] * ded[0] );
734 [ + - ]: 51006 : m_dsupint[0].push_back( sig[ed] * ded[1] );
735 [ + - ]: 51006 : m_dsupint[0].push_back( sig[ed] * ded[2] );
736 [ + - ]: 51006 : m_dsupint[0].push_back( ded[3] );
737 [ + - ]: 51006 : m_domedgeint.erase( d[ed] );
738 : : }
739 : : }
740 : : }
741 : :
742 [ + + ]: 96386 : for (const auto& N : untri) {
743 : 96001 : int f = 0;
744 : : tk::real sig[3];
745 [ + - ][ + + ]: 384004 : decltype(m_domedgeint)::const_iterator d[3];
746 [ + + ]: 160888 : for (const auto& [p,q] : tk::lpoet) {
747 : 152754 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
748 [ + + ]: 152754 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
749 [ + - ]: 152754 : d[f] = m_domedgeint.find( ed );
750 [ + + ]: 152754 : if (d[f] == end(m_domedgeint)) break; else ++f;
751 : : }
752 [ + + ]: 96001 : if (f == 3) {
753 [ + - ]: 8134 : m_dsupedge[1].push_back( N[0] );
754 [ + - ]: 8134 : m_dsupedge[1].push_back( N[1] );
755 [ + - ]: 8134 : m_dsupedge[1].push_back( N[2] );
756 [ + + ]: 32536 : for (int ed=0; ed<3; ++ed) {
757 : 24402 : const auto& ded = d[ed]->second;
758 [ + - ]: 24402 : m_dsupint[1].push_back( sig[ed] * ded[0] );
759 [ + - ]: 24402 : m_dsupint[1].push_back( sig[ed] * ded[1] );
760 [ + - ]: 24402 : m_dsupint[1].push_back( sig[ed] * ded[2] );
761 [ + - ]: 24402 : m_dsupint[1].push_back( ded[3] );
762 [ + - ]: 24402 : m_domedgeint.erase( d[ed] );
763 : : }
764 : : }
765 : : }
766 : :
767 [ + - ]: 385 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
768 [ + - ]: 385 : m_dsupint[2].resize( m_domedgeint.size()*4 );
769 : 385 : std::size_t k = 0;
770 [ + + ]: 23777 : for (const auto& [ed,d] : m_domedgeint) {
771 : 23392 : auto e = m_dsupedge[2].data() + k*2;
772 [ + - ]: 23392 : e[0] = tk::cref_find( lid, ed[0] );
773 [ + - ]: 23392 : e[1] = tk::cref_find( lid, ed[1] );
774 : 23392 : auto i = m_dsupint[2].data() + k*4;
775 : 23392 : i[0] = d[0];
776 : 23392 : i[1] = d[1];
777 : 23392 : i[2] = d[2];
778 : 23392 : i[3] = d[3];
779 : 23392 : ++k;
780 : : }
781 : :
782 : 385 : tk::destroy( m_domedgeint );
783 : :
784 : : //std::cout << std::setprecision(2)
785 : : // << "superedges: ntet:" << m_dsupedge[0].size()/4 << "(nedge:"
786 : : // << m_dsupedge[0].size()/4*6 << ","
787 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[0].size()/4*6 ) /
788 : : // static_cast< tk::real >( nedge )
789 : : // << "%) + ntri:" << m_dsupedge[1].size()/3
790 : : // << "(nedge:" << m_dsupedge[1].size() << ","
791 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[1].size() ) /
792 : : // static_cast< tk::real >( nedge )
793 : : // << "%) + nedge:"
794 : : // << m_dsupedge[2].size()/2 << "("
795 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[2].size()/2 ) /
796 : : // static_cast< tk::real >( nedge )
797 : : // << "%) = " << m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
798 : : // m_dsupedge[2].size()/2 << " of "<< nedge << " total edges\n";
799 : :
800 [ - + ][ - - ]: 385 : Assert( m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
[ - - ][ - - ]
801 : : m_dsupedge[2].size()/2 == nedge,
802 : : "Not all edges accounted for in superedge groups" );
803 : 385 : }
804 : :
805 : : void
806 : : // cppcheck-suppress unusedFunction
807 : 385 : LohCG::merge()
808 : : // *****************************************************************************
809 : : // Combine own and communicated portions of the integrals
810 : : // *****************************************************************************
811 : : {
812 : 385 : auto d = Disc();
813 : :
814 : : // Combine own and communicated contributions to boundary point normals
815 : 385 : bnorm();
816 : :
817 : : // Convert integrals into streamable data structures
818 : 385 : streamable();
819 : :
820 : : // Enforce boundary conditions on initial conditions
821 : 385 : auto t = d->T() + d->Dt();
822 : 385 : physics::dirbc( m_u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
823 : 385 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
824 : 385 : physics::noslipbc( m_u, m_noslipbcnodes, /*pos=*/1 );
825 : :
826 : : // Start measuring initial div-free time
827 : 385 : m_timer.emplace_back();
828 : :
829 : : // Compute initial momentum flux
830 [ + - ][ + - ]: 385 : thisProxy[ thisIndex ].wait4div();
831 [ + - ][ + - ]: 385 : thisProxy[ thisIndex ].wait4sgrad();
832 : 385 : div( m_u, /*pos=*/1 );
833 : 385 : }
834 : :
835 : : void
836 : 2195 : LohCG::fingrad( tk::Fields& grad,
837 : : std::unordered_map< std::size_t, std::vector< tk::real > >& gradc )
838 : : // *****************************************************************************
839 : : // Finalize computing gradient
840 : : //! \param[in,out] grad Gradient to finalize
841 : : //! \param[in,out] gradc Gradient communication buffer to finalize
842 : : // *****************************************************************************
843 : : {
844 [ + - ]: 2195 : auto d = Disc();
845 [ + - ]: 2195 : const auto lid = d->Lid();
846 : :
847 : : // Combine own and communicated contributions
848 [ + + ]: 91003 : for (const auto& [g,r] : gradc) {
849 [ + - ]: 88808 : auto i = tk::cref_find( lid, g );
850 [ + - ][ + + ]: 975248 : for (std::size_t c=0; c<r.size(); ++c) grad(i,c) += r[c];
851 : : }
852 : 2195 : tk::destroy(gradc);
853 : :
854 : : // Divide weak result by nodal volume
855 : 2195 : const auto& vol = d->Vol();
856 [ + + ]: 463926 : for (std::size_t p=0; p<grad.nunk(); ++p) {
857 [ + + ]: 5801506 : for (std::size_t c=0; c<grad.nprop(); ++c) {
858 [ + - ]: 5339775 : grad(p,c) /= vol[p];
859 : : }
860 : : }
861 : 2195 : }
862 : :
863 : : void
864 : 770 : LohCG::div( const tk::Fields& u, std::size_t pos )
865 : : // *****************************************************************************
866 : : // Start computing divergence
867 : : //! \param[in] u Vector field whose divergence to compute
868 : : //! \param[in] pos Position at which the three vector components are in u
869 : : // *****************************************************************************
870 : : {
871 [ + - ]: 770 : auto d = Disc();
872 [ + - ]: 770 : const auto lid = d->Lid();
873 : :
874 : : // Finalize momentum flux communications if needed
875 [ + + ]: 770 : if (m_np == 1) {
876 [ + - ]: 385 : fingrad( m_flux, m_fluxc );
877 [ + - ]: 385 : physics::symbc( m_flux, m_symbcnodes, m_symbcnorms, /*pos=*/0 );
878 : : }
879 : :
880 : : // Compute velocity divergence
881 [ + - ]: 770 : std::fill( begin(m_div), end(m_div), 0.0 );
882 [ + - ]: 770 : lohner::div( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, u, m_div, pos );
883 : :
884 : : // Communicate velocity divergence to other chares on chare-boundary
885 [ + + ]: 770 : if (d->NodeCommMap().empty()) {
886 [ + - ]: 18 : comdiv_complete();
887 : : } else {
888 [ + + ]: 8460 : for (const auto& [c,n] : d->NodeCommMap()) {
889 : 7708 : decltype(m_divc) exp;
890 [ + - ][ + - ]: 42152 : for (auto g : n) exp[g] = m_div[ tk::cref_find(lid,g) ];
[ + + ]
891 [ + - ][ + - ]: 7708 : thisProxy[c].comdiv( exp );
892 : 7708 : }
893 : : }
894 [ + - ]: 770 : owndiv_complete();
895 : 770 : }
896 : :
897 : : void
898 : 7708 : LohCG::comdiv( const std::unordered_map< std::size_t, tk::real >& indiv )
899 : : // *****************************************************************************
900 : : // Receive contributions to velocity divergence on chare-boundaries
901 : : //! \param[in] indiv Partial contributions of velocity divergence to
902 : : //! chare-boundary nodes. Key: global mesh node IDs, value: contribution.
903 : : //! \details This function receives contributions to m_div, which stores the
904 : : //! velocity divergence at mesh nodes. While m_div stores own contributions,
905 : : //! m_divc collects the neighbor chare contributions during communication.
906 : : //! This way work on m_div and m_divc is overlapped.
907 : : // *****************************************************************************
908 : : {
909 [ + - ][ + + ]: 42152 : for (const auto& [g,r] : indiv) m_divc[g] += r;
910 : :
911 : : // When we have heard from all chares we communicate with, this chare is done
912 [ + + ]: 7708 : if (++m_ndiv == Disc()->NodeCommMap().size()) {
913 : 752 : m_ndiv = 0;
914 : 752 : comdiv_complete();
915 : : }
916 : 7708 : }
917 : :
918 : : void
919 : 385 : LohCG::velgrad()
920 : : // *****************************************************************************
921 : : // Start computing velocity gradient
922 : : // *****************************************************************************
923 : : {
924 : 385 : auto d = Disc();
925 : :
926 : : // Compute momentum flux
927 : 385 : m_vgrad.fill( 0.0 );
928 : 385 : lohner::vgrad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad );
929 : :
930 : : // Communicate velocity divergence to other chares on chare-boundary
931 : 385 : const auto& lid = d->Lid();
932 [ + + ]: 385 : if (d->NodeCommMap().empty()) {
933 : 9 : comvgrad_complete();
934 : : } else {
935 [ + + ]: 4230 : for (const auto& [c,n] : d->NodeCommMap()) {
936 : 3854 : decltype(m_vgradc) exp;
937 [ + - ][ + - ]: 21076 : for (auto g : n) exp[g] = m_vgrad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
938 [ + - ][ + - ]: 3854 : thisProxy[c].comvgrad( exp );
939 : 3854 : }
940 : : }
941 : 385 : ownvgrad_complete();
942 : 385 : }
943 : :
944 : : void
945 : 3854 : LohCG::comvgrad(
946 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
947 : : // *****************************************************************************
948 : : // Receive contributions to velocity gradients on chare-boundaries
949 : : //! \param[in] ingrad Partial contributions of momentum flux to
950 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
951 : : //! \details This function receives contributions to m_vgrad, which stores the
952 : : //! velocity gradients at mesh nodes. While m_vgrad stores own contributions,
953 : : //! m_vgradc collects the neighbor chare contributions during communication.
954 : : //! This way work on m_vgrad and m_vgradc is overlapped.
955 : : // *****************************************************************************
956 : : {
957 : : using tk::operator+=;
958 [ + - ][ + - ]: 21076 : for (const auto& [g,r] : ingrad) m_vgradc[g] += r;
[ + + ]
959 : :
960 : : // When we have heard from all chares we communicate with, this chare is done
961 [ + + ]: 3854 : if (++m_nvgrad == Disc()->NodeCommMap().size()) {
962 : 376 : m_nvgrad = 0;
963 : 376 : comvgrad_complete();
964 : : }
965 : 3854 : }
966 : :
967 : : void
968 : 385 : LohCG::flux()
969 : : // *****************************************************************************
970 : : // Start computing momentum flux
971 : : // *****************************************************************************
972 : : {
973 : 385 : auto d = Disc();
974 : :
975 : : // Finalize computing velocity gradients
976 : 385 : fingrad( m_vgrad, m_vgradc );
977 : :
978 : : // Compute momentum flux
979 : 385 : m_flux.fill( 0.0 );
980 : 385 : lohner::flux( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_vgrad,
981 : 385 : m_flux );
982 : :
983 : : // Communicate velocity divergence to other chares on chare-boundary
984 : 385 : const auto& lid = d->Lid();
985 [ + + ]: 385 : if (d->NodeCommMap().empty()) {
986 : 9 : comflux_complete();
987 : : } else {
988 [ + + ]: 4230 : for (const auto& [c,n] : d->NodeCommMap()) {
989 : 3854 : decltype(m_fluxc) exp;
990 [ + - ][ + - ]: 21076 : for (auto g : n) exp[g] = m_flux[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
991 [ + - ][ + - ]: 3854 : thisProxy[c].comflux( exp );
992 : 3854 : }
993 : : }
994 : 385 : ownflux_complete();
995 : 385 : }
996 : :
997 : : void
998 : 3854 : LohCG::comflux(
999 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& influx )
1000 : : // *****************************************************************************
1001 : : // Receive contributions to momentum flux on chare-boundaries
1002 : : //! \param[in] influx Partial contributions of momentum flux to
1003 : : //! chare-boundary nodes. Key: global mesh node IDs, values: contributions.
1004 : : //! \details This function receives contributions to m_flux, which stores the
1005 : : //! momentum flux at mesh nodes. While m_flux stores own contributions,
1006 : : //! m_fluxc collects the neighbor chare contributions during communication.
1007 : : //! This way work on m_flux and m_fluxc is overlapped.
1008 : : // *****************************************************************************
1009 : : {
1010 : : using tk::operator+=;
1011 [ + - ][ + - ]: 21076 : for (const auto& [g,r] : influx) m_fluxc[g] += r;
[ + + ]
1012 : :
1013 : : // When we have heard from all chares we communicate with, this chare is done
1014 [ + + ]: 3854 : if (++m_nflux == Disc()->NodeCommMap().size()) {
1015 : 376 : m_nflux = 0;
1016 : 376 : comflux_complete();
1017 : : }
1018 : 3854 : }
1019 : :
1020 : : void
1021 : 770 : LohCG::pinit()
1022 : : // *****************************************************************************
1023 : : // Initialize Poisson solve
1024 : : // *****************************************************************************
1025 : : {
1026 [ + - ]: 770 : auto d = Disc();
1027 [ + - ]: 770 : const auto lid = d->Lid();
1028 : 770 : const auto& coord = d->Coord();
1029 : 770 : const auto& x = coord[0];
1030 : 770 : const auto& y = coord[1];
1031 : 770 : const auto& z = coord[2];
1032 : :
1033 : : // Combine own and communicated contributions to velocity divergence
1034 [ + - ][ + + ]: 17842 : for (const auto& [g,r] : m_divc) m_div[ tk::cref_find(lid,g) ] += r;
1035 : 770 : tk::destroy(m_divc);
1036 : :
1037 : : // Configure Dirichlet BCs
1038 : : std::unordered_map< std::size_t,
1039 : 770 : std::vector< std::pair< int, tk::real > > > dirbc;
1040 [ + + ]: 770 : if (!g_cfg.get< tag::pre_bc_dir >().empty()) {
1041 [ + - ]: 70 : auto ic = problems::PRESSURE_IC();
1042 : 70 : std::size_t nmask = 1 + 1;
1043 [ - + ][ - - ]: 70 : Assert( m_dirbcmaskp.size() % nmask == 0, "Size mismatch" );
[ - - ][ - - ]
1044 [ + + ]: 500 : for (std::size_t i=0; i<m_dirbcmaskp.size()/nmask; ++i) {
1045 : 430 : auto p = m_dirbcmaskp[i*nmask+0]; // local node id
1046 : 430 : auto mask = m_dirbcmaskp[i*nmask+1];
1047 [ + + ]: 430 : if (mask == 1) { // mask == 1: IC value
1048 [ + - ]: 190 : auto val = ic( x[p], y[p], z[p] );
1049 [ + - ][ + - ]: 190 : dirbc[p] = {{ { 1, val } }};
1050 [ + - ][ + - ]: 240 : } else if (mask == 2 && !m_dirbcvalp.empty()) { // mask == 2: BC value
[ + - ]
1051 : 240 : auto val = m_dirbcvalp[i*nmask+1];
1052 [ + - ][ + - ]: 240 : dirbc[p] = {{ { 1, val } }};
1053 : : }
1054 : : }
1055 : 70 : }
1056 : :
1057 : : // Configure Neumann BCs
1058 : 770 : std::vector< tk::real > neubc;
1059 [ + - ]: 770 : auto pg = problems::PRESSURE_GRAD();
1060 [ - + ]: 770 : if (pg) {
1061 : : // Collect Neumann BC elements
1062 [ - - ]: 0 : std::vector< std::uint8_t > besym( m_triinpoel.size(), 0 );
1063 [ - - ]: 0 : for (auto s : g_cfg.get< tag::pre_bc_sym >()) {
1064 [ - - ]: 0 : auto k = m_bface.find(s);
1065 [ - - ][ - - ]: 0 : if (k != end(m_bface)) for (auto f : k->second) besym[f] = 1;
1066 : : }
1067 : : // Setup Neumann BCs
1068 [ - - ]: 0 : neubc.resize( x.size(), 0.0 );
1069 [ - - ]: 0 : for (std::size_t e=0; e<m_triinpoel.size()/3; ++e) {
1070 [ - - ]: 0 : if (besym[e]) {
1071 : 0 : const auto N = m_triinpoel.data() + e*3;
1072 : : tk::real n[3];
1073 : 0 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
1074 : 0 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
1075 : : n[0], n[1], n[2] );
1076 [ - - ]: 0 : auto g = pg( x[N[0]], y[N[0]], z[N[0]] );
1077 : 0 : neubc[ N[0] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1078 [ - - ]: 0 : g = pg( x[N[1]], y[N[1]], z[N[1]] );
1079 : 0 : neubc[ N[1] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1080 [ - - ]: 0 : g = pg( x[N[2]], y[N[2]], z[N[2]] );
1081 : 0 : neubc[ N[2] ] -= n[0]*g[0] + n[1]*g[1] + n[2]*g[2];
1082 : : }
1083 : : }
1084 : 0 : }
1085 : :
1086 : : // Set hydrostat
1087 : 770 : auto h = g_cfg.get< tag::pre_hydrostat >();
1088 [ + + ]: 770 : if (h != std::numeric_limits< uint64_t >::max()) {
1089 [ + - ]: 116 : auto pi = lid.find( h );
1090 [ + + ]: 116 : if (pi != end(lid)) {
1091 : 16 : auto p = pi->second;
1092 [ + - ]: 16 : auto ic = problems::PRESSURE_IC();
1093 [ + - ][ + - ]: 16 : auto val = m_np>1 ? 0.0 : ic( x[p], y[p], z[p] );
1094 [ + - ]: 16 : auto& b = dirbc[p];
1095 [ + - ][ + - ]: 16 : if (b.empty()) b = {{ { 1, val }} };
1096 : 16 : }
1097 : : }
1098 : :
1099 : : // Configure right hand side
1100 [ + - ]: 770 : auto pr = problems::PRESSURE_RHS();
1101 [ - + ]: 770 : if (pr) {
1102 : 0 : const auto& vol = d->Vol();
1103 [ - - ]: 0 : for (std::size_t i=0; i<x.size(); ++i) {
1104 [ - - ]: 0 : m_div[i] = pr( x[i], y[i], z[i] ) * vol[i];
1105 : : }
1106 : : }
1107 : :
1108 : : // Initialize Poisson solve
1109 : 770 : const auto& pc = g_cfg.get< tag::pre_pc >();
1110 [ + - ][ + - ]: 1540 : m_cgpre[ thisIndex ].ckLocal()->init( {}, m_div, neubc, dirbc, true, pc,
[ + - ]
1111 [ + - ][ + - ]: 1540 : CkCallback( CkIndex_LohCG::psolve(), thisProxy[thisIndex] ) );
[ + - ]
1112 : 770 : }
1113 : :
1114 : : void
1115 : 770 : LohCG::psolve()
1116 : : // *****************************************************************************
1117 : : // Solve Poisson equation
1118 : : // *****************************************************************************
1119 : : {
1120 : 770 : auto iter = g_cfg.get< tag::pre_iter >();
1121 : 770 : auto tol = g_cfg.get< tag::pre_tol >();
1122 : 770 : auto verbose = g_cfg.get< tag::pre_verbose >();
1123 : :
1124 : 770 : auto c = m_np != 1 ?
1125 : 385 : CkCallback( CkIndex_LohCG::sgrad(), thisProxy[thisIndex] ) :
1126 [ + + ][ + - ]: 770 : CkCallback( CkIndex_LohCG::psolved(), thisProxy[thisIndex] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ][ - - ]
[ - - ]
1127 : :
1128 [ + - ][ + - ]: 770 : m_cgpre[ thisIndex ].ckLocal()->solve( iter, tol, thisIndex, verbose, c );
[ + - ]
1129 : 770 : }
1130 : :
1131 : : void
1132 : 385 : LohCG::sgrad()
1133 : : // *****************************************************************************
1134 : : // Compute recent conjugate gradients solution gradient
1135 : : // *****************************************************************************
1136 : : {
1137 [ + - ]: 385 : auto d = Disc();
1138 : :
1139 [ + - ][ + - ]: 385 : auto sol = m_cgpre[ thisIndex ].ckLocal()->solution();
[ + - ]
1140 [ + - ]: 385 : m_sgrad.fill( 0.0 );
1141 [ + - ]: 385 : lohner::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, sol, m_sgrad );
1142 : :
1143 : : // Send gradient contributions to neighbor chares
1144 [ + + ]: 385 : if (d->NodeCommMap().empty()) {
1145 [ + - ]: 9 : comsgrad_complete();
1146 : : } else {
1147 : 376 : const auto& lid = d->Lid();
1148 [ + + ]: 4230 : for (const auto& [c,n] : d->NodeCommMap()) {
1149 : 3854 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1150 [ + - ][ + - ]: 21076 : for (auto g : n) exp[g] = m_sgrad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1151 [ + - ][ + - ]: 3854 : thisProxy[c].comsgrad( exp );
1152 : 3854 : }
1153 : : }
1154 [ + - ]: 385 : ownsgrad_complete();
1155 : 385 : }
1156 : :
1157 : : void
1158 : 3854 : LohCG::comsgrad(
1159 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
1160 : : // *****************************************************************************
1161 : : // Receive contributions to conjugrate gradients solution gradient
1162 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
1163 : : //! global mesh node IDs, value: contributions for all scalar components.
1164 : : //! \details This function receives contributions to m_sgrad, which stores the
1165 : : //! gradients at mesh nodes. While m_sgrad stores own contributions, m_sgradc
1166 : : //! collects the neighbor chare contributions during communication. This way
1167 : : //! work on m_sgrad and m_sgradc is overlapped.
1168 : : // *****************************************************************************
1169 : : {
1170 : : using tk::operator+=;
1171 [ + - ][ + - ]: 21076 : for (const auto& [g,r] : ingrad) m_sgradc[g] += r;
[ + + ]
1172 : :
1173 : : // When we have heard from all chares we communicate with, this chare is done
1174 [ + + ]: 3854 : if (++m_nsgrad == Disc()->NodeCommMap().size()) {
1175 : 376 : m_nsgrad = 0;
1176 : 376 : comsgrad_complete();
1177 : : }
1178 : 3854 : }
1179 : :
1180 : : void
1181 : 770 : LohCG::psolved()
1182 : : // *****************************************************************************
1183 : : // Continue setup after Poisson solve and gradient computation
1184 : : // *****************************************************************************
1185 : : {
1186 : 770 : auto d = Disc();
1187 : :
1188 [ + + ][ + - ]: 770 : if (thisIndex == 0) d->pit( m_cgpre[ thisIndex ].ckLocal()->it() );
[ + - ][ + - ]
1189 : :
1190 [ + + ]: 770 : if (m_np != 1) {
1191 : : // Finalize gradient communications
1192 : 385 : fingrad( m_sgrad, m_sgradc );
1193 : : // Project velocity to divergence-free subspace
1194 [ + + ]: 23682 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1195 : 23297 : m_u(i,1) -= m_sgrad(i,0);
1196 : 23297 : m_u(i,2) -= m_sgrad(i,1);
1197 : 23297 : m_u(i,3) -= m_sgrad(i,2);
1198 : : }
1199 : : // Enforce boundary conditions
1200 : 385 : auto t = d->T() + d->Dt();
1201 : 385 : physics::dirbc( m_u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval);
1202 : 385 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
1203 : 385 : physics::noslipbc( m_u, m_noslipbcnodes, /*pos=*/1 );
1204 : : }
1205 : :
1206 [ - + ]: 770 : if (g_cfg.get< tag::nstep >() == 1) { // test first Poisson solve only
1207 : :
1208 [ - - ][ - - ]: 0 : auto p = m_cgpre[ thisIndex ].ckLocal()->solution();
[ - - ]
1209 [ - - ][ - - ]: 0 : for (std::size_t i=0; i<m_u.nunk(); ++i) m_u(i,0) = p[i];
1210 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4step();
1211 [ - - ][ - - ]: 0 : writeFields( CkCallback(CkIndex_LohCG::diag(), thisProxy[thisIndex]) );
[ - - ][ - - ]
1212 : :
1213 : 0 : } else {
1214 : :
1215 [ + + ]: 770 : if (++m_np < 2) {
1216 : : // Compute momentum flux for initial pressure-Poisson rhs
1217 [ + - ][ + - ]: 385 : thisProxy[ thisIndex ].wait4vgrad();
1218 [ + - ][ + - ]: 385 : thisProxy[ thisIndex ].wait4flux();
1219 [ + - ][ + - ]: 385 : thisProxy[ thisIndex ].wait4div();
1220 : 385 : velgrad();
1221 : : } else {
1222 [ + + ]: 385 : if (thisIndex == 0) {
1223 [ + - ]: 99 : tk::Print() << "Initial div-free time: " << m_timer[0].dsec()
1224 [ + - ][ + - ]: 66 : << " sec\n";
[ + - ]
1225 : : }
1226 : : // Assign initial pressure and start timestepping
1227 [ + - ][ + - ]: 385 : auto p = m_cgpre[ thisIndex ].ckLocal()->solution();
[ + - ]
1228 [ + - ][ + + ]: 23682 : for (std::size_t i=0; i<m_u.nunk(); ++i) m_u(i,0) = p[i];
1229 [ + - ][ + - ]: 385 : writeFields( CkCallback(CkIndex_LohCG::start(), thisProxy[thisIndex]) );
[ + - ][ + - ]
1230 : 385 : }
1231 : :
1232 : : }
1233 : 770 : }
1234 : :
1235 : : void
1236 : 4730 : LohCG::diag()
1237 : : // *****************************************************************************
1238 : : // Compute diagnostics
1239 : : // *****************************************************************************
1240 : : {
1241 : 4730 : auto d = Disc();
1242 : :
1243 : : // Increase number of iterations and physical time
1244 : 4730 : d->next();
1245 : :
1246 : : // Compute diagnostics, e.g., residuals
1247 : 4730 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1248 : 4730 : auto diagnostics = m_diag.accompute( *d, m_u, m_un, diag_iter );
1249 : :
1250 : : // Evaluate residuals
1251 [ + + ][ + - ]: 4730 : if (!diagnostics) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ + - ]
1252 : 4730 : }
1253 : :
1254 : : void
1255 : 385 : LohCG::start()
1256 : : // *****************************************************************************
1257 : : // Start time stepping
1258 : : // *****************************************************************************
1259 : : {
1260 : : // Set flag that indicates that we are now during time stepping
1261 : 385 : Disc()->Initial( 0 );
1262 : : // Start timer measuring time stepping wall clock time
1263 : 385 : Disc()->Timer().zero();
1264 : : // Zero grind-timer
1265 : 385 : Disc()->grindZero();
1266 : : // Continue to first time step
1267 : 385 : dt();
1268 : 385 : }
1269 : :
1270 : : void
1271 : 4730 : LohCG::dt()
1272 : : // *****************************************************************************
1273 : : // Compute time step size
1274 : : // *****************************************************************************
1275 : : {
1276 [ + - ]: 4730 : auto d = Disc();
1277 : 4730 : const auto& vol = d->Vol();
1278 : :
1279 : 4730 : tk::real mindt = std::numeric_limits< tk::real >::max();
1280 : 4730 : auto const_dt = g_cfg.get< tag::dt >();
1281 : 4730 : auto eps = std::numeric_limits< tk::real >::epsilon();
1282 : :
1283 : : // use constant dt if configured
1284 [ + + ]: 4730 : if (std::abs(const_dt) > eps) {
1285 : :
1286 : : // cppcheck-suppress redundantInitialization
1287 : 2920 : mindt = const_dt;
1288 : :
1289 : : } else {
1290 : :
1291 : 1810 : auto cfl = g_cfg.get< tag::cfl >();
1292 : 1810 : auto large = std::numeric_limits< tk::real >::max();
1293 : 1810 : auto c = g_cfg.get< tag::soundspeed >();
1294 : 1810 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1295 : 1810 : auto dif = g_cfg.get< tag::mat_dyn_diffusivity >();
1296 : 1810 : dif = std::max( mu, dif );
1297 : :
1298 [ + + ]: 308220 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1299 [ + - ]: 306410 : auto u = m_u(i,1);
1300 [ + - ]: 306410 : auto v = m_u(i,2);
1301 [ + - ]: 306410 : auto w = m_u(i,3);
1302 : 306410 : auto vel = std::sqrt( u*u + v*v + w*w );
1303 : 306410 : auto L = std::cbrt( vol[i] );
1304 : 306410 : auto euler_dt = L / std::max( vel+c, 1.0e-8 );
1305 : 306410 : mindt = std::min( mindt, euler_dt );
1306 [ + + ]: 306410 : auto visc_dt = dif > eps ? L * L / dif : large;
1307 : 306410 : mindt = std::min( mindt, visc_dt );
1308 : : }
1309 : 1810 : mindt *= cfl;
1310 : :
1311 : : }
1312 : :
1313 : : // Actiavate SDAG waits for next time step stage
1314 [ + - ][ + - ]: 4730 : thisProxy[ thisIndex ].wait4step();
1315 : :
1316 : : // Contribute to minimum dt across all chares and advance to next step
1317 [ + - ]: 4730 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
1318 [ + - ][ + - ]: 9460 : CkCallback(CkReductionTarget(LohCG,advance), thisProxy) );
1319 : 4730 : }
1320 : :
1321 : : void
1322 : 4730 : LohCG::advance( tk::real newdt )
1323 : : // *****************************************************************************
1324 : : // Advance equations to next time step
1325 : : //! \param[in] newdt The smallest dt across the whole problem
1326 : : // *****************************************************************************
1327 : : {
1328 : : // Detect blowup
1329 : 4730 : auto eps = std::numeric_limits< tk::real >::epsilon();
1330 [ - + ]: 4730 : if (newdt < eps) m_finished = 1;
1331 : :
1332 : : // Set new time step size
1333 : 4730 : Disc()->setdt( newdt );
1334 : :
1335 : : // Start next time step stage
1336 : 4730 : stage();
1337 : 4730 : }
1338 : :
1339 : : void
1340 : 8600 : LohCG::stage()
1341 : : // *****************************************************************************
1342 : : // Start next time step stage
1343 : : // *****************************************************************************
1344 : : {
1345 : : // Activate SDAG waits for next time step stage
1346 [ + - ][ + - ]: 8600 : thisProxy[ thisIndex ].wait4grad();
1347 [ + - ][ + - ]: 8600 : thisProxy[ thisIndex ].wait4rhs();
1348 : :
1349 : : // Compute gradients
1350 : 8600 : grad();
1351 : 8600 : }
1352 : :
1353 : : void
1354 : 8600 : LohCG::grad()
1355 : : // *****************************************************************************
1356 : : // Compute gradients
1357 : : // *****************************************************************************
1358 : : {
1359 : 8600 : auto d = Disc();
1360 : :
1361 [ + + ]: 8600 : if (m_grad.nprop()) {
1362 : 1040 : m_grad.fill( 0.0 );
1363 : 1040 : lohner::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_grad );
1364 : : }
1365 : :
1366 : : // Send gradient contributions to neighbor chares
1367 [ + + ]: 8600 : if (d->NodeCommMap().empty()) {
1368 : 460 : comgrad_complete();
1369 : : } else {
1370 : 8140 : const auto& lid = d->Lid();
1371 [ + + ]: 66060 : for (const auto& [c,n] : d->NodeCommMap()) {
1372 : 57920 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1373 [ + - ][ + - ]: 454600 : for (auto g : n) exp[g] = m_grad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1374 [ + - ][ + - ]: 57920 : thisProxy[c].comgrad( exp );
1375 : 57920 : }
1376 : : }
1377 : 8600 : owngrad_complete();
1378 : 8600 : }
1379 : :
1380 : : void
1381 : 57920 : LohCG::comgrad(
1382 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
1383 : : // *****************************************************************************
1384 : : // Receive contributions to gradients on chare-boundaries
1385 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
1386 : : //! global mesh node IDs, value: contributions for all scalar components.
1387 : : //! \details This function receives contributions to m_grad, which stores the
1388 : : //! gradients at mesh nodes. While m_grad stores own contributions, m_gradc
1389 : : //! collects the neighbor chare contributions during communication. This way
1390 : : //! work on m_grad and m_gradc is overlapped.
1391 : : // *****************************************************************************
1392 : : {
1393 : : using tk::operator+=;
1394 [ + - ][ + - ]: 454600 : for (const auto& [g,r] : ingrad) m_gradc[g] += r;
[ + + ]
1395 : :
1396 : : // When we have heard from all chares we communicate with, this chare is done
1397 [ + + ]: 57920 : if (++m_ngrad == Disc()->NodeCommMap().size()) {
1398 : 8140 : m_ngrad = 0;
1399 : 8140 : comgrad_complete();
1400 : : }
1401 : 57920 : }
1402 : :
1403 : : void
1404 : 8600 : LohCG::rhs()
1405 : : // *****************************************************************************
1406 : : // Compute right-hand side of transport equations
1407 : : // *****************************************************************************
1408 : : {
1409 : 8600 : auto d = Disc();
1410 : 8600 : const auto& lid = d->Lid();
1411 : :
1412 : : // Combine own and communicated contributions to gradients
1413 [ + + ]: 8600 : if (m_grad.nprop()) fingrad( m_grad, m_gradc );
1414 : :
1415 : : // Compute own portion of right-hand side for all equations
1416 : 8600 : lohner::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, d->V(), d->T(),
1417 : 8600 : m_u, m_grad, m_rhs );
1418 : :
1419 : : // Communicate rhs to other chares on chare-boundary
1420 [ + + ]: 8600 : if (d->NodeCommMap().empty()) {
1421 : 460 : comrhs_complete();
1422 : : } else {
1423 [ + + ]: 66060 : for (const auto& [c,n] : d->NodeCommMap()) {
1424 : 57920 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1425 [ + - ][ + - ]: 454600 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1426 [ + - ][ + - ]: 57920 : thisProxy[c].comrhs( exp );
1427 : 57920 : }
1428 : : }
1429 : 8600 : ownrhs_complete();
1430 : 8600 : }
1431 : :
1432 : : void
1433 : 57920 : LohCG::comrhs(
1434 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1435 : : // *****************************************************************************
1436 : : // Receive contributions to right-hand side vector on chare-boundaries
1437 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1438 : : //! global mesh node IDs, value: contributions for all scalar components.
1439 : : //! \details This function receives contributions to m_rhs, which stores the
1440 : : //! right hand side vector at mesh nodes. While m_rhs stores own
1441 : : //! contributions, m_rhsc collects the neighbor chare contributions during
1442 : : //! communication. This way work on m_rhs and m_rhsc is overlapped.
1443 : : // *****************************************************************************
1444 : : {
1445 : : using tk::operator+=;
1446 [ + - ][ + - ]: 454600 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
[ + + ]
1447 : :
1448 : : // When we have heard from all chares we communicate with, this chare is done
1449 [ + + ]: 57920 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
1450 : 8140 : m_nrhs = 0;
1451 : 8140 : comrhs_complete();
1452 : : }
1453 : 57920 : }
1454 : :
1455 : : void
1456 : : // cppcheck-suppress unusedFunction
1457 : 8600 : LohCG::solve()
1458 : : // *****************************************************************************
1459 : : // Advance systems of equations
1460 : : // *****************************************************************************
1461 : : {
1462 [ + - ]: 8600 : auto d = Disc();
1463 : 8600 : const auto npoin = m_u.nunk();
1464 : 8600 : const auto ncomp = m_u.nprop();
1465 : 8600 : const auto& vol = d->Vol();
1466 : :
1467 : : // Combine own and communicated contributions to rhs
1468 [ + - ]: 8600 : const auto lid = d->Lid();
1469 [ + + ]: 266660 : for (const auto& [g,r] : m_rhsc) {
1470 [ + - ]: 258060 : auto i = tk::cref_find( lid, g );
1471 [ + - ][ + + ]: 1395720 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1472 : : }
1473 : 8600 : tk::destroy(m_rhsc);
1474 : :
1475 [ + + ][ + - ]: 8600 : if (m_stage == 0) m_un = m_u;
1476 : :
1477 : : // Advance system
1478 : 8600 : auto dt = m_rkcoef[m_stage] * d->Dt();
1479 [ + + ]: 1098000 : for (std::size_t i=0; i<npoin; ++i) {
1480 [ + + ]: 5816140 : for (std::size_t c=0; c<ncomp; ++c) {
1481 [ + - ][ + - ]: 4726740 : m_u(i,c) = m_un(i,c) - dt*m_rhs(i,c)/vol[i];
[ + - ]
1482 : : }
1483 : : }
1484 : :
1485 : : // Configure and apply scalar source to solution (if defined)
1486 [ + - ]: 8600 : auto src = problems::PHYS_SRC();
1487 [ - + ][ - - ]: 8600 : if (src) src( d->Coord(), d->T(), m_u );
1488 : :
1489 : : // Enforce boundary conditions
1490 : 8600 : auto t = d->T() + m_rkcoef[m_stage] * d->Dt();
1491 [ + - ]: 8600 : physics::dirbc( m_u, t, d->Coord(), d->BoxNodes(), m_dirbcmask, m_dirbcval );
1492 [ + - ]: 8600 : physics::dirbcp( m_u, d->Coord(), m_dirbcmaskp, m_dirbcvalp );
1493 [ + - ]: 8600 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
1494 [ + - ]: 8600 : physics::noslipbc( m_u, m_noslipbcnodes, /*pos=*/1 );
1495 : :
1496 [ + + ]: 8600 : if (++m_stage < m_rkcoef.size()) {
1497 : :
1498 : : // Start next time step stage
1499 [ + - ]: 3870 : stage();
1500 : :
1501 : : } else {
1502 : :
1503 : : // Reset Runge-Kutta stage counter
1504 : 4730 : m_stage = 0;
1505 : : // Compute diagnostics, e.g., residuals
1506 [ + - ]: 4730 : diag();
1507 : :
1508 : : }
1509 : 8600 : }
1510 : :
1511 : : void
1512 : 4730 : LohCG::evalres( const std::vector< tk::real >& )
1513 : : // *****************************************************************************
1514 : : // Evaluate residuals
1515 : : // *****************************************************************************
1516 : : {
1517 : 4730 : refine();
1518 : 4730 : }
1519 : :
1520 : : void
1521 : 4730 : LohCG::refine()
1522 : : // *****************************************************************************
1523 : : // Optionally refine/derefine mesh
1524 : : // *****************************************************************************
1525 : : {
1526 : 4730 : auto d = Disc();
1527 : :
1528 : : // See if this is the last time step
1529 [ + + ]: 4730 : if (d->finished()) m_finished = 1;
1530 : :
1531 : 4730 : auto dtref = g_cfg.get< tag::href_dt >();
1532 : 4730 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
1533 : :
1534 : : // if t>0 refinement enabled and we hit the frequency
1535 [ - + ][ - - ]: 4730 : if (dtref && !(d->It() % dtfreq)) { // refine
[ - + ]
1536 : :
1537 : 0 : d->refined() = 1;
1538 : 0 : d->startvol();
1539 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1540 : :
1541 : : // Activate SDAG waits for re-computing the integrals
1542 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1543 : :
1544 : : } else { // do not refine
1545 : :
1546 : 4730 : d->refined() = 0;
1547 : 4730 : feop_complete();
1548 : 4730 : resize_complete();
1549 : :
1550 : : }
1551 : 4730 : }
1552 : :
1553 : : void
1554 : 0 : LohCG::resizePostAMR(
1555 : : const std::vector< std::size_t >& /*ginpoel*/,
1556 : : const tk::UnsMesh::Chunk& chunk,
1557 : : const tk::UnsMesh::Coords& coord,
1558 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1559 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1560 : : const std::set< std::size_t >& removedNodes,
1561 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1562 : : nodeCommMap,
1563 : : const std::map< int, std::vector< std::size_t > >& bface,
1564 : : const std::map< int, std::vector< std::size_t > >& bnode,
1565 : : const std::vector< std::size_t >& triinpoel )
1566 : : // *****************************************************************************
1567 : : // Receive new mesh from Refiner
1568 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1569 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1570 : : //! \param[in] coord New mesh node coordinates
1571 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1572 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1573 : : //! \param[in] removedNodes Newly removed mesh node local ids
1574 : : //! \param[in] nodeCommMap New node communication map
1575 : : //! \param[in] bface Boundary-faces mapped to side set ids
1576 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1577 : : //! \param[in] triinpoel Boundary-face connectivity
1578 : : // *****************************************************************************
1579 : : {
1580 [ - - ]: 0 : auto d = Disc();
1581 : :
1582 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1583 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1584 : :
1585 : : // Resize mesh data structures after mesh refinement
1586 [ - - ]: 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1587 : :
1588 [ - - ][ - - ]: 0 : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1589 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1590 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1591 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1592 : :
1593 : : // Remove newly removed nodes from solution vectors
1594 [ - - ]: 0 : m_u.rm( removedNodes );
1595 [ - - ]: 0 : m_rhs.rm( removedNodes );
1596 : :
1597 : : // Resize auxiliary solution vectors
1598 : 0 : auto npoin = coord[0].size();
1599 [ - - ]: 0 : m_u.resize( npoin );
1600 [ - - ]: 0 : m_rhs.resize( npoin );
1601 : :
1602 : : // Update solution on new mesh
1603 [ - - ]: 0 : for (const auto& n : addedNodes)
1604 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1605 [ - - ][ - - ]: 0 : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
[ - - ][ - - ]
1606 [ - - ][ - - ]: 0 : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
[ - - ][ - - ]
[ - - ]
1607 : : "Indices of parent-edge nodes out of bounds post-AMR");
1608 [ - - ][ - - ]: 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
[ - - ]
1609 : : }
1610 : :
1611 : : // Update physical-boundary node-, face-, and element lists
1612 [ - - ]: 0 : m_bnode = bnode;
1613 [ - - ]: 0 : m_bface = bface;
1614 [ - - ]: 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1615 : :
1616 : 0 : auto meshid = d->MeshId();
1617 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1618 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1619 : 0 : }
1620 : :
1621 : : void
1622 : 914 : LohCG::writeFields( CkCallback cb )
1623 : : // *****************************************************************************
1624 : : // Output mesh-based fields to file
1625 : : //! \param[in] cb Function to continue with after the write
1626 : : // *****************************************************************************
1627 : : {
1628 [ + + ][ + - ]: 914 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1629 : :
1630 [ + - ]: 330 : auto d = Disc();
1631 : :
1632 : : // Field output
1633 : :
1634 : : std::vector< std::string > nodefieldnames{
1635 [ + - ][ + + ]: 1650 : "pressure", "velocityx", "velocityy", "velocityz" };
[ - - ]
1636 : :
1637 : 330 : std::vector< std::vector< tk::real > > nodefields;
1638 : :
1639 [ + - ][ + - ]: 330 : nodefields.push_back( m_u.extract(0) );
1640 [ + - ][ + - ]: 330 : nodefields.push_back( m_u.extract(1) );
1641 [ + - ][ + - ]: 330 : nodefields.push_back( m_u.extract(2) );
1642 [ + - ][ + - ]: 330 : nodefields.push_back( m_u.extract(3) );
1643 : :
1644 : 330 : auto ncomp = m_u.nprop();
1645 [ + + ]: 430 : for (std::size_t c=0; c<ncomp-4; ++c) {
1646 [ + - ][ + - ]: 100 : nodefieldnames.push_back( "c" + std::to_string(c) );
[ + - ]
1647 [ + - ][ + - ]: 100 : nodefields.push_back( m_u.extract(4+c) );
1648 : : }
1649 : :
1650 : : // query function to evaluate analytic solution (if defined)
1651 [ + - ]: 330 : auto sol = problems::SOL();
1652 : :
1653 [ + + ]: 330 : if (sol) {
1654 : 100 : const auto& coord = d->Coord();
1655 : 100 : const auto& x = coord[0];
1656 : 100 : const auto& y = coord[1];
1657 : 100 : const auto& z = coord[2];
1658 [ + - ]: 100 : auto an = m_u;
1659 [ + + ]: 11604 : for (std::size_t i=0; i<an.nunk(); ++i) {
1660 [ + - ]: 11504 : auto s = sol( x[i], y[i], z[i], d->T() );
1661 [ + - ][ + + ]: 69024 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
1662 : 11504 : }
1663 [ + - ][ + - ]: 100 : nodefieldnames.push_back( "pressure_analytic" );
1664 [ + - ][ + - ]: 100 : nodefields.push_back( an.extract(0) );
1665 [ + - ][ + - ]: 100 : nodefieldnames.push_back( "velocity_analyticx" );
1666 [ + - ][ + - ]: 100 : nodefields.push_back( an.extract(1) );
1667 [ + - ][ + - ]: 100 : nodefieldnames.push_back( "velocity_analyticy" );
1668 [ + - ][ + - ]: 100 : nodefields.push_back( an.extract(2) );
1669 [ + - ][ + - ]: 100 : nodefieldnames.push_back( "velocity_analyticz" );
1670 [ + - ][ + - ]: 100 : nodefields.push_back( an.extract(3) );
1671 [ + + ]: 200 : for (std::size_t c=0; c<ncomp-4; ++c) {
1672 [ + - ][ + - ]: 100 : nodefieldnames.push_back( nodefieldnames[4+c] + "_analytic" );
1673 [ + - ][ + - ]: 100 : nodefields.push_back( an.extract(4+c) );
1674 : : }
1675 : 100 : }
1676 : :
1677 [ - + ][ - - ]: 330 : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
1678 : :
1679 : : // Surface output
1680 : :
1681 : 330 : std::vector< std::string > nodesurfnames;
1682 : 330 : std::vector< std::vector< tk::real > > nodesurfs;
1683 : :
1684 : 330 : const auto& f = g_cfg.get< tag::fieldout >();
1685 : :
1686 [ + + ]: 330 : if (!f.empty()) {
1687 : 22 : std::size_t nc = 4;
1688 [ + - ][ + - ]: 22 : nodesurfnames.push_back( "pressure" );
1689 [ + - ][ + - ]: 22 : nodesurfnames.push_back( "velocityx" );
1690 [ + - ][ + - ]: 22 : nodesurfnames.push_back( "velocityy" );
1691 [ + - ][ + - ]: 22 : nodesurfnames.push_back( "velocityz" );
1692 : :
1693 [ + - ]: 22 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
1694 [ + - ]: 22 : std::set< int > outsets( begin(f), end(f) );
1695 [ + + ]: 44 : for (auto sideset : outsets) {
1696 [ + - ]: 22 : auto b = bnode.find(sideset);
1697 [ + + ]: 22 : if (b == end(bnode)) continue;
1698 : 16 : const auto& nodes = b->second;
1699 : 16 : auto i = nodesurfs.size();
1700 [ + - ]: 16 : nodesurfs.insert( end(nodesurfs), nc,
1701 [ + - ]: 32 : std::vector< tk::real >( nodes.size() ) );
1702 : 16 : std::size_t j = 0;
1703 [ + + ]: 1932 : for (auto n : nodes) {
1704 [ + - ]: 1916 : const auto s = m_u[n];
1705 : 1916 : std::size_t p = 0;
1706 : 1916 : nodesurfs[i+(p++)][j] = s[0];
1707 : 1916 : nodesurfs[i+(p++)][j] = s[1];
1708 : 1916 : nodesurfs[i+(p++)][j] = s[2];
1709 : 1916 : nodesurfs[i+(p++)][j] = s[3];
1710 : 1916 : ++j;
1711 : 1916 : }
1712 : : }
1713 : 22 : }
1714 : :
1715 : : // Send mesh and fields data (solution dump) for output to file
1716 [ + - ][ + - ]: 660 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
1717 : 330 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
1718 : : {}, nodefields, {}, nodesurfs, cb );
1719 [ + - ][ + - ]: 990 : }
[ + - ][ + - ]
[ - - ][ - - ]
1720 : :
1721 : : void
1722 : 4730 : LohCG::out()
1723 : : // *****************************************************************************
1724 : : // Output mesh field data
1725 : : // *****************************************************************************
1726 : : {
1727 : 4730 : auto d = Disc();
1728 : :
1729 : : // Time history
1730 [ + - ][ + - ]: 4730 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ][ - + ]
1731 : 0 : auto ncomp = m_u.nprop();
1732 : 0 : const auto& inpoel = d->Inpoel();
1733 [ - - ]: 0 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
1734 : 0 : std::size_t j = 0;
1735 [ - - ]: 0 : for (const auto& p : d->Hist()) {
1736 : 0 : auto e = p.get< tag::elem >(); // host element id
1737 : 0 : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
1738 [ - - ]: 0 : hist[j].resize( ncomp+1, 0.0 );
1739 [ - - ]: 0 : for (std::size_t i=0; i<4; ++i) {
1740 [ - - ]: 0 : const auto u = m_u[ inpoel[e*4+i] ];
1741 : 0 : hist[j][0] += n[i] * u[0];
1742 : 0 : hist[j][1] += n[i] * u[1]/u[0];
1743 : 0 : hist[j][2] += n[i] * u[2]/u[0];
1744 : 0 : hist[j][3] += n[i] * u[3]/u[0];
1745 : 0 : hist[j][4] += n[i] * u[4]/u[0];
1746 : 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];
1747 : 0 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
1748 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
1749 : 0 : }
1750 : 0 : ++j;
1751 : : }
1752 [ - - ]: 0 : d->history( std::move(hist) );
1753 : 0 : }
1754 : :
1755 : : // Field data
1756 [ + + ][ + - ]: 4730 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
[ + + ]
1757 [ + - ][ + - ]: 529 : writeFields( CkCallback(CkIndex_LohCG::integrals(), thisProxy[thisIndex]) );
[ + - ][ + - ]
1758 : : } else {
1759 : 4201 : integrals();
1760 : : }
1761 : 4730 : }
1762 : :
1763 : : void
1764 : 4730 : LohCG::integrals()
1765 : : // *****************************************************************************
1766 : : // Compute integral quantities for output
1767 : : // *****************************************************************************
1768 : : {
1769 : 4730 : auto d = Disc();
1770 : :
1771 [ + + ][ + - ]: 4730 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ][ + + ]
1772 : :
1773 : : using namespace integrals;
1774 [ + - ]: 22 : std::vector< std::map< int, tk::real > > ints( NUMINT );
1775 : :
1776 : : // Prepend integral vector with metadata on the current time step:
1777 : : // current iteration count, current physical time, time step size
1778 [ + - ]: 22 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
1779 [ + - ]: 22 : ints[ TIME ][ 0 ] = d->T();
1780 [ + - ]: 22 : ints[ DT ][ 0 ] = d->Dt();
1781 : :
1782 : : // Compute integrals requested for surfaces requested
1783 : 22 : const auto& reqv = g_cfg.get< tag::integout_integrals >();
1784 [ + - ]: 22 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
1785 [ + - ][ + - ]: 22 : if (req.count("mass_flow_rate")) {
[ - + ]
1786 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
1787 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
1788 : 0 : const auto& nodes = sint.first;
1789 : 0 : const auto& ndA = sint.second;
1790 : 0 : auto n = ndA.data();
1791 [ - - ]: 0 : for (auto p : nodes) {
1792 [ - - ][ - - ]: 0 : mfr += n[0]*m_u(p,1) + n[1]*m_u(p,2) + n[2]*m_u(p,3);
[ - - ]
1793 : 0 : n += 3;
1794 : : }
1795 : : }
1796 : : }
1797 [ + - ][ + - ]: 22 : if (req.count("force")) {
[ + - ]
1798 : 22 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
1799 [ + + ]: 38 : for (const auto& [s,sint] : m_surfint) {
1800 [ + - ]: 16 : auto& fx = ints[ FORCE_X ][ s ];
1801 [ + - ]: 16 : auto& fy = ints[ FORCE_Y ][ s ];
1802 [ + - ]: 16 : auto& fz = ints[ FORCE_Z ][ s ];
1803 : 16 : const auto& nodes = sint.first;
1804 : 16 : const auto& ndA = sint.second;
1805 : 16 : auto n = ndA.data();
1806 [ + + ]: 1932 : for (auto p : nodes) {
1807 : : // pressure force
1808 [ + - ]: 1916 : fx -= n[0]*m_u(p,0);
1809 [ + - ]: 1916 : fy -= n[1]*m_u(p,0);
1810 [ + - ]: 1916 : fz -= n[2]*m_u(p,0);
1811 : : // viscous force
1812 [ + - ][ + - ]: 1916 : fx += mu * (m_grad(p,3)*n[0] + m_grad(p,4)*n[1] + m_grad(p,5)*n[2]);
[ + - ]
1813 [ + - ][ + - ]: 1916 : fy += mu * (m_grad(p,6)*n[0] + m_grad(p,7)*n[1] + m_grad(p,8)*n[2]);
[ + - ]
1814 [ + - ][ + - ]: 1916 : fz += mu * (m_grad(p,9)*n[0] + m_grad(p,10)*n[1] + m_grad(p,11)*n[2]);
[ + - ]
1815 : 1916 : n += 3;
1816 : : }
1817 : : }
1818 : : }
1819 : :
1820 [ + - ]: 22 : auto stream = serialize( d->MeshId(), ints );
1821 [ + - ]: 22 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
1822 [ + - ][ + - ]: 44 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
1823 : :
1824 : 22 : } else {
1825 : :
1826 : 4708 : step();
1827 : :
1828 : : }
1829 : 4730 : }
1830 : :
1831 : : void
1832 : 4345 : LohCG::evalLB( int nrestart )
1833 : : // *****************************************************************************
1834 : : // Evaluate whether to do load balancing
1835 : : //! \param[in] nrestart Number of times restarted
1836 : : // *****************************************************************************
1837 : : {
1838 : 4345 : auto d = Disc();
1839 : :
1840 : : // Detect if just returned from a checkpoint and if so, zero timers and
1841 : : // finished flag
1842 [ + + ]: 4345 : if (d->restarted( nrestart )) m_finished = 0;
1843 : :
1844 : 4345 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
1845 : 4345 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
1846 : :
1847 : : // Load balancing if user frequency is reached or after the second time-step
1848 [ + + ][ + + ]: 4345 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
[ + + ]
1849 : :
1850 : 3525 : AtSync();
1851 [ - + ]: 3525 : if (nonblocking) dt();
1852 : :
1853 : : } else {
1854 : :
1855 : 820 : dt();
1856 : :
1857 : : }
1858 : 4345 : }
1859 : :
1860 : : void
1861 : 4340 : LohCG::evalRestart()
1862 : : // *****************************************************************************
1863 : : // Evaluate whether to save checkpoint/restart
1864 : : // *****************************************************************************
1865 : : {
1866 : 4340 : auto d = Disc();
1867 : :
1868 : 4340 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1869 : 4340 : const auto benchmark = g_cfg.get< tag::benchmark >();
1870 : :
1871 [ + + ][ - + ]: 4340 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
1872 : :
1873 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1874 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1875 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
1876 : :
1877 : 0 : } else {
1878 : :
1879 : 4340 : evalLB( /* nrestart = */ -1 );
1880 : :
1881 : : }
1882 : 4340 : }
1883 : :
1884 : : void
1885 : 4730 : LohCG::step()
1886 : : // *****************************************************************************
1887 : : // Evaluate whether to continue with next time step
1888 : : // *****************************************************************************
1889 : : {
1890 : 4730 : auto d = Disc();
1891 : :
1892 : : // Output one-liner status report to screen
1893 [ + + ]: 4730 : if(thisIndex == 0) d->status();
1894 : :
1895 [ + + ]: 4730 : if (not m_finished) {
1896 : :
1897 : 4340 : evalRestart();
1898 : :
1899 : : } else {
1900 : :
1901 : 390 : auto meshid = d->MeshId();
1902 [ + - ]: 390 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1903 [ + - ][ + - ]: 780 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1904 : :
1905 : : }
1906 : 4730 : }
1907 : :
1908 : : #include "NoWarning/lohcg.def.h"
|
