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