Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/LaxCG.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 LaxCG: Time-derivative preconditioning for all Ma
10 : : \see Luo, Baum, Lohner, "Extension of Harten-Lax-van Leer Scheme for
11 : : Flows at All Speeds", AIAA Journal, Vol. 43, No. 6, 2005
12 : : \see Weiss & Smith, "Preconditioning Applied to Variable and Constant
13 : : Density Time-Accurate Flows on Unstructured Meshes", AIAA Journal,
14 : : Vol. 33, No. 11, 1995, pp. 2050-2057.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include "XystBuildConfig.hpp"
19 : : #include "LaxCG.hpp"
20 : : #include "Vector.hpp"
21 : : #include "Reader.hpp"
22 : : #include "ContainerUtil.hpp"
23 : : #include "UnsMesh.hpp"
24 : : #include "ExodusIIMeshWriter.hpp"
25 : : #include "InciterConfig.hpp"
26 : : #include "DerivedData.hpp"
27 : : #include "Discretization.hpp"
28 : : #include "DiagReducer.hpp"
29 : : #include "IntegralReducer.hpp"
30 : : #include "Integrals.hpp"
31 : : #include "Refiner.hpp"
32 : : #include "Reorder.hpp"
33 : : #include "Around.hpp"
34 : : #include "Lax.hpp"
35 : : #include "Problems.hpp"
36 : : #include "EOS.hpp"
37 : : #include "BC.hpp"
38 : :
39 : : namespace inciter {
40 : :
41 : : extern ctr::Config g_cfg;
42 : :
43 : : static CkReduction::reducerType IntegralsMerger;
44 : :
45 : : //! Runge-Kutta coefficients
46 : : static const std::array< tk::real, 3 > rkcoef{{ 1.0/3.0, 1.0/2.0, 1.0 }};
47 : :
48 : : } // inciter::
49 : :
50 : : using inciter::g_cfg;
51 : : using inciter::LaxCG;
52 : :
53 : 321 : LaxCG::LaxCG( const CProxy_Discretization& disc,
54 : : const std::map< int, std::vector< std::size_t > >& bface,
55 : : const std::map< int, std::vector< std::size_t > >& bnode,
56 : 321 : const std::vector< std::size_t >& triinpoel ) :
57 [ + - ]: 321 : m_disc( disc ),
58 : 321 : m_nrhs( 0 ),
59 : 321 : m_nnorm( 0 ),
60 : 321 : m_nbpint( 0 ),
61 : 321 : m_nbeint( 0 ),
62 : 321 : m_ndeint( 0 ),
63 : 321 : m_ngrad( 0 ),
64 [ + - ]: 321 : m_bnode( bnode ),
65 [ + - ]: 321 : m_bface( bface ),
66 [ + - ][ + - ]: 321 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
67 [ + - ][ + - ]: 321 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
68 [ + - ]: 321 : m_un( m_u.nunk(), m_u.nprop() ),
69 [ + - ]: 321 : m_rhs( m_u.nunk(), m_u.nprop() ),
70 [ + - ]: 321 : m_grad( m_u.nunk(), m_u.nprop()*3 ),
71 : 321 : m_stage( 0 ),
72 [ + - ]: 321 : m_dtp( m_u.nunk(), 0.0 ),
73 [ + - ]: 321 : m_tp( m_u.nunk(), g_cfg.get< tag::t0 >() ),
74 : 963 : m_finished( 0 )
75 : : // *****************************************************************************
76 : : // Constructor
77 : : //! \param[in] disc Discretization proxy
78 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
79 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
80 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
81 : : // *****************************************************************************
82 : : {
83 : 321 : usesAtSync = true; // enable migration at AtSync
84 : :
85 [ + - ]: 321 : auto d = Disc();
86 : :
87 : : // Create new local ids based on mesh locality
88 : 321 : std::unordered_map< std::size_t, std::size_t > map;
89 : 321 : std::size_t n = 0;
90 : :
91 [ + - ][ + - ]: 321 : auto psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
92 [ + + ]: 71647 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
93 [ + - ][ + + ]: 71326 : if (!map.count(p)) map[p] = n++;
[ + - ]
94 [ + + ]: 900764 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
95 [ + - ][ + + ]: 829438 : if (!map.count(q)) map[q] = n++;
[ + - ]
96 : : }
97 : : }
98 : :
99 [ - + ][ - - ]: 321 : Assert( map.size() == d->Gid().size(),
[ - - ][ - - ]
100 : : "Mesh-locality reorder map size mismatch" );
101 : :
102 : : // Remap data in bound Discretization object
103 [ + - ]: 321 : d->remap( map );
104 : : // Remap boundary triangle face connectivity
105 [ + - ]: 321 : tk::remap( m_triinpoel, map );
106 : :
107 : : // Compute total box IC volume
108 [ + - ]: 321 : d->boxvol();
109 : :
110 : : // Activate SDAG wait for initially computing integrals
111 [ + - ][ + - ]: 321 : thisProxy[ thisIndex ].wait4int();
112 : 321 : }
113 : :
114 : : void
115 : 10350 : LaxCG::primitive( tk::Fields& U )
116 : : // *****************************************************************************
117 : : // Convert from conservative to primitive variables
118 : : //! \param[in,out] U Unknown/solution vector to convert
119 : : //! \details On input U is assumed to contain the conservative variables r, ru,
120 : : //! rv, rw, rE, and on output the primitive variables p, u, v, w, T.
121 : : // *****************************************************************************
122 : : {
123 : 10350 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
124 : :
125 [ + + ]: 3330990 : for (std::size_t i=0; i<U.nunk(); ++i) {
126 : 3320640 : auto r = U(i,0);
127 : 3320640 : auto u = U(i,1)/r;
128 : 3320640 : auto v = U(i,2)/r;
129 : 3320640 : auto w = U(i,3)/r;
130 : 3320640 : auto p = eos::pressure( U(i,4) - 0.5*r*(u*u + v*v + w*w) );
131 : 3320640 : U(i,0) = p;
132 : 3320640 : U(i,1) = u;
133 : 3320640 : U(i,2) = v;
134 : 3320640 : U(i,3) = w;
135 : 3320640 : U(i,4) = p/r/rgas;
136 : : }
137 : 10350 : }
138 : :
139 : : void
140 : 13800 : LaxCG::conservative( tk::Fields& U )
141 : : // *****************************************************************************
142 : : // Convert from primitive to conservative variables
143 : : //! \param[in,out] U Unknown/solution vector to convert
144 : : //! \details On input U is assumed to contain the primitive variables p, u, v,
145 : : //! w, T and on output the conservative variables r, ru, rv, rw, rE.
146 : : // *****************************************************************************
147 : : {
148 : 13800 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
149 : 13800 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
150 : :
151 [ + + ]: 4441320 : for (std::size_t i=0; i<U.nunk(); ++i) {
152 : 4427520 : auto p = U(i,0);
153 : 4427520 : auto u = U(i,1);
154 : 4427520 : auto v = U(i,2);
155 : 4427520 : auto w = U(i,3);
156 : 4427520 : auto T = U(i,4);
157 : 4427520 : auto r = p/T/rgas;
158 : 4427520 : U(i,0) = r;
159 : 4427520 : U(i,1) = r*u;
160 : 4427520 : U(i,2) = r*v;
161 : 4427520 : U(i,3) = r*w;
162 : 4427520 : U(i,4) = p/(g-1.0) + 0.5*r*(u*u + v*v + w*w);
163 : : }
164 : 13800 : }
165 : :
166 : : std::array< tk::real, 5*5 >
167 : 3320640 : LaxCG::precond( const tk::Fields& U, std::size_t i )
168 : : // *****************************************************************************
169 : : // Compute the inverse of the time-derivative preconditioning matrix
170 : : //! \param[in] U Unknown/solution vector to use
171 : : //! \param[in] i Mesh point index
172 : : //! \return Preconditioning matrix inverse
173 : : //! \see Nishikawa, Weiss-Smith Local-Preconditioning Matrix is a Diagonal
174 : : //! Matrix in the Symmetric Form of the Euler Equations, 2021.
175 : : // *****************************************************************************
176 : : {
177 : 3320640 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
178 : 3320640 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
179 : :
180 : 3320640 : auto p = U(i,0);
181 : 3320640 : auto u = U(i,1);
182 : 3320640 : auto v = U(i,2);
183 : 3320640 : auto w = U(i,3);
184 : 3320640 : auto T = U(i,4);
185 : 3320640 : auto r = p/T/rgas;
186 : 3320640 : auto cp = g*rgas/(g-1.0);
187 : 3320640 : auto k = u*u + v*v + w*w;
188 : 3320640 : auto vr = lax::refvel( r, p, std::sqrt(k) );
189 : 3320640 : auto vr2 = vr*vr;
190 : 3320640 : auto rt = -r/T;
191 : 3320640 : auto H = cp*T + k/2.0;
192 : 3320640 : auto theta = 1.0/vr2 - rt/r/cp;
193 : 3320640 : auto coef = r*cp*theta + rt;
194 : :
195 : : return {
196 : 3320640 : (rt*(H - k) + r*cp)/coef,
197 : 3320640 : rt*u/coef,
198 : 3320640 : rt*v/coef,
199 : 3320640 : rt*w/coef,
200 : 3320640 : -rt/coef,
201 : :
202 : 3320640 : -u/r,
203 : 3320640 : 1.0/r,
204 : : 0.0,
205 : : 0.0,
206 : : 0.0,
207 : :
208 : 3320640 : -v/r,
209 : : 0.0,
210 : 3320640 : 1.0/r,
211 : : 0.0,
212 : : 0.0,
213 : :
214 : 3320640 : -w/r,
215 : : 0.0,
216 : : 0.0,
217 : 3320640 : 1.0/r,
218 : : 0.0,
219 : :
220 : 3320640 : -(theta*(H - k) - 1.0)/coef,
221 : 3320640 : -theta*u/coef,
222 : 3320640 : -theta*v/coef,
223 : 3320640 : -theta*w/coef,
224 : 3320640 : theta/coef
225 : 3320640 : };
226 : : }
227 : :
228 : : tk::real
229 : 1029900 : LaxCG::charvel( std::size_t i )
230 : : // *****************************************************************************
231 : : // Compute characteristic velocity of the preconditioned system at a point
232 : : //! \param[in] i Mesh point index
233 : : //! \return Maximum eigenvalue: abs(v') + c'
234 : : // *****************************************************************************
235 : : {
236 : 1029900 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
237 : 1029900 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
238 : 1029900 : auto cp = g*rgas/(g-1.0);
239 : :
240 : 1029900 : auto r = m_u(i,0);
241 : 1029900 : auto u = m_u(i,1) / r;
242 : 1029900 : auto v = m_u(i,2) / r;
243 : 1029900 : auto w = m_u(i,3) / r;
244 : 1029900 : auto k = u*u + v*v + w*w;
245 : 1029900 : auto e = m_u(i,4)/r - k/2.0;
246 : 1029900 : auto p = eos::pressure( r*e );
247 : 1029900 : auto T = p/r/rgas;
248 : 1029900 : auto rp = r/p;
249 : 1029900 : auto rt = -r/T;
250 : 1029900 : auto vel = std::sqrt( k );
251 : 1029900 : auto vr = lax::refvel( r, p, vel );
252 : 1029900 : auto vr2 = vr*vr;
253 : 1029900 : auto beta = rp + rt/r/cp;
254 : 1029900 : auto alpha = 0.5*(1.0 - beta*vr2);
255 : 1029900 : auto vpri = vel*(1.0 - alpha);
256 : 1029900 : auto cpri = std::sqrt( alpha*alpha*k + vr2 );
257 : :
258 : 1029900 : return std::abs(vpri) + cpri;
259 : : }
260 : :
261 : : void
262 : 321 : LaxCG::setupBC()
263 : : // *****************************************************************************
264 : : // Prepare boundary condition data structures
265 : : // *****************************************************************************
266 : : {
267 : : // Query Dirichlet BC nodes associated to side sets
268 : 321 : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
269 [ - + ]: 321 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
270 [ - - ]: 0 : auto k = m_bface.find(s[0]);
271 [ - - ]: 0 : if (k != end(m_bface)) {
272 [ - - ]: 0 : auto& n = dir[ k->first ];
273 [ - - ]: 0 : for (auto f : k->second) {
274 [ - - ]: 0 : n.insert( m_triinpoel[f*3+0] );
275 [ - - ]: 0 : n.insert( m_triinpoel[f*3+1] );
276 [ - - ]: 0 : n.insert( m_triinpoel[f*3+2] );
277 : : }
278 : : }
279 : : }
280 : :
281 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
282 [ + - ]: 321 : const auto& lid = Disc()->Lid();
283 [ - + ]: 321 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
284 [ - - ]: 0 : auto k = m_bnode.find(s[0]);
285 [ - - ]: 0 : if (k != end(m_bnode)) {
286 [ - - ]: 0 : auto& n = dir[ k->first ];
287 [ - - ]: 0 : for (auto g : k->second) {
288 [ - - ][ - - ]: 0 : n.insert( tk::cref_find(lid,g) );
289 : : }
290 : : }
291 : : }
292 : :
293 : : // Collect unique set of nodes + Dirichlet BC components mask
294 : 321 : auto ncomp = m_u.nprop();
295 : 321 : auto nmask = ncomp + 1;
296 : 321 : const auto& dbc = g_cfg.get< tag::bc_dir >();
297 : 321 : std::unordered_map< std::size_t, std::vector< int > > dirbcset;
298 [ - + ]: 321 : for (const auto& mask : dbc) {
299 [ - - ][ - - ]: 0 : ErrChk( mask.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
300 [ - - ]: 0 : auto n = dir.find( mask[0] );
301 [ - - ]: 0 : if (n != end(dir))
302 [ - - ]: 0 : for (auto p : n->second) {
303 [ - - ]: 0 : auto& m = dirbcset[p];
304 [ - - ][ - - ]: 0 : if (m.empty()) m.resize( ncomp, 0 );
305 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
306 [ - - ]: 0 : if (!m[c]) m[c] = mask[c+1]; // overwrite mask if 0 -> 1
307 : : }
308 : : }
309 : : // Compile streamable list of nodes + Dirichlet BC components mask
310 : 321 : tk::destroy( m_dirbcmasks );
311 [ - + ]: 321 : for (const auto& [p,mask] : dirbcset) {
312 [ - - ]: 0 : m_dirbcmasks.push_back( p );
313 [ - - ]: 0 : m_dirbcmasks.insert( end(m_dirbcmasks), begin(mask), end(mask) );
314 : : }
315 [ - + ][ - - ]: 321 : ErrChk( m_dirbcmasks.size() % nmask == 0, "Dirichlet BC masks incomplete" );
[ - - ][ - - ]
316 : :
317 : : // Query pressure BC nodes associated to side sets
318 : 321 : const auto& pbc = g_cfg.get< tag::bc_pre >();
319 : 321 : const auto& pbc_sets = pbc.get< tag::sidesets >();
320 : 321 : std::unordered_map< int, std::unordered_set< std::size_t > > pre;
321 [ - + ]: 321 : for (const auto& ss : pbc_sets) {
322 [ - - ]: 0 : for (const auto& s : ss) {
323 [ - - ]: 0 : auto k = m_bface.find(s);
324 [ - - ]: 0 : if (k != end(m_bface)) {
325 [ - - ]: 0 : auto& n = pre[ k->first ];
326 [ - - ]: 0 : for (auto f : k->second) {
327 [ - - ]: 0 : n.insert( m_triinpoel[f*3+0] );
328 [ - - ]: 0 : n.insert( m_triinpoel[f*3+1] );
329 [ - - ]: 0 : n.insert( m_triinpoel[f*3+2] );
330 : : }
331 : : }
332 : : }
333 : : }
334 : :
335 : : // Prepare density and pressure values for pressure BC nodes
336 [ - + ]: 321 : if (!pbc_sets.empty()) {
337 : 0 : const auto& pbc_r = pbc.get< tag::density >();
338 [ - - ][ - - ]: 0 : ErrChk( pbc_r.size() == pbc_sets.size(), "Pressure BC density unspecified" );
[ - - ][ - - ]
339 : 0 : const auto& pbc_p = pbc.get< tag::pressure >();
340 [ - - ][ - - ]: 0 : ErrChk( pbc_p.size() == pbc_sets.size(), "Pressure BC pressure unspecified" );
[ - - ][ - - ]
341 : 0 : tk::destroy( m_prebcnodes );
342 : 0 : tk::destroy( m_prebcvals );
343 [ - - ]: 0 : for (const auto& [s,n] : pre) {
344 [ - - ]: 0 : m_prebcnodes.insert( end(m_prebcnodes), begin(n), end(n) );
345 [ - - ]: 0 : for (std::size_t p=0; p<pbc_sets.size(); ++p) {
346 [ - - ]: 0 : for (auto u : pbc_sets[p]) {
347 [ - - ]: 0 : if (s == u) {
348 [ - - ]: 0 : for (std::size_t i=0; i<n.size(); ++i) {
349 [ - - ]: 0 : m_prebcvals.push_back( pbc_r[p] );
350 [ - - ]: 0 : m_prebcvals.push_back( pbc_p[p] );
351 : : }
352 : : }
353 : : }
354 : : }
355 : : }
356 [ - - ][ - - ]: 0 : ErrChk( m_prebcnodes.size()*2 == m_prebcvals.size(),
[ - - ][ - - ]
357 : : "Pressure BC data incomplete" );
358 : : }
359 : :
360 : : // Query symmetry BC nodes associated to side sets
361 : 321 : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
362 [ + + ]: 642 : for (auto s : g_cfg.get< tag::bc_sym >()) {
363 [ + - ]: 321 : auto k = m_bface.find(s);
364 [ + + ]: 321 : if (k != end(m_bface)) {
365 [ + - ]: 139 : auto& n = sym[ k->first ];
366 [ + + ]: 32273 : for (auto f : k->second) {
367 [ + - ]: 32134 : n.insert( m_triinpoel[f*3+0] );
368 [ + - ]: 32134 : n.insert( m_triinpoel[f*3+1] );
369 [ + - ]: 32134 : n.insert( m_triinpoel[f*3+2] );
370 : : }
371 : : }
372 : : }
373 : :
374 : : // Query farfield BC nodes associated to side sets
375 : 321 : std::unordered_map< int, std::unordered_set< std::size_t > > far;
376 [ + + ]: 350 : for (auto s : g_cfg.get< tag::bc_far, tag::sidesets >()) {
377 [ + - ]: 29 : auto k = m_bface.find(s);
378 [ + + ]: 29 : if (k != end(m_bface)) {
379 [ + - ]: 10 : auto& n = far[ k->first ];
380 [ + + ]: 1190 : for (auto f : k->second) {
381 [ + - ]: 1180 : n.insert( m_triinpoel[f*3+0] );
382 [ + - ]: 1180 : n.insert( m_triinpoel[f*3+1] );
383 [ + - ]: 1180 : n.insert( m_triinpoel[f*3+2] );
384 : : }
385 : : }
386 : : }
387 : :
388 : : // Generate unique set of symmetry BC nodes
389 : 321 : tk::destroy( m_symbcnodeset );
390 [ + - ][ + + ]: 460 : for (const auto& [s,n] : sym) m_symbcnodeset.insert( begin(n), end(n) );
391 : : // Generate unique set of farfield BC nodes
392 : 321 : tk::destroy( m_farbcnodeset );
393 [ + - ][ + + ]: 331 : for (const auto& [s,n] : far) m_farbcnodeset.insert( begin(n), end(n) );
394 : :
395 : : // If farfield BC is set on a node, will not also set symmetry BC
396 [ + - ][ + + ]: 1125 : for (auto i : m_farbcnodeset) m_symbcnodeset.erase(i);
397 : 321 : }
398 : :
399 : : void
400 : 321 : LaxCG::feop()
401 : : // *****************************************************************************
402 : : // Start (re-)computing finite element domain and boundary operators
403 : : // *****************************************************************************
404 : : {
405 : 321 : auto d = Disc();
406 : :
407 : : // Prepare boundary conditions data structures
408 : 321 : setupBC();
409 : :
410 : : // Compute local contributions to boundary normals and integrals
411 : 321 : bndint();
412 : : // Compute local contributions to domain edge integrals
413 : 321 : domint();
414 : :
415 : : // Send boundary point normal contributions to neighbor chares
416 [ + + ]: 321 : if (d->NodeCommMap().empty()) {
417 : 4 : comnorm_complete();
418 : : } else {
419 [ + + ]: 3907 : for (const auto& [c,nodes] : d->NodeCommMap()) {
420 : 3590 : decltype(m_bnorm) exp;
421 [ + + ]: 27792 : for (auto i : nodes) {
422 [ + + ]: 62926 : for (const auto& [s,b] : m_bnorm) {
423 [ + - ]: 38724 : auto k = b.find(i);
424 [ + + ][ + - ]: 38724 : if (k != end(b)) exp[s][i] = k->second;
[ + - ]
425 : : }
426 : : }
427 [ + - ][ + - ]: 3590 : thisProxy[c].comnorm( exp );
428 : 3590 : }
429 : : }
430 : 321 : ownnorm_complete();
431 : 321 : }
432 : :
433 : : void
434 : 321 : LaxCG::bndint()
435 : : // *****************************************************************************
436 : : //! Compute local contributions to boundary normals and integrals
437 : : // *****************************************************************************
438 : : {
439 [ + - ]: 321 : auto d = Disc();
440 : 321 : const auto& coord = d->Coord();
441 : 321 : const auto& gid = d->Gid();
442 : 321 : const auto& x = coord[0];
443 : 321 : const auto& y = coord[1];
444 : 321 : const auto& z = coord[2];
445 : :
446 : : // Lambda to compute the inverse distance squared between boundary face
447 : : // centroid and boundary point. Here p is the global node id and c is the
448 : : // the boundary face centroid.
449 : 118302 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
450 : 118302 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
451 : 118302 : (c[1] - y[p]) * (c[1] - y[p]) +
452 : 118302 : (c[2] - z[p]) * (c[2] - z[p]) );
453 : 321 : };
454 : :
455 : 321 : tk::destroy( m_bnorm );
456 : 321 : tk::destroy( m_bndpoinint );
457 : :
458 [ + + ]: 970 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
459 [ + + ]: 40083 : for (auto f : faceids) { // for all side set triangles
460 : 39434 : const auto N = m_triinpoel.data() + f*3;
461 : : const std::array< tk::real, 3 >
462 : 39434 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
463 : 39434 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
464 : 39434 : auto n = tk::cross( ba, ca );
465 : 39434 : auto A2 = tk::length( n );
466 : 39434 : n[0] /= A2;
467 : 39434 : n[1] /= A2;
468 : 39434 : n[2] /= A2;
469 : : const tk::real centroid[3] = {
470 : 39434 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
471 : 39434 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
472 : 78868 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
473 [ + + ]: 157736 : for (const auto& [i,j] : tk::lpoet) {
474 : 118302 : auto p = N[i];
475 : 118302 : tk::real r = invdistsq( centroid, p );
476 [ + - ]: 118302 : auto& v = m_bnorm[setid]; // associate side set id
477 [ + - ]: 118302 : auto& bpn = v[gid[p]]; // associate global node id of bnd pnt
478 : 118302 : bpn[0] += r * n[0]; // inv.dist.sq-weighted normal
479 : 118302 : bpn[1] += r * n[1];
480 : 118302 : bpn[2] += r * n[2];
481 : 118302 : bpn[3] += r; // inv.dist.sq of node from centroid
482 [ + - ]: 118302 : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
483 : 118302 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
484 : 118302 : b[1] += n[1] * A2 / 6.0;
485 : 118302 : b[2] += n[2] * A2 / 6.0;
486 : : }
487 : : }
488 : : }
489 : 321 : }
490 : :
491 : : void
492 : 321 : LaxCG::domint()
493 : : // *****************************************************************************
494 : : //! Compute local contributions to domain edge integrals
495 : : // *****************************************************************************
496 : : {
497 : 321 : auto d = Disc();
498 : :
499 : 321 : const auto& gid = d->Gid();
500 : 321 : const auto& inpoel = d->Inpoel();
501 : :
502 : 321 : const auto& coord = d->Coord();
503 : 321 : const auto& x = coord[0];
504 : 321 : const auto& y = coord[1];
505 : 321 : const auto& z = coord[2];
506 : :
507 : 321 : tk::destroy( m_domedgeint );
508 : :
509 [ + + ]: 311891 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
510 : 311570 : const auto N = inpoel.data() + e*4;
511 : : const std::array< tk::real, 3 >
512 : 311570 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
513 : 311570 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
514 : 311570 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
515 : : std::array< std::array< tk::real, 3 >, 4 > grad;
516 : 311570 : grad[1] = tk::cross( ca, da );
517 : 311570 : grad[2] = tk::cross( da, ba );
518 : 311570 : grad[3] = tk::cross( ba, ca );
519 [ + + ]: 1246280 : for (std::size_t i=0; i<3; ++i)
520 : 934710 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
521 [ + + ]: 2180990 : for (const auto& [p,q] : tk::lpoed) {
522 : 1869420 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
523 : 1869420 : tk::real sig = 1.0;
524 [ + + ]: 1869420 : if (ed[0] > ed[1]) {
525 : 919156 : std::swap( ed[0], ed[1] );
526 : 919156 : sig = -1.0;
527 : : }
528 [ + - ]: 1869420 : auto& n = m_domedgeint[ ed ];
529 : 1869420 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
530 : 1869420 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
531 : 1869420 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
532 : : }
533 : : }
534 : 321 : }
535 : :
536 : : void
537 : 3590 : LaxCG::comnorm( const decltype(m_bnorm)& inbnd )
538 : : // *****************************************************************************
539 : : // Receive contributions to boundary point normals on chare-boundaries
540 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
541 : : // *****************************************************************************
542 : : {
543 : : // Buffer up incoming boundary point normal vector contributions
544 [ + + ]: 6150 : for (const auto& [s,b] : inbnd) {
545 [ + - ]: 2560 : auto& bndnorm = m_bnormc[s];
546 [ + + ]: 10921 : for (const auto& [p,n] : b) {
547 [ + - ]: 8361 : auto& norm = bndnorm[p];
548 : 8361 : norm[0] += n[0];
549 : 8361 : norm[1] += n[1];
550 : 8361 : norm[2] += n[2];
551 : 8361 : norm[3] += n[3];
552 : : }
553 : : }
554 : :
555 [ + + ]: 3590 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
556 : 317 : m_nnorm = 0;
557 : 317 : comnorm_complete();
558 : : }
559 : 3590 : }
560 : :
561 : : void
562 : 804 : LaxCG::registerReducers()
563 : : // *****************************************************************************
564 : : // Configure Charm++ reduction types initiated from this chare array
565 : : //! \details Since this is a [initnode] routine, the runtime system executes the
566 : : //! routine exactly once on every logical node early on in the Charm++ init
567 : : //! sequence. Must be static as it is called without an object. See also:
568 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
569 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
570 : : // *****************************************************************************
571 : : {
572 : 804 : NodeDiagnostics::registerReducers();
573 : 804 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
574 : 804 : }
575 : :
576 : : void
577 : : // cppcheck-suppress unusedFunction
578 : 2844 : LaxCG::ResumeFromSync()
579 : : // *****************************************************************************
580 : : // Return from migration
581 : : //! \details This is called when load balancing (LB) completes. The presence of
582 : : //! this function does not affect whether or not we block on LB.
583 : : // *****************************************************************************
584 : : {
585 [ - + ][ - - ]: 2844 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
586 : :
587 [ + - ]: 2844 : if (!g_cfg.get< tag::nonblocking >()) dt();
588 : 2844 : }
589 : :
590 : : void
591 : 321 : LaxCG::setup( tk::real v )
592 : : // *****************************************************************************
593 : : // Start setup for solution
594 : : //! \param[in] v Total volume within user-specified box
595 : : // *****************************************************************************
596 : : {
597 : 321 : auto d = Disc();
598 : :
599 : : // Store user-defined box IC volume
600 : 321 : Disc()->Boxvol() = v;
601 : :
602 : : // Set initial conditions
603 : 321 : problems::initialize( d->Coord(), m_u, d->T(), /*meshid=*/0, d->BoxNodes() );
604 : :
605 : : // Query time history field output labels from all PDEs integrated
606 [ - + ]: 321 : if (!g_cfg.get< tag::histout, tag::points >().empty()) {
607 : : std::vector< std::string > var
608 [ - - ][ - - ]: 0 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ]
609 : 0 : auto ncomp = m_u.nprop();
610 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c)
611 [ - - ][ - - ]: 0 : var.push_back( "c" + std::to_string(c-5) );
[ - - ]
612 [ - - ]: 0 : d->histheader( std::move(var) );
613 : 0 : }
614 : :
615 : : // Compute finite element operators
616 : 321 : feop();
617 [ - - ][ - - ]: 321 : }
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
618 : :
619 : : void
620 : 321 : LaxCG::start()
621 : : // *****************************************************************************
622 : : // Start time stepping
623 : : // *****************************************************************************
624 : : {
625 : : // Set flag that indicates that we are now during time stepping
626 : 321 : Disc()->Initial( 0 );
627 : : // Start timer measuring time stepping wall clock time
628 : 321 : Disc()->Timer().zero();
629 : : // Zero grind-timer
630 : 321 : Disc()->grindZero();
631 : : // Continue to first time step
632 : 321 : dt();
633 : 321 : }
634 : :
635 : : void
636 : 321 : LaxCG::bnorm()
637 : : // *****************************************************************************
638 : : // Combine own and communicated portions of the boundary point normals
639 : : // *****************************************************************************
640 : : {
641 [ + - ]: 321 : const auto& lid = Disc()->Lid();
642 : :
643 : : // Combine own and communicated contributions to boundary point normals
644 [ + + ]: 998 : for (const auto& [s,b] : m_bnormc) {
645 [ + - ]: 677 : auto& bndnorm = m_bnorm[s];
646 [ + + ]: 6906 : for (const auto& [g,n] : b) {
647 [ + - ]: 6229 : auto& norm = bndnorm[g];
648 : 6229 : norm[0] += n[0];
649 : 6229 : norm[1] += n[1];
650 : 6229 : norm[2] += n[2];
651 : 6229 : norm[3] += n[3];
652 : : }
653 : : }
654 : 321 : tk::destroy( m_bnormc );
655 : :
656 : : // Divide summed point normals by the sum of the inverse distance squared
657 [ + + ]: 1024 : for (auto& [s,b] : m_bnorm) {
658 [ + + ]: 25488 : for (auto& [g,n] : b) {
659 : 24785 : n[0] /= n[3];
660 : 24785 : n[1] /= n[3];
661 : 24785 : n[2] /= n[3];
662 [ - + ][ - - ]: 24785 : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
[ - - ][ - - ]
663 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
664 : : "Non-unit normal" );
665 : : }
666 : : }
667 : :
668 : : // Replace global->local ids associated to boundary point normals
669 : 321 : decltype(m_bnorm) loc;
670 [ + + ]: 1024 : for (auto& [s,b] : m_bnorm) {
671 [ + - ]: 703 : auto& bnd = loc[s];
672 [ + + ]: 25488 : for (auto&& [g,n] : b) {
673 [ + - ][ + - ]: 24785 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
674 : : }
675 : : }
676 : 321 : m_bnorm = std::move(loc);
677 : 321 : }
678 : :
679 : : void
680 : 321 : LaxCG::streamable()
681 : : // *****************************************************************************
682 : : // Convert integrals into streamable data structures
683 : : // *****************************************************************************
684 : : {
685 : : // Generate boundary element symmetry BC flags
686 [ + - ]: 321 : m_besym.resize( m_triinpoel.size() );
687 : 321 : std::size_t i = 0;
688 [ + + ]: 118623 : for (auto p : m_triinpoel) {
689 [ + - ]: 118302 : m_besym[i++] = static_cast< std::uint8_t >(m_symbcnodeset.count(p));
690 : : }
691 : :
692 : : // Query surface integral output nodes
693 : 321 : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
694 : 321 : const auto& is = g_cfg.get< tag::integout, tag::sidesets >();
695 [ + - ]: 321 : std::set< int > outsets( begin(is), end(is) );
696 [ - + ]: 321 : for (auto s : outsets) {
697 [ - - ]: 0 : auto m = m_bface.find(s);
698 [ - - ]: 0 : if (m != end(m_bface)) {
699 [ - - ]: 0 : auto& n = surfintnodes[ m->first ]; // associate set id
700 [ - - ]: 0 : for (auto f : m->second) { // face ids on side set
701 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+0] ); // nodes on side set
702 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+1] );
703 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+2] );
704 : : }
705 : : }
706 : : }
707 [ - - ][ - + ]: 321 : for (auto& [s,n] : surfintnodes) tk::unique( n );
708 : : // Prepare surface integral data
709 : 321 : tk::destroy( m_surfint );
710 [ + - ]: 321 : const auto& gid = Disc()->Gid();
711 [ - + ]: 321 : for (auto&& [s,n] : surfintnodes) {
712 [ - - ]: 0 : auto& sint = m_surfint[s]; // associate set id
713 : 0 : auto& nodes = sint.first;
714 : 0 : auto& ndA = sint.second;
715 : 0 : nodes = std::move(n);
716 [ - - ]: 0 : ndA.resize( nodes.size()*3 );
717 : 0 : std::size_t a = 0;
718 [ - - ]: 0 : for (auto p : nodes) {
719 [ - - ]: 0 : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
720 : 0 : ndA[a*3+0] = b[0]; // store ni * dA
721 : 0 : ndA[a*3+1] = b[1];
722 : 0 : ndA[a*3+2] = b[2];
723 : 0 : ++a;
724 : : }
725 : : }
726 : 321 : tk::destroy( m_bndpoinint );
727 : :
728 : : // Generate domain superedges
729 [ + - ]: 321 : domsuped();
730 : 321 : tk::destroy( m_domedgeint );
731 : :
732 : : // Convert symmetry BC data to streamable data structures
733 : 321 : tk::destroy( m_symbcnodes );
734 : 321 : tk::destroy( m_symbcnorms );
735 [ + + ]: 17839 : for (auto p : m_symbcnodeset) {
736 [ + + ]: 35036 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
737 [ + - ]: 17518 : auto m = m_bnorm.find(s);
738 [ + - ]: 17518 : if (m != end(m_bnorm)) {
739 [ + - ]: 17518 : auto r = m->second.find(p);
740 [ + - ]: 17518 : if (r != end(m->second)) {
741 [ + - ]: 17518 : m_symbcnodes.push_back( p );
742 [ + - ]: 17518 : m_symbcnorms.push_back( r->second[0] );
743 [ + - ]: 17518 : m_symbcnorms.push_back( r->second[1] );
744 [ + - ]: 17518 : m_symbcnorms.push_back( r->second[2] );
745 : : }
746 : : }
747 : : }
748 : : }
749 : 321 : tk::destroy( m_symbcnodeset );
750 : :
751 : : // Convert farfield BC data to streamable data structures
752 : 321 : tk::destroy( m_farbcnodes );
753 : 321 : tk::destroy( m_farbcnorms );
754 [ + + ]: 1125 : for (auto p : m_farbcnodeset) {
755 [ + + ]: 1608 : for (const auto& s : g_cfg.get< tag::bc_far, tag::sidesets >()) {
756 [ + - ]: 804 : auto n = m_bnorm.find(s);
757 [ + - ]: 804 : if (n != end(m_bnorm)) {
758 [ + - ]: 804 : auto a = n->second.find(p);
759 [ + - ]: 804 : if (a != end(n->second)) {
760 [ + - ]: 804 : m_farbcnodes.push_back( p );
761 [ + - ]: 804 : m_farbcnorms.push_back( a->second[0] );
762 [ + - ]: 804 : m_farbcnorms.push_back( a->second[1] );
763 [ + - ]: 804 : m_farbcnorms.push_back( a->second[2] );
764 : : }
765 : : }
766 : : }
767 : : }
768 : 321 : tk::destroy( m_farbcnodeset );
769 : 321 : tk::destroy( m_bnorm );
770 : 321 : }
771 : :
772 : : void
773 : 321 : LaxCG::domsuped()
774 : : // *****************************************************************************
775 : : // Generate superedge-groups for domain-edge loops
776 : : //! \see See Lohner, Sec. 15.1.6.2, An Introduction to Applied CFD Techniques,
777 : : //! Wiley, 2008.
778 : : // *****************************************************************************
779 : : {
780 [ - + ][ - - ]: 321 : Assert( !m_domedgeint.empty(), "No domain edges to group" );
[ - - ][ - - ]
781 : :
782 : : #ifndef NDEBUG
783 : 321 : auto nedge = m_domedgeint.size();
784 : : #endif
785 : :
786 [ + - ]: 321 : const auto& inpoel = Disc()->Inpoel();
787 [ + - ]: 321 : const auto& lid = Disc()->Lid();
788 [ + - ]: 321 : const auto& gid = Disc()->Gid();
789 : :
790 : 321 : tk::destroy( m_dsupedge[0] );
791 : 321 : tk::destroy( m_dsupedge[1] );
792 : 321 : tk::destroy( m_dsupedge[2] );
793 : :
794 : 321 : tk::destroy( m_dsupint[0] );
795 : 321 : tk::destroy( m_dsupint[1] );
796 : 321 : tk::destroy( m_dsupint[2] );
797 : :
798 : 321 : tk::UnsMesh::FaceSet untri;
799 [ + + ]: 311891 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
800 : : std::size_t N[4] = {
801 : 311570 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
802 [ + - ][ + + ]: 1557850 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
803 : : }
804 : :
805 [ + + ]: 311891 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
806 : : std::size_t N[4] = {
807 : 311570 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
808 : 311570 : int f = 0;
809 : : tk::real sig[6];
810 [ + - ][ + + ]: 2180990 : decltype(m_domedgeint)::const_iterator d[6];
811 [ + + ]: 968794 : for (const auto& [p,q] : tk::lpoed) {
812 : 924520 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
813 [ + + ]: 924520 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
814 [ + - ]: 924520 : d[f] = m_domedgeint.find( ed );
815 [ + + ]: 924520 : if (d[f] == end(m_domedgeint)) break; else ++f;
816 : : }
817 [ + + ]: 311570 : if (f == 6) {
818 [ + - ]: 44274 : m_dsupedge[0].push_back( N[0] );
819 [ + - ]: 44274 : m_dsupedge[0].push_back( N[1] );
820 [ + - ]: 44274 : m_dsupedge[0].push_back( N[2] );
821 [ + - ]: 44274 : m_dsupedge[0].push_back( N[3] );
822 [ + - ][ + + ]: 221370 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
823 [ + + ]: 309918 : for (int ed=0; ed<6; ++ed) {
824 [ + - ]: 265644 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[0] );
825 [ + - ]: 265644 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[1] );
826 [ + - ]: 265644 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[2] );
827 [ + - ]: 265644 : m_domedgeint.erase( d[ed] );
828 : : }
829 : : }
830 : : }
831 : :
832 [ + + ]: 478509 : for (const auto& N : untri) {
833 : 478188 : int f = 0;
834 : : tk::real sig[3];
835 [ + - ][ + + ]: 1912752 : decltype(m_domedgeint)::const_iterator d[3];
836 [ + + ]: 731990 : for (const auto& [p,q] : tk::lpoet) {
837 : 713155 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
838 [ + + ]: 713155 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
839 [ + - ]: 713155 : d[f] = m_domedgeint.find( ed );
840 [ + + ]: 713155 : if (d[f] == end(m_domedgeint)) break; else ++f;
841 : : }
842 [ + + ]: 478188 : if (f == 3) {
843 [ + - ]: 18835 : m_dsupedge[1].push_back( N[0] );
844 [ + - ]: 18835 : m_dsupedge[1].push_back( N[1] );
845 [ + - ]: 18835 : m_dsupedge[1].push_back( N[2] );
846 [ + + ]: 75340 : for (int ed=0; ed<3; ++ed) {
847 [ + - ]: 56505 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[0] );
848 [ + - ]: 56505 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[1] );
849 [ + - ]: 56505 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[2] );
850 [ + - ]: 56505 : m_domedgeint.erase( d[ed] );
851 : : }
852 : : }
853 : : }
854 : :
855 [ + - ]: 321 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
856 [ + - ]: 321 : m_dsupint[2].resize( m_domedgeint.size()*3 );
857 : 321 : std::size_t k = 0;
858 [ + + ]: 92891 : for (const auto& [ed,d] : m_domedgeint) {
859 : 92570 : auto e = m_dsupedge[2].data() + k*2;
860 [ + - ]: 92570 : e[0] = tk::cref_find( lid, ed[0] );
861 [ + - ]: 92570 : e[1] = tk::cref_find( lid, ed[1] );
862 : 92570 : auto i = m_dsupint[2].data() + k*3;
863 : 92570 : i[0] = d[0];
864 : 92570 : i[1] = d[1];
865 : 92570 : i[2] = d[2];
866 : 92570 : ++k;
867 : : }
868 : :
869 : : //std::cout << std::setprecision(2)
870 : : // << "superedges: ntet:" << m_dsupedge[0].size()/4 << "(nedge:"
871 : : // << m_dsupedge[0].size()/4*6 << ","
872 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[0].size()/4*6 ) /
873 : : // static_cast< tk::real >( nedge )
874 : : // << "%) + ntri:" << m_dsupedge[1].size()/3
875 : : // << "(nedge:" << m_dsupedge[1].size() << ","
876 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[1].size() ) /
877 : : // static_cast< tk::real >( nedge )
878 : : // << "%) + nedge:"
879 : : // << m_dsupedge[2].size()/2 << "("
880 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[2].size()/2 ) /
881 : : // static_cast< tk::real >( nedge )
882 : : // << "%) = " << m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
883 : : // m_dsupedge[2].size()/2 << " of "<< nedge << " total edges\n";
884 : :
885 [ - + ][ - - ]: 321 : Assert( m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
[ - - ][ - - ]
886 : : m_dsupedge[2].size()/2 == nedge,
887 : : "Not all edges accounted for in superedge groups" );
888 : 321 : }
889 : :
890 : : void
891 : : // cppcheck-suppress unusedFunction
892 : 321 : LaxCG::merge()
893 : : // *****************************************************************************
894 : : // Combine own and communicated portions of the integrals
895 : : // *****************************************************************************
896 : : {
897 : 321 : auto d = Disc();
898 : :
899 : : // Combine own and communicated contributions to boundary point normals
900 : 321 : bnorm();
901 : :
902 : : // Convert integrals into streamable data structures
903 : 321 : streamable();
904 : :
905 : : // Enforce boundary conditions using (re-)computed boundary data
906 : 321 : BC( d->T() );
907 : :
908 [ + - ]: 321 : if (d->Initial()) {
909 : : // Output initial conditions to file
910 [ + - ][ + - ]: 321 : writeFields( CkCallback(CkIndex_LaxCG::start(), thisProxy[thisIndex]) );
[ + - ][ + - ]
911 : : } else {
912 : 0 : feop_complete();
913 : : }
914 : 321 : }
915 : :
916 : : void
917 : 10671 : LaxCG::BC( tk::real t )
918 : : // *****************************************************************************
919 : : // Apply boundary conditions
920 : : //! \param[in] t Physical time
921 : : // *****************************************************************************
922 : : {
923 : 10671 : auto d = Disc();
924 : :
925 : : // Apply Dirichlet BCs
926 [ + - ]: 21342 : physics::dirbc( d->MeshId(), m_u, t, d->Coord(), d->BoxNodes(),
927 : 10671 : m_dirbcmasks );
928 : :
929 : : // Apply symmetry BCs
930 : 10671 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
931 : :
932 : : // Apply farfield BCs
933 : 10671 : physics::farbc( m_u, m_farbcnodes, m_farbcnorms );
934 : :
935 : : // Apply pressure BCs
936 : 10671 : physics::prebc( m_u, m_prebcnodes, m_prebcvals );
937 : 10671 : }
938 : :
939 : : void
940 : 3450 : LaxCG::dt()
941 : : // *****************************************************************************
942 : : // Compute time step size
943 : : // *****************************************************************************
944 : : {
945 : 3450 : tk::real mindt = std::numeric_limits< tk::real >::max();
946 : :
947 : 3450 : auto const_dt = g_cfg.get< tag::dt >();
948 : 3450 : auto eps = std::numeric_limits< tk::real >::epsilon();
949 [ + - ]: 3450 : auto d = Disc();
950 : :
951 : : // use constant dt if configured
952 [ + + ]: 3450 : if (std::abs(const_dt) > eps) {
953 : :
954 : : // cppcheck-suppress redundantInitialization
955 : 2920 : mindt = const_dt;
956 : :
957 : : } else {
958 : :
959 : 530 : const auto& vol = d->Vol();
960 : 530 : auto cfl = g_cfg.get< tag::cfl >();
961 : :
962 [ + - ]: 530 : if (g_cfg.get< tag::steady >()) {
963 : :
964 [ + + ]: 1030430 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
965 [ + - ]: 1029900 : auto v = charvel( i );
966 : 1029900 : auto L = std::cbrt( vol[i] );
967 : 1029900 : m_dtp[i] = L / std::max( v, 1.0e-8 ) * cfl;
968 : : }
969 [ + - ]: 530 : mindt = *std::min_element( begin(m_dtp), end(m_dtp) );
970 : :
971 : : } else {
972 : :
973 [ - - ]: 0 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
974 [ - - ]: 0 : auto v = charvel( i );
975 : 0 : auto L = std::cbrt( vol[i] );
976 : 0 : auto euler_dt = L / std::max( v, 1.0e-8 );
977 : 0 : mindt = std::min( mindt, euler_dt );
978 : : }
979 : 0 : mindt *= cfl;
980 : :
981 : : }
982 : :
983 : : }
984 : :
985 : : // Actiavate SDAG waits for next time step stage
986 [ + - ][ + - ]: 3450 : thisProxy[ thisIndex ].wait4grad();
987 [ + - ][ + - ]: 3450 : thisProxy[ thisIndex ].wait4rhs();
988 : :
989 : : // Contribute to minimum dt across all chares and advance to next step
990 [ + - ]: 3450 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
991 [ + - ][ + - ]: 6900 : CkCallback(CkReductionTarget(LaxCG,advance), thisProxy) );
992 : 3450 : }
993 : :
994 : : void
995 : 3450 : LaxCG::advance( tk::real newdt )
996 : : // *****************************************************************************
997 : : // Advance equations to next time step
998 : : //! \param[in] newdt The smallest dt across the whole problem
999 : : // *****************************************************************************
1000 : : {
1001 : : // Detect blowup
1002 : 3450 : auto eps = std::numeric_limits< tk::real >::epsilon();
1003 [ - + ]: 3450 : if (newdt < eps) m_finished = 1;
1004 : :
1005 : : // Set new time step size
1006 [ + - ]: 3450 : if (m_stage == 0) Disc()->setdt( newdt );
1007 : :
1008 : 3450 : grad();
1009 : 3450 : }
1010 : :
1011 : : void
1012 : 10350 : LaxCG::grad()
1013 : : // *****************************************************************************
1014 : : // Compute gradients for next time step
1015 : : // *****************************************************************************
1016 : : {
1017 : 10350 : auto d = Disc();
1018 : :
1019 : : // Convert unknowns: r,ru,rv,rw,rE -> p,u,v,w,T
1020 : 10350 : primitive( m_u );
1021 : :
1022 : 10350 : lax::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_grad );
1023 : :
1024 : : // Send gradient contributions to neighbor chares
1025 [ + + ]: 10350 : if (d->NodeCommMap().empty()) {
1026 : 150 : comgrad_complete();
1027 : : } else {
1028 : 10200 : const auto& lid = d->Lid();
1029 [ + + ]: 120780 : for (const auto& [c,n] : d->NodeCommMap()) {
1030 : 110580 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1031 [ + - ][ + - ]: 1116060 : for (auto g : n) exp[g] = m_grad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1032 [ + - ][ + - ]: 110580 : thisProxy[c].comgrad( exp );
1033 : 110580 : }
1034 : : }
1035 : 10350 : owngrad_complete();
1036 : 10350 : }
1037 : :
1038 : : void
1039 : 110580 : LaxCG::comgrad(
1040 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
1041 : : // *****************************************************************************
1042 : : // Receive contributions to node gradients on chare-boundaries
1043 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
1044 : : //! global mesh node IDs, value: contributions for all scalar components.
1045 : : //! \details This function receives contributions to m_grad, which stores the
1046 : : //! gradients at mesh nodes. While m_grad stores own contributions, m_gradc
1047 : : //! collects the neighbor chare contributions during communication. This way
1048 : : //! work on m_grad and m_gradc is overlapped. The two are combined in rhs().
1049 : : // *****************************************************************************
1050 : : {
1051 : : using tk::operator+=;
1052 [ + - ][ + - ]: 1116060 : for (const auto& [g,r] : ingrad) m_gradc[g] += r;
[ + + ]
1053 : :
1054 : : // When we have heard from all chares we communicate with, this chare is done
1055 [ + + ]: 110580 : if (++m_ngrad == Disc()->NodeCommMap().size()) {
1056 : 10200 : m_ngrad = 0;
1057 : 10200 : comgrad_complete();
1058 : : }
1059 : 110580 : }
1060 : :
1061 : : void
1062 : 10350 : LaxCG::rhs()
1063 : : // *****************************************************************************
1064 : : // Compute right-hand side of transport equations
1065 : : // *****************************************************************************
1066 : : {
1067 : 10350 : auto d = Disc();
1068 : 10350 : const auto& lid = d->Lid();
1069 : 10350 : const auto steady = g_cfg.get< tag::steady >();
1070 : :
1071 : : // Combine own and communicated contributions to gradients
1072 [ + + ]: 720630 : for (const auto& [g,r] : m_gradc) {
1073 [ + - ]: 710280 : auto i = tk::cref_find( lid, g );
1074 [ + - ][ + + ]: 11364480 : for (std::size_t c=0; c<r.size(); ++c) m_grad(i,c) += r[c];
1075 : : }
1076 : 10350 : tk::destroy(m_gradc);
1077 : :
1078 : : // divide weak result in gradients by nodal volume
1079 : 10350 : const auto& vol = d->Vol();
1080 [ + + ]: 3330990 : for (std::size_t p=0; p<m_grad.nunk(); ++p)
1081 [ + + ]: 53130240 : for (std::size_t c=0; c<m_grad.nprop(); ++c)
1082 : 49809600 : m_grad(p,c) /= vol[p];
1083 : :
1084 : : // Compute own portion of right-hand side for all equations
1085 [ + + ]: 10350 : auto prev_rkcoef = m_stage == 0 ? 0.0 : rkcoef[m_stage-1];
1086 : :
1087 [ + + ]: 10350 : if (steady) {
1088 [ + + ]: 3091290 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] += prev_rkcoef * m_dtp[p];
1089 : : }
1090 : :
1091 : 20700 : lax::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_besym, m_grad,
1092 : 20700 : m_u, d->V(), d->T(), m_tp, m_rhs );
1093 : :
1094 [ + + ]: 10350 : if (steady) {
1095 [ + + ]: 3091290 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] -= prev_rkcoef * m_dtp[p];
1096 : : }
1097 : :
1098 : : // Communicate rhs to other chares on chare-boundary
1099 [ + + ]: 10350 : if (d->NodeCommMap().empty()) {
1100 : 150 : comrhs_complete();
1101 : : } else {
1102 [ + + ]: 120780 : for (const auto& [c,n] : d->NodeCommMap()) {
1103 : 110580 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1104 [ + - ][ + - ]: 1116060 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1105 [ + - ][ + - ]: 110580 : thisProxy[c].comrhs( exp );
1106 : 110580 : }
1107 : : }
1108 : 10350 : ownrhs_complete();
1109 : 10350 : }
1110 : :
1111 : : void
1112 : 110580 : LaxCG::comrhs(
1113 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1114 : : // *****************************************************************************
1115 : : // Receive contributions to right-hand side vector on chare-boundaries
1116 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1117 : : //! global mesh node IDs, value: contributions for all scalar components.
1118 : : //! \details This function receives contributions to m_rhs, which stores the
1119 : : //! right hand side vector at mesh nodes. While m_rhs stores own
1120 : : //! contributions, m_rhsc collects the neighbor chare contributions during
1121 : : //! communication. This way work on m_rhs and m_rhsc is overlapped. The two
1122 : : //! are combined in solve().
1123 : : // *****************************************************************************
1124 : : {
1125 : : using tk::operator+=;
1126 [ + - ][ + - ]: 1116060 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
[ + + ]
1127 : :
1128 : : // When we have heard from all chares we communicate with, this chare is done
1129 [ + + ]: 110580 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
1130 : 10200 : m_nrhs = 0;
1131 : 10200 : comrhs_complete();
1132 : : }
1133 : 110580 : }
1134 : :
1135 : : void
1136 : : // cppcheck-suppress unusedFunction
1137 : 10350 : LaxCG::solve()
1138 : : // *****************************************************************************
1139 : : // Advance systems of equations
1140 : : // *****************************************************************************
1141 : : {
1142 [ + - ]: 10350 : auto d = Disc();
1143 [ + - ]: 10350 : const auto lid = d->Lid();
1144 : 10350 : const auto steady = g_cfg.get< tag::steady >();
1145 : :
1146 : : // Combine own and communicated contributions to rhs
1147 [ + + ]: 720630 : for (const auto& [g,r] : m_rhsc) {
1148 [ + - ]: 710280 : auto i = tk::cref_find( lid, g );
1149 [ + - ][ + + ]: 4261680 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1150 : : }
1151 : 10350 : tk::destroy(m_rhsc);
1152 : :
1153 : : // Update state at time n
1154 [ + + ][ + - ]: 10350 : if (m_stage == 0) m_un = m_u;
1155 : :
1156 : : // Advance solution
1157 : 10350 : auto dt = d->Dt();
1158 : 10350 : const auto& vol = d->Vol();
1159 : 10350 : auto ncomp = m_u.nprop();
1160 [ + + ]: 3330990 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1161 [ + + ]: 3320640 : if (steady) dt = m_dtp[i];
1162 : 3320640 : auto R = -rkcoef[m_stage] * dt / vol[i];
1163 : : // flow
1164 [ + - ]: 3320640 : auto P = precond( m_u, i );
1165 [ + - ]: 3320640 : tk::real r[] = { R*m_rhs(i,0), R*m_rhs(i,1), R*m_rhs(i,2),
1166 [ + - ][ + - ]: 3320640 : R*m_rhs(i,3), R*m_rhs(i,4) };
[ + - ][ + - ]
1167 : 3320640 : auto p = P.data();
1168 [ + + ]: 19923840 : for (std::size_t c=0; c<5; ++c, p+=5) {
1169 [ + - ]: 33206400 : m_u(i,c) = m_un(i,c)
1170 [ + - ]: 16603200 : + p[0]*r[0] + p[1]*r[1] + p[2]*r[2] + p[3]*r[3] + p[4]*r[4];
1171 : : }
1172 : : // scalar
1173 [ - - ][ - - ]: 3320640 : for (std::size_t c=5; c<ncomp; ++c) m_u(i,c) = m_un(i,c) + R*m_rhs(i,c);
[ - - ][ - + ]
1174 : : }
1175 : :
1176 : : // Convert unknowns: p,u,v,w,T -> r,ru,rv,rw,rE
1177 [ + - ]: 10350 : conservative( m_u );
1178 : :
1179 : : // Configure and apply scalar source to solution (if defined)
1180 [ + - ]: 10350 : auto src = problems::PHYS_SRC();
1181 [ - + ][ - - ]: 10350 : if (src) src( d->Coord(), d->T(), m_u );
1182 : :
1183 : : // Enforce boundary conditions
1184 [ + - ]: 10350 : BC( d->T() + rkcoef[m_stage] * d->Dt() );
1185 : :
1186 [ + + ]: 10350 : if (m_stage < 2) {
1187 : :
1188 : : // Activate SDAG wait for next time step stage
1189 [ + - ][ + - ]: 6900 : thisProxy[ thisIndex ].wait4grad();
1190 [ + - ][ + - ]: 6900 : thisProxy[ thisIndex ].wait4rhs();
1191 : :
1192 : : // start next time step stage
1193 [ + - ]: 6900 : stage();
1194 : :
1195 : : } else {
1196 : :
1197 : : // Activate SDAG waits for finishing this time step stage
1198 [ + - ][ + - ]: 3450 : thisProxy[ thisIndex ].wait4stage();
1199 : : // Compute diagnostics, e.g., residuals
1200 [ + - ]: 3450 : conservative( m_un );
1201 : 3450 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1202 [ + - ]: 3450 : auto diag = m_diag.rhocompute( *d, m_u, m_un, diag_iter );
1203 : : // Increase number of iterations and physical time
1204 [ + - ]: 3450 : d->next();
1205 : : // Advance physical time for local time stepping
1206 [ + + ]: 3450 : if (steady) {
1207 : : using tk::operator+=;
1208 [ + - ]: 530 : m_tp += m_dtp;
1209 : : }
1210 : : // Evaluate residuals
1211 [ - + ][ - - ]: 3450 : if (!diag) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ - - ]
1212 : :
1213 : : }
1214 : 10350 : }
1215 : :
1216 : : void
1217 : 3450 : LaxCG::evalres( const std::vector< tk::real >& l2res )
1218 : : // *****************************************************************************
1219 : : // Evaluate residuals
1220 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1221 : : //! computed across the whole problem
1222 : : // *****************************************************************************
1223 : : {
1224 [ + + ]: 3450 : if (g_cfg.get< tag::steady >()) {
1225 : 530 : const auto rc = g_cfg.get< tag::rescomp >() - 1;
1226 : 530 : Disc()->residual( l2res[rc] );
1227 : : }
1228 : :
1229 : 3450 : refine();
1230 : 3450 : }
1231 : :
1232 : : void
1233 : 3450 : LaxCG::refine()
1234 : : // *****************************************************************************
1235 : : // Optionally refine/derefine mesh
1236 : : // *****************************************************************************
1237 : : {
1238 : 3450 : auto d = Disc();
1239 : :
1240 : : // See if this is the last time step
1241 [ + + ]: 3450 : if (d->finished()) m_finished = 1;
1242 : :
1243 : 3450 : const auto& ht = g_cfg.get< tag::href >();
1244 : 3450 : auto dtref = ht.get< tag::dt >();
1245 : 3450 : auto dtfreq = ht.get< tag::dtfreq >();
1246 : :
1247 : : // if t>0 refinement enabled and we hit the frequency
1248 [ - + ][ - - ]: 3450 : if (dtref && !(d->It() % dtfreq)) { // refine
[ - + ]
1249 : :
1250 : 0 : d->refined() = 1;
1251 : 0 : d->startvol();
1252 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1253 : :
1254 : : // Activate SDAG waits for re-computing the integrals
1255 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1256 : :
1257 : : } else { // do not refine
1258 : :
1259 : 3450 : d->refined() = 0;
1260 : 3450 : feop_complete();
1261 : 3450 : resize_complete();
1262 : :
1263 : : }
1264 : 3450 : }
1265 : :
1266 : : void
1267 : 0 : LaxCG::resizePostAMR(
1268 : : const std::vector< std::size_t >& /*ginpoel*/,
1269 : : const tk::UnsMesh::Chunk& chunk,
1270 : : const tk::UnsMesh::Coords& coord,
1271 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1272 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1273 : : const std::set< std::size_t >& removedNodes,
1274 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1275 : : nodeCommMap,
1276 : : const std::map< int, std::vector< std::size_t > >& bface,
1277 : : const std::map< int, std::vector< std::size_t > >& bnode,
1278 : : const std::vector< std::size_t >& triinpoel )
1279 : : // *****************************************************************************
1280 : : // Receive new mesh from Refiner
1281 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1282 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1283 : : //! \param[in] coord New mesh node coordinates
1284 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1285 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1286 : : //! \param[in] removedNodes Newly removed mesh node local ids
1287 : : //! \param[in] nodeCommMap New node communication map
1288 : : //! \param[in] bface Boundary-faces mapped to side set ids
1289 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1290 : : //! \param[in] triinpoel Boundary-face connectivity
1291 : : // *****************************************************************************
1292 : : {
1293 [ - - ]: 0 : auto d = Disc();
1294 : :
1295 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1296 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1297 : :
1298 : : // Resize mesh data structures after mesh refinement
1299 [ - - ]: 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1300 : :
1301 [ - - ][ - - ]: 0 : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1302 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1303 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1304 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1305 : :
1306 : : // Remove newly removed nodes from solution vectors
1307 [ - - ]: 0 : m_u.rm( removedNodes );
1308 [ - - ]: 0 : m_un.rm( removedNodes );
1309 [ - - ]: 0 : m_rhs.rm( removedNodes );
1310 [ - - ]: 0 : m_grad.rm( removedNodes );
1311 : :
1312 : : // Resize auxiliary solution vectors
1313 : 0 : auto npoin = coord[0].size();
1314 [ - - ]: 0 : m_u.resize( npoin );
1315 [ - - ]: 0 : m_un.resize( npoin );
1316 [ - - ]: 0 : m_rhs.resize( npoin );
1317 [ - - ]: 0 : m_grad.resize( npoin );
1318 : :
1319 : : // Update solution on new mesh
1320 [ - - ]: 0 : for (const auto& n : addedNodes)
1321 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1322 [ - - ][ - - ]: 0 : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
[ - - ][ - - ]
1323 [ - - ][ - - ]: 0 : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
[ - - ][ - - ]
[ - - ]
1324 : : "Indices of parent-edge nodes out of bounds post-AMR");
1325 [ - - ][ - - ]: 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
[ - - ]
1326 : : }
1327 : :
1328 : : // Update physical-boundary node-, face-, and element lists
1329 [ - - ]: 0 : m_bnode = bnode;
1330 [ - - ]: 0 : m_bface = bface;
1331 [ - - ]: 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1332 : :
1333 : 0 : auto meshid = d->MeshId();
1334 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1335 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1336 : 0 : }
1337 : :
1338 : : void
1339 : 711 : LaxCG::writeFields( CkCallback cb )
1340 : : // *****************************************************************************
1341 : : // Output mesh-based fields to file
1342 : : //! \param[in] cb Function to continue with after the write
1343 : : // *****************************************************************************
1344 : : {
1345 [ + + ][ + - ]: 711 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1346 : :
1347 [ + - ]: 127 : auto d = Disc();
1348 : 127 : auto ncomp = m_u.nprop();
1349 : 127 : auto steady = g_cfg.get< tag::steady >();
1350 : :
1351 : : // Field output
1352 : :
1353 : : std::vector< std::string > nodefieldnames
1354 [ + - ][ + + ]: 889 : {"density", "velocityx", "velocityy", "velocityz", "energy", "pressure"};
[ - - ]
1355 [ + - ][ + - ]: 127 : if (steady) nodefieldnames.push_back( "mach" );
1356 : :
1357 : : using tk::operator/=;
1358 [ + - ]: 127 : auto r = m_u.extract(0);
1359 [ + - ][ + - ]: 127 : auto u = m_u.extract(1); u /= r;
1360 [ + - ][ + - ]: 127 : auto v = m_u.extract(2); v /= r;
1361 [ + - ][ + - ]: 127 : auto w = m_u.extract(3); w /= r;
1362 [ + - ][ + - ]: 127 : auto e = m_u.extract(4); e /= r;
1363 [ + - ]: 127 : std::vector< tk::real > pr( m_u.nunk() ), ma;
1364 [ + - ][ + - ]: 127 : if (steady) ma.resize( m_u.nunk() );
1365 [ + + ]: 385139 : for (std::size_t i=0; i<pr.size(); ++i) {
1366 : 385012 : auto vv = u[i]*u[i] + v[i]*v[i] + w[i]*w[i];
1367 : 385012 : pr[i] = eos::pressure( r[i]*(e[i] - 0.5*vv) );
1368 [ + - ]: 385012 : if (steady) ma[i] = std::sqrt(vv) / eos::soundspeed( r[i], pr[i] );
1369 : : }
1370 : :
1371 : : std::vector< std::vector< tk::real > > nodefields{
1372 : 635 : std::move(r), std::move(u), std::move(v), std::move(w), std::move(e),
1373 [ + - ][ + + ]: 889 : std::move(pr) };
[ - - ]
1374 [ + - ][ + - ]: 127 : if (steady) nodefields.push_back( std::move(ma) );
1375 : :
1376 [ - + ]: 127 : for (std::size_t c=0; c<ncomp-5; ++c) {
1377 [ - - ][ - - ]: 0 : nodefieldnames.push_back( "c" + std::to_string(c) );
[ - - ]
1378 [ - - ][ - - ]: 0 : nodefields.push_back( m_u.extract(5+c) );
1379 : : }
1380 : :
1381 : : // query function to evaluate analytic solution (if defined)
1382 [ + - ]: 127 : auto sol = problems::SOL();
1383 : :
1384 [ - + ]: 127 : if (sol) {
1385 : 0 : const auto& coord = d->Coord();
1386 : 0 : const auto& x = coord[0];
1387 : 0 : const auto& y = coord[1];
1388 : 0 : const auto& z = coord[2];
1389 [ - - ]: 0 : auto an = m_u;
1390 [ - - ]: 0 : std::vector< tk::real > ap( m_u.nunk() );
1391 [ - - ]: 0 : for (std::size_t i=0; i<an.nunk(); ++i) {
1392 [ - - ]: 0 : auto s = sol( x[i], y[i], z[i], d->T(), /*meshid=*/0 );
1393 : 0 : s[1] /= s[0];
1394 : 0 : s[2] /= s[0];
1395 : 0 : s[3] /= s[0];
1396 : 0 : s[4] /= s[0];
1397 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
1398 : 0 : s[4] -= 0.5*(s[1]*s[1] + s[2]*s[2] + s[3]*s[3]);
1399 : 0 : ap[i] = eos::pressure( s[0]*s[4] );
1400 : 0 : }
1401 [ - - ]: 0 : for (std::size_t c=0; c<5; ++c) {
1402 [ - - ][ - - ]: 0 : nodefieldnames.push_back( nodefieldnames[c] + "_analytic" );
1403 [ - - ][ - - ]: 0 : nodefields.push_back( an.extract(c) );
1404 : : }
1405 [ - - ][ - - ]: 0 : nodefieldnames.push_back( nodefieldnames[5] + "_analytic" );
1406 [ - - ]: 0 : nodefields.push_back( std::move(ap) );
1407 [ - - ]: 0 : for (std::size_t c=0; c<ncomp-5; ++c) {
1408 [ - - ][ - - ]: 0 : nodefieldnames.push_back( nodefieldnames[6+c] + "_analytic" );
1409 [ - - ][ - - ]: 0 : nodefields.push_back( an.extract(5+c) );
1410 : : }
1411 : 0 : }
1412 : :
1413 [ - + ][ - - ]: 127 : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
1414 : :
1415 : : // Surface output
1416 : :
1417 : 127 : std::vector< std::string > nodesurfnames;
1418 : 127 : std::vector< std::vector< tk::real > > nodesurfs;
1419 : :
1420 : 127 : const auto& f = g_cfg.get< tag::fieldout, tag::sidesets >();
1421 : :
1422 [ - + ]: 127 : if (!f.empty()) {
1423 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "density" );
1424 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "velocityx" );
1425 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "velocityy" );
1426 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "velocityz" );
1427 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "energy" );
1428 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "pressure" );
1429 : :
1430 [ - - ]: 0 : for (std::size_t c=0; c<ncomp-5; ++c) {
1431 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "c" + std::to_string(c) );
[ - - ]
1432 : : }
1433 : :
1434 [ - - ][ - - ]: 0 : if (steady) nodesurfnames.push_back( "mach" );
[ - - ]
1435 : :
1436 [ - - ]: 0 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
1437 [ - - ]: 0 : std::set< int > outsets( begin(f), end(f) );
1438 [ - - ]: 0 : for (auto sideset : outsets) {
1439 [ - - ]: 0 : auto b = bnode.find(sideset);
1440 [ - - ]: 0 : if (b == end(bnode)) continue;
1441 : 0 : const auto& nodes = b->second;
1442 : 0 : auto i = nodesurfs.size();
1443 : 0 : auto ns = ncomp + 1;
1444 [ - - ]: 0 : if (steady) ++ns;
1445 [ - - ]: 0 : nodesurfs.insert( end(nodesurfs), ns,
1446 [ - - ]: 0 : std::vector< tk::real >( nodes.size() ) );
1447 : 0 : std::size_t j = 0;
1448 [ - - ]: 0 : for (auto n : nodes) {
1449 [ - - ]: 0 : const auto s = m_u[n];
1450 : 0 : std::size_t p = 0;
1451 : 0 : nodesurfs[i+(p++)][j] = s[0];
1452 : 0 : nodesurfs[i+(p++)][j] = s[1]/s[0];
1453 : 0 : nodesurfs[i+(p++)][j] = s[2]/s[0];
1454 : 0 : nodesurfs[i+(p++)][j] = s[3]/s[0];
1455 : 0 : nodesurfs[i+(p++)][j] = s[4]/s[0];
1456 : 0 : auto vv = (s[1]*s[1] + s[2]*s[2] + s[3]*s[3])/s[0]/s[0];
1457 : 0 : auto ei = s[4]/s[0] - 0.5*vv;
1458 : 0 : auto sp = eos::pressure( s[0]*ei );
1459 : 0 : nodesurfs[i+(p++)][j] = sp;
1460 [ - - ]: 0 : for (std::size_t c=0; c<ncomp-5; ++c) nodesurfs[i+(p++)+c][j] = s[5+c];
1461 [ - - ]: 0 : if (steady) {
1462 : 0 : nodesurfs[i+(p++)][j] = std::sqrt(vv) / eos::soundspeed( s[0], sp );
1463 : : }
1464 : 0 : ++j;
1465 : 0 : }
1466 : : }
1467 : 0 : }
1468 : :
1469 : : // Send mesh and fields data (solution dump) for output to file
1470 [ + - ][ + - ]: 254 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
1471 : 127 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
1472 : : {}, nodefields, {}, nodesurfs, cb );
1473 [ + - ][ + - ]: 508 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ - - ]
[ - - ][ - - ]
[ - - ]
1474 : :
1475 : : void
1476 : 3450 : LaxCG::out()
1477 : : // *****************************************************************************
1478 : : // Output mesh field data
1479 : : // *****************************************************************************
1480 : : {
1481 : 3450 : auto d = Disc();
1482 : :
1483 : : // Time history
1484 [ + - ][ + - ]: 3450 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ][ - + ]
1485 : 0 : auto ncomp = m_u.nprop();
1486 : 0 : const auto& inpoel = d->Inpoel();
1487 [ - - ]: 0 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
1488 : 0 : std::size_t j = 0;
1489 [ - - ]: 0 : for (const auto& p : d->Hist()) {
1490 : 0 : auto e = p.get< tag::elem >(); // host element id
1491 : 0 : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
1492 [ - - ]: 0 : hist[j].resize( ncomp+1, 0.0 );
1493 [ - - ]: 0 : for (std::size_t i=0; i<4; ++i) {
1494 [ - - ]: 0 : const auto u = m_u[ inpoel[e*4+i] ];
1495 : 0 : hist[j][0] += n[i] * u[0];
1496 : 0 : hist[j][1] += n[i] * u[1]/u[0];
1497 : 0 : hist[j][2] += n[i] * u[2]/u[0];
1498 : 0 : hist[j][3] += n[i] * u[3]/u[0];
1499 : 0 : hist[j][4] += n[i] * u[4]/u[0];
1500 : 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];
1501 : 0 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
1502 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
1503 : 0 : }
1504 : 0 : ++j;
1505 : : }
1506 [ - - ]: 0 : d->history( std::move(hist) );
1507 : 0 : }
1508 : :
1509 : : // Field data
1510 [ + + ][ + - ]: 3450 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
[ + + ]
1511 [ + - ][ + - ]: 390 : writeFields( CkCallback(CkIndex_LaxCG::integrals(), thisProxy[thisIndex]) );
[ + - ][ + - ]
1512 : : } else {
1513 : 3060 : integrals();
1514 : : }
1515 : 3450 : }
1516 : :
1517 : : void
1518 : 3450 : LaxCG::integrals()
1519 : : // *****************************************************************************
1520 : : // Compute integral quantities for output
1521 : : // *****************************************************************************
1522 : : {
1523 : 3450 : auto d = Disc();
1524 : :
1525 [ + - ][ + - ]: 3450 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ][ - + ]
1526 : :
1527 : : using namespace integrals;
1528 [ - - ]: 0 : std::vector< std::map< int, tk::real > > ints( NUMINT );
1529 : :
1530 : : // Prepend integral vector with metadata on the current time step:
1531 : : // current iteration count, current physical time, time step size
1532 [ - - ]: 0 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
1533 [ - - ]: 0 : ints[ TIME ][ 0 ] = d->T();
1534 [ - - ]: 0 : ints[ DT ][ 0 ] = d->Dt();
1535 : :
1536 : : // Compute integrals requested for surfaces requested
1537 : 0 : const auto& reqv = g_cfg.get< tag::integout, tag::integrals >();
1538 [ - - ]: 0 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
1539 [ - - ][ - - ]: 0 : if (req.count("mass_flow_rate")) {
[ - - ]
1540 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
1541 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
1542 : 0 : const auto& nodes = sint.first;
1543 : 0 : const auto& ndA = sint.second;
1544 : 0 : auto n = ndA.data();
1545 [ - - ]: 0 : for (auto p : nodes) {
1546 [ - - ][ - - ]: 0 : mfr += n[0]*m_u(p,1) + n[1]*m_u(p,2) + n[2]*m_u(p,3);
[ - - ]
1547 : 0 : n += 3;
1548 : : }
1549 : : }
1550 : : }
1551 : :
1552 [ - - ]: 0 : auto stream = serialize( d->MeshId(), ints );
1553 [ - - ]: 0 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
1554 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
1555 : :
1556 : 0 : } else {
1557 : :
1558 : 3450 : step();
1559 : :
1560 : : }
1561 : 3450 : }
1562 : :
1563 : : void
1564 : 10350 : LaxCG::stage()
1565 : : // *****************************************************************************
1566 : : // Evaluate whether to continue with next time step stage
1567 : : // *****************************************************************************
1568 : : {
1569 : : // Increment Runge-Kutta stage counter
1570 : 10350 : ++m_stage;
1571 : :
1572 : : // If not all Runge-Kutta stages complete, continue to next time stage,
1573 : : // otherwise output field data to file(s)
1574 [ + + ]: 10350 : if (m_stage < 3) grad(); else out();
1575 : 10350 : }
1576 : :
1577 : : void
1578 : 3129 : LaxCG::evalLB( int nrestart )
1579 : : // *****************************************************************************
1580 : : // Evaluate whether to do load balancing
1581 : : //! \param[in] nrestart Number of times restarted
1582 : : // *****************************************************************************
1583 : : {
1584 : 3129 : auto d = Disc();
1585 : :
1586 : : // Detect if just returned from a checkpoint and if so, zero timers and
1587 : : // finished flag
1588 [ + + ]: 3129 : if (d->restarted( nrestart )) m_finished = 0;
1589 : :
1590 : 3129 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
1591 : 3129 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
1592 : :
1593 : : // Load balancing if user frequency is reached or after the second time-step
1594 [ + + ][ + + ]: 3129 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
[ + + ]
1595 : :
1596 : 2844 : AtSync();
1597 [ - + ]: 2844 : if (nonblocking) dt();
1598 : :
1599 : : } else {
1600 : :
1601 : 285 : dt();
1602 : :
1603 : : }
1604 : 3129 : }
1605 : :
1606 : : void
1607 : 3124 : LaxCG::evalRestart()
1608 : : // *****************************************************************************
1609 : : // Evaluate whether to save checkpoint/restart
1610 : : // *****************************************************************************
1611 : : {
1612 : 3124 : auto d = Disc();
1613 : :
1614 : 3124 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1615 : 3124 : const auto benchmark = g_cfg.get< tag::benchmark >();
1616 : :
1617 [ + + ][ - + ]: 3124 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
1618 : :
1619 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1620 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1621 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
1622 : :
1623 : 0 : } else {
1624 : :
1625 : 3124 : evalLB( /* nrestart = */ -1 );
1626 : :
1627 : : }
1628 : 3124 : }
1629 : :
1630 : : void
1631 : 3450 : LaxCG::step()
1632 : : // *****************************************************************************
1633 : : // Evaluate whether to continue with next time step
1634 : : // *****************************************************************************
1635 : : {
1636 : 3450 : auto d = Disc();
1637 : :
1638 : : // Output one-liner status report to screen
1639 : 3450 : d->status();
1640 : : // Reset Runge-Kutta stage counter
1641 : 3450 : m_stage = 0;
1642 : :
1643 [ + + ]: 3450 : if (not m_finished) {
1644 : :
1645 : 3124 : evalRestart();
1646 : :
1647 : : } else {
1648 : :
1649 : 326 : auto meshid = d->MeshId();
1650 [ + - ]: 326 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1651 [ + - ][ + - ]: 652 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1652 : :
1653 : : }
1654 : 3450 : }
1655 : :
1656 : : #include "NoWarning/laxcg.def.h"
|
