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