Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ALECG.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief ALECG for a PDE system with continuous Galerkin + ALE + RK
9 : : \details ALECG advances a system of partial differential equations (PDEs)
10 : : using a continuous Galerkin (CG) finite element (FE) spatial discretization
11 : : (using linear shapefunctions on tetrahedron elements) combined with a
12 : : Runge-Kutta (RK) time stepping scheme in the arbitrary Eulerian-Lagrangian
13 : : reference frame.
14 : : \see The documentation in ALECG.hpp.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include "QuinoaBuildConfig.hpp"
19 : : #include "ALECG.hpp"
20 : : #include "Vector.hpp"
21 : : #include "Reader.hpp"
22 : : #include "ContainerUtil.hpp"
23 : : #include "UnsMesh.hpp"
24 : : #include "ExodusIIMeshWriter.hpp"
25 : : #include "Inciter/InputDeck/InputDeck.hpp"
26 : : #include "DerivedData.hpp"
27 : : #include "CGPDE.hpp"
28 : : #include "Discretization.hpp"
29 : : #include "DiagReducer.hpp"
30 : : #include "NodeBC.hpp"
31 : : #include "Refiner.hpp"
32 : : #include "Reorder.hpp"
33 : : #include "Around.hpp"
34 : : #include "CGPDE.hpp"
35 : : #include "Integrate/Mass.hpp"
36 : : #include "FieldOutput.hpp"
37 : :
38 : : #ifdef HAS_ROOT
39 : : #include "RootMeshWriter.hpp"
40 : : #endif
41 : :
42 : : namespace inciter {
43 : :
44 : : extern ctr::InputDeck g_inputdeck;
45 : : extern ctr::InputDeck g_inputdeck_defaults;
46 : : extern std::vector< CGPDE > g_cgpde;
47 : :
48 : : //! Runge-Kutta coefficients
49 : : static const std::array< tk::real, 3 > rkcoef{{ 1.0/3.0, 1.0/2.0, 1.0 }};
50 : :
51 : : } // inciter::
52 : :
53 : : using inciter::ALECG;
54 : :
55 : 534 : ALECG::ALECG( const CProxy_Discretization& disc,
56 : : const std::map< int, std::vector< std::size_t > >& bface,
57 : : const std::map< int, std::vector< std::size_t > >& bnode,
58 : 534 : const std::vector< std::size_t >& triinpoel ) :
59 : : m_disc( disc ),
60 : : m_nsol( 0 ),
61 : : m_ngrad( 0 ),
62 : : m_nrhs( 0 ),
63 : : m_nbnorm( 0 ),
64 : : m_ndfnorm( 0 ),
65 : : m_bnode( bnode ),
66 : : m_bface( bface ),
67 : : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
68 : : m_bndel( Disc()->bndel() ),
69 : : m_dfnorm(),
70 : : m_dfnormc(),
71 : : m_dfn(),
72 : : m_esup( tk::genEsup( Disc()->Inpoel(), 4 ) ),
73 : 534 : m_psup( tk::genPsup( Disc()->Inpoel(), 4, m_esup ) ),
74 [ + - ]: 534 : m_u( Disc()->Gid().size(), g_inputdeck.get< tag::component >().nprop() ),
75 : : m_un( m_u.nunk(), m_u.nprop() ),
76 : : m_rhs( m_u.nunk(), m_u.nprop() ),
77 : : m_rhsc(),
78 : 534 : m_chBndGrad( Disc()->Bid().size(), m_u.nprop()*3 ),
79 : : m_dirbc(),
80 : : m_chBndGradc(),
81 : : m_diag(),
82 : : m_bnorm(),
83 : : m_bnormc(),
84 : : m_symbcnodes(),
85 : : m_farfieldbcnodes(),
86 : : m_symbctri(),
87 : : m_spongenodes(),
88 : : m_timedepbcnodes(),
89 : : m_timedepbcFn(),
90 : : m_stage( 0 ),
91 : : m_boxnodes(),
92 : : m_edgenode(),
93 : : m_edgeid(),
94 : 0 : m_dtp( m_u.nunk(), 0.0 ),
95 : 534 : m_tp( m_u.nunk(), g_inputdeck.get< tag::discr, tag::t0 >() ),
96 : : m_finished( 0 ),
97 [ + - ][ + - ]: 2136 : m_newmesh( 0 )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
98 : : // *****************************************************************************
99 : : // Constructor
100 : : //! \param[in] disc Discretization proxy
101 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
102 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
103 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
104 : : // *****************************************************************************
105 : : //! [Constructor]
106 : : {
107 : 534 : usesAtSync = true; // enable migration at AtSync
108 : :
109 [ + - ]: 534 : auto d = Disc();
110 : :
111 : : // Perform optional operator-access-pattern mesh node reordering
112 [ + + ]: 534 : if (g_inputdeck.get< tag::discr, tag::operator_reorder >()) {
113 : :
114 : : // Create new local ids based on access pattern of PDE operators
115 : 20 : std::unordered_map< std::size_t, std::size_t > map;
116 : 10 : std::size_t n = 0;
117 : :
118 [ + + ]: 1104 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
119 [ + - ][ + + ]: 1094 : if (map.find(p) == end(map)) map[p] = n++;
[ + - ]
120 [ + + ]: 11110 : for (auto q : tk::Around(m_psup,p)) { // for each edge p-q
121 [ + - ][ + + ]: 10016 : if (map.find(q) == end(map)) map[q] = n++;
[ + - ]
122 : : }
123 : : }
124 : :
125 [ - + ][ - - ]: 10 : Assert( map.size() == d->Gid().size(), "Map size mismatch" );
[ - - ][ - - ]
126 : :
127 : : // Remap data in bound Discretization object
128 [ + - ]: 10 : d->remap( map );
129 : : // Recompute elements surrounding points
130 [ + - ]: 10 : m_esup = tk::genEsup( d->Inpoel(), 4 );
131 : : // Recompute points surrounding points
132 [ + - ]: 10 : m_psup = tk::genPsup( d->Inpoel(), 4, m_esup );
133 : : // Remap boundary triangle face connectivity
134 [ + - ]: 10 : tk::remap( m_triinpoel, map );
135 : : }
136 : :
137 : : // Query/update boundary-conditions-related data structures from user input
138 [ + - ]: 534 : queryBnd();
139 : :
140 : : // Activate SDAG wait for initially computing normals
141 [ + - ][ + - ]: 534 : thisProxy[ thisIndex ].wait4norm();
142 : :
143 : : // Generate callbacks for solution transfers we are involved in
144 : :
145 : : // Always add a callback to be used when we are not involved in any transfers
146 : 1068 : std::vector< CkCallback > cb;
147 [ + - ][ + - ]: 1068 : auto c = CkCallback(CkIndex_ALECG::transfer_complete(), thisProxy[thisIndex]);
[ + - ]
148 [ + - ]: 534 : cb.push_back( c );
149 : :
150 : : // Generate a callback for each transfer we are involved in (either as a
151 : : // source or a destination)
152 : 534 : auto meshid = d->MeshId();
153 [ - + ]: 534 : for (const auto& t : d->Transfers())
154 [ - - ][ - - ]: 0 : if (meshid == t.src || meshid == t.dst)
155 [ - - ]: 0 : cb.push_back( c );
156 : :
157 : : // Send callbacks to base
158 [ + - ]: 534 : d->transferCallback( cb );
159 : 534 : }
160 : : //! [Constructor]
161 : :
162 : : void
163 : 47340 : ALECG::queryBnd()
164 : : // *****************************************************************************
165 : : // Query/update boundary-conditions-related data structures from user input
166 : : // *****************************************************************************
167 : : {
168 [ + - ]: 47340 : auto d = Disc();
169 : :
170 : : // Query and match user-specified Dirichlet boundary conditions to side sets
171 : 47340 : const auto steady = g_inputdeck.get< tag::discr, tag::steady_state >();
172 [ + + ][ + + ]: 114858 : if (steady) for (auto& deltat : m_dtp) deltat *= rkcoef[m_stage];
173 [ + - ]: 94680 : m_dirbc = match( m_u.nprop(), d->T(), rkcoef[m_stage] * d->Dt(),
174 : 47340 : m_tp, m_dtp, d->Coord(), d->Lid(), m_bnode,
175 : 47340 : /* increment = */ false );
176 [ + + ][ + + ]: 114858 : if (steady) for (auto& deltat : m_dtp) deltat /= rkcoef[m_stage];
177 : :
178 : : // Prepare unique set of symmetry BC nodes
179 [ + - ]: 94680 : auto sym = d->bcnodes< tag::bc, tag::bcsym >( m_bface, m_triinpoel );
180 [ + + ]: 74192 : for (const auto& [s,nodes] : sym)
181 [ + - ]: 26852 : m_symbcnodes.insert( begin(nodes), end(nodes) );
182 : :
183 : : // Prepare unique set of farfield BC nodes
184 [ + - ]: 94680 : auto far = d->bcnodes< tag::bc, tag::bcfarfield >( m_bface, m_triinpoel );
185 [ + + ]: 48365 : for (const auto& [s,nodes] : far)
186 [ + - ]: 1025 : m_farfieldbcnodes.insert( begin(nodes), end(nodes) );
187 : :
188 : : // If farfield BC is set on a node, will not also set symmetry BC
189 [ + + ][ + - ]: 129135 : for (auto fn : m_farfieldbcnodes) m_symbcnodes.erase(fn);
190 : :
191 : : // Prepare boundary nodes contiguously accessible from a triangle-face loop
192 [ + - ]: 47340 : m_symbctri.resize( m_triinpoel.size()/3, 0 );
193 [ + + ]: 5599902 : for (std::size_t e=0; e<m_triinpoel.size()/3; ++e)
194 [ + - ][ + + ]: 5552562 : if (m_symbcnodes.find(m_triinpoel[e*3+0]) != end(m_symbcnodes))
195 : 4408182 : m_symbctri[e] = 1;
196 : :
197 : : // Prepare unique set of sponge nodes
198 [ + - ]: 94680 : auto sponge = d->bcnodes< tag::sponge, tag::sideset >( m_bface, m_triinpoel );
199 [ + + ]: 47402 : for (const auto& [s,nodes] : sponge)
200 [ + - ]: 62 : m_spongenodes.insert( begin(nodes), end(nodes) );
201 : :
202 : : // Prepare unique set of time dependent BC nodes
203 : 47340 : m_timedepbcnodes.clear();
204 : 47340 : m_timedepbcFn.clear();
205 : : const auto& timedep =
206 : 47340 : g_inputdeck.template get< tag::param, tag::compflow, tag::bctimedep >();
207 [ + + ]: 47340 : if (!timedep.empty()) {
208 [ + - ]: 25081 : m_timedepbcnodes.resize(timedep[0].size());
209 [ + - ]: 25081 : m_timedepbcFn.resize(timedep[0].size());
210 : 25081 : std::size_t ib=0;
211 [ + + ]: 25286 : for (const auto& bndry : timedep[0]) {
212 : 410 : std::unordered_set< std::size_t > nodes;
213 [ + + ]: 410 : for (const auto& s : bndry.template get< tag::sideset >()) {
214 [ + - ][ + - ]: 205 : auto k = m_bnode.find( std::stoi(s) );
215 [ + - ]: 205 : if (k != end(m_bnode)) {
216 [ + + ]: 2460 : for (auto g : k->second) { // global node ids on side set
217 [ + - ][ + - ]: 2255 : nodes.insert( tk::cref_find(d->Lid(),g) );
218 : : }
219 : : }
220 : : }
221 [ + - ]: 205 : m_timedepbcnodes[ib].insert( begin(nodes), end(nodes) );
222 : :
223 : : // Store user defined discrete function in time. This is done in the same
224 : : // loop as the BC nodes, so that the indices for the two vectors
225 : : // m_timedepbcnodes and m_timedepbcFn are consistent with each other
226 [ + - ]: 205 : auto fn = bndry.template get< tag::fn >();
227 [ + + ]: 820 : for (std::size_t ir=0; ir<fn.size()/6; ++ir) {
228 [ + - ]: 1230 : m_timedepbcFn[ib].push_back({{ fn[ir*6+0], fn[ir*6+1], fn[ir*6+2],
229 : 615 : fn[ir*6+3], fn[ir*6+4], fn[ir*6+5] }});
230 : : }
231 : 205 : ++ib;
232 : : }
233 : : }
234 : :
235 [ - + ][ - - ]: 47340 : Assert(m_timedepbcFn.size() == m_timedepbcnodes.size(), "Incorrect number of "
[ - - ][ - - ]
236 : : "time dependent functions.");
237 : :
238 : : // Query ALE mesh velocity boundary condition node lists and node lists at
239 : : // which ALE moves boundaries
240 [ + - ]: 47340 : d->meshvelBnd( m_bface, m_bnode, m_triinpoel );
241 : 47340 : }
242 : :
243 : : void
244 : 1614 : ALECG::norm()
245 : : // *****************************************************************************
246 : : // Start (re-)computing boundary point-, and dual-face normals
247 : : // *****************************************************************************
248 : : {
249 [ + - ]: 1614 : auto d = Disc();
250 : :
251 : : // Query nodes at which symmetry BCs are specified
252 [ + - ]: 3228 : auto bn = d->bcnodes< tag::bc, tag::bcsym >( m_bface, m_triinpoel );
253 : :
254 : : // Query nodes at which farfield BCs are specified
255 [ + - ]: 3228 : auto far = d->bcnodes< tag::bc, tag::bcfarfield >( m_bface, m_triinpoel );
256 : : // Merge BC data where boundary-point normals are required
257 [ + + ][ + - ]: 1619 : for (const auto& [s,n] : far) bn[s].insert( begin(n), end(n) );
[ + - ]
258 : :
259 : : // Query nodes at which mesh velocity symmetry BCs are specified
260 : 3228 : std::unordered_map<int, std::unordered_set< std::size_t >> ms;
261 [ + + ]: 3846 : for (const auto& s : g_inputdeck.template get< tag::ale, tag::bcsym >()) {
262 [ + - ][ + - ]: 2232 : auto k = m_bface.find( std::stoi(s) );
263 [ + + ]: 2232 : if (k != end(m_bface)) {
264 [ + - ]: 1550 : auto& n = ms[ k->first ];
265 [ + + ]: 392150 : for (auto f : k->second) {
266 [ + - ]: 390600 : n.insert( m_triinpoel[f*3+0] );
267 [ + - ]: 390600 : n.insert( m_triinpoel[f*3+1] );
268 [ + - ]: 390600 : n.insert( m_triinpoel[f*3+2] );
269 : : }
270 : : }
271 : : }
272 : : // Merge BC data where boundary-point normals are required
273 [ + + ][ + - ]: 3164 : for (const auto& [s,n] : ms) bn[s].insert( begin(n), end(n) );
[ + - ]
274 : :
275 : : // Compute boundary point normals
276 [ + - ]: 1614 : bnorm( bn );
277 : :
278 : : // Compute dual-face normals associated to edges
279 [ + - ]: 1614 : dfnorm();
280 : 1614 : }
281 : :
282 : : std::array< tk::real, 3 >
283 : 8206028 : ALECG::edfnorm( const tk::UnsMesh::Edge& edge,
284 : : const std::unordered_map< tk::UnsMesh::Edge,
285 : : std::vector< std::size_t >,
286 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> >& esued )
287 : : const
288 : : // *****************************************************************************
289 : : // Compute normal of dual-mesh associated to edge
290 : : //! \param[in] edge Edge whose dual-face normal to compute given by local ids
291 : : //! \param[in] esued Elements surrounding edges
292 : : //! \return Dual-face normal for edge
293 : : // *****************************************************************************
294 : : {
295 : 8206028 : auto d = Disc();
296 : 8206028 : const auto& inpoel = d->Inpoel();
297 : 8206028 : const auto& coord = d->Coord();
298 : 8206028 : const auto& x = coord[0];
299 : 8206028 : const auto& y = coord[1];
300 : 8206028 : const auto& z = coord[2];
301 : :
302 : 8206028 : std::array< tk::real, 3 > n{ 0.0, 0.0, 0.0 };
303 : :
304 [ + - ][ + + ]: 42527150 : for (auto e : tk::cref_find(esued,edge)) {
305 : : // access node IDs
306 : : const std::array< std::size_t, 4 >
307 : 34321122 : N{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
308 : : // compute element Jacobi determinant
309 : : const std::array< tk::real, 3 >
310 : 34321122 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
311 : 34321122 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
312 : 34321122 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
313 : 34321122 : const auto J = tk::triple( ba, ca, da ); // J = 6V
314 [ - + ][ - - ]: 34321122 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
315 : : // shape function derivatives, nnode*ndim [4][3]
316 : : std::array< std::array< tk::real, 3 >, 4 > grad;
317 : 34321122 : grad[1] = tk::crossdiv( ca, da, J );
318 : 34321122 : grad[2] = tk::crossdiv( da, ba, J );
319 : 34321122 : grad[3] = tk::crossdiv( ba, ca, J );
320 [ + + ]: 137284488 : for (std::size_t i=0; i<3; ++i)
321 : 102963366 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
322 : : // sum normal contributions
323 : : // The constant 1/48: Eq (12) from Waltz et al. Computers & fluids (92) 2014
324 : : // The result of the integral of shape function N on a tet is V/4.
325 : : // This can be written as J/(6*4). Eq (12) has a 1/2 multiplier.
326 : : // This leads to J/48.
327 : 34321122 : auto J48 = J/48.0;
328 [ + + ]: 240247854 : for (const auto& [a,b] : tk::lpoed) {
329 [ + - ]: 205926732 : auto s = tk::orient( {N[a],N[b]}, edge );
330 [ + + ]: 823706928 : for (std::size_t j=0; j<3; ++j)
331 : 617780196 : n[j] += J48 * s * (grad[a][j] - grad[b][j]);
332 : : }
333 : : }
334 : :
335 : 8206028 : return n;
336 : : }
337 : :
338 : : void
339 : 1614 : ALECG::dfnorm()
340 : : // *****************************************************************************
341 : : // Compute dual-face normals associated to edges
342 : : // *****************************************************************************
343 : : {
344 [ + - ]: 1614 : auto d = Disc();
345 : 1614 : const auto& inpoel = d->Inpoel();
346 : 1614 : const auto& gid = d->Gid();
347 : :
348 : : // compute derived data structures
349 [ + - ][ + - ]: 3228 : auto esued = tk::genEsued( inpoel, 4, tk::genEsup( inpoel, 4 ) );
350 : :
351 : : // Compute dual-face normals for domain edges
352 [ + + ]: 1497960 : for (std::size_t p=0; p<gid.size(); ++p) // for each point p
353 [ + + ]: 17908402 : for (auto q : tk::Around(m_psup,p)) // for each edge p-q
354 [ + + ]: 16412056 : if (gid[p] < gid[q])
355 [ + - ][ + - ]: 8206028 : m_dfnorm[{gid[p],gid[q]}] = edfnorm( {p,q}, esued );
356 : :
357 : : // Send our dual-face normal contributions to neighbor chares
358 [ + + ]: 1614 : if (d->EdgeCommMap().empty())
359 [ + - ]: 170 : comdfnorm_complete();
360 : : else {
361 [ + + ]: 9446 : for (const auto& [c,edges] : d->EdgeCommMap()) {
362 : 8002 : decltype(m_dfnorm) exp;
363 [ + + ][ + - ]: 973100 : for (const auto& e : edges) exp[e] = tk::cref_find(m_dfnorm,e);
[ + - ]
364 [ + - ][ + - ]: 8002 : thisProxy[c].comdfnorm( exp );
365 : : }
366 : : }
367 : :
368 [ + - ]: 1614 : owndfnorm_complete();
369 : 1614 : }
370 : :
371 : : void
372 : 8002 : ALECG::comdfnorm( const std::unordered_map< tk::UnsMesh::Edge,
373 : : std::array< tk::real, 3 >,
374 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> >& dfnorm )
375 : : // *****************************************************************************
376 : : // Receive contributions to dual-face normals on chare-boundaries
377 : : //! \param[in] dfnorm Incoming partial sums of dual-face normals associated to
378 : : //! chare-boundary edges
379 : : // *****************************************************************************
380 : : {
381 : : // Buffer up inccoming contributions to dual-face normals
382 [ + + ]: 973100 : for (const auto& [e,n] : dfnorm) {
383 [ + - ]: 965098 : auto& dfn = m_dfnormc[e];
384 : 965098 : dfn[0] += n[0];
385 : 965098 : dfn[1] += n[1];
386 : 965098 : dfn[2] += n[2];
387 : : }
388 : :
389 [ + + ]: 8002 : if (++m_ndfnorm == Disc()->EdgeCommMap().size()) {
390 : 1444 : m_ndfnorm = 0;
391 : 1444 : comdfnorm_complete();
392 : : }
393 : 8002 : }
394 : :
395 : : void
396 : 1614 : ALECG::bnorm( const std::unordered_map< int,
397 : : std::unordered_set< std::size_t > >& bcnodes )
398 : : // *****************************************************************************
399 : : // Compute boundary point normals
400 : : //! \param[in] bcnodes Local node ids associated to side set ids at which BCs
401 : : //! are set that require normals
402 : : //*****************************************************************************
403 : : {
404 : 1614 : auto d = Disc();
405 : :
406 : 1614 : m_bnorm = cg::bnorm( m_bface, m_triinpoel, d->Coord(), d->Gid(), bcnodes );
407 : :
408 : : // Send our nodal normal contributions to neighbor chares
409 [ + + ]: 1614 : if (d->NodeCommMap().empty())
410 : 170 : comnorm_complete();
411 : : else
412 [ + + ]: 9446 : for (const auto& [ neighborchare, sharednodes ] : d->NodeCommMap()) {
413 : : std::unordered_map< int,
414 : 8002 : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > exp;
415 [ + + ]: 410228 : for (auto i : sharednodes) {
416 [ + + ]: 1650474 : for (const auto& [s,norms] : m_bnorm) {
417 [ + - ]: 1248248 : auto j = norms.find(i);
418 [ + + ][ + - ]: 1248248 : if (j != end(norms)) exp[s][i] = j->second;
[ + - ]
419 : : }
420 : : }
421 [ + - ][ + - ]: 8002 : thisProxy[ neighborchare ].comnorm( exp );
422 : : }
423 : :
424 : 1614 : ownnorm_complete();
425 : 1614 : }
426 : :
427 : : void
428 : 8002 : ALECG::comnorm( const std::unordered_map< int,
429 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& innorm )
430 : : // *****************************************************************************
431 : : // Receive boundary point normals on chare-boundaries
432 : : //! \param[in] innorm Incoming partial sums of boundary point normal
433 : : //! contributions to normals (first 3 components), inverse distance squared
434 : : //! (4th component), associated to side set ids
435 : : // *****************************************************************************
436 : : {
437 : : // Buffer up incoming boundary-point normal vector contributions
438 [ + + ]: 14192 : for (const auto& [s,norms] : innorm) {
439 [ + - ]: 6190 : auto& bnorms = m_bnormc[s];
440 [ + + ]: 173730 : for (const auto& [p,n] : norms) {
441 [ + - ]: 167540 : auto& bnorm = bnorms[p];
442 : 167540 : bnorm[0] += n[0];
443 : 167540 : bnorm[1] += n[1];
444 : 167540 : bnorm[2] += n[2];
445 : 167540 : bnorm[3] += n[3];
446 : : }
447 : : }
448 : :
449 [ + + ]: 8002 : if (++m_nbnorm == Disc()->NodeCommMap().size()) {
450 : 1444 : m_nbnorm = 0;
451 : 1444 : comnorm_complete();
452 : : }
453 : 8002 : }
454 : :
455 : : void
456 : 666 : ALECG::registerReducers()
457 : : // *****************************************************************************
458 : : // Configure Charm++ reduction types initiated from this chare array
459 : : //! \details Since this is a [initnode] routine, the runtime system executes the
460 : : //! routine exactly once on every logical node early on in the Charm++ init
461 : : //! sequence. Must be static as it is called without an object. See also:
462 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
463 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
464 : : // *****************************************************************************
465 : : {
466 : 666 : NodeDiagnostics::registerReducers();
467 : 666 : }
468 : :
469 : : void
470 : 10306 : ALECG::ResumeFromSync()
471 : : // *****************************************************************************
472 : : // Return from migration
473 : : //! \details This is called when load balancing (LB) completes. The presence of
474 : : //! this function does not affect whether or not we block on LB.
475 : : // *****************************************************************************
476 : : {
477 [ - + ][ - - ]: 10306 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
478 : :
479 [ + - ]: 10306 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
480 : 10306 : }
481 : :
482 : : //! [setup]
483 : : void
484 : 534 : ALECG::setup()
485 : : // *****************************************************************************
486 : : // Start setup for solution
487 : : // *****************************************************************************
488 : : {
489 : 534 : auto d = Disc();
490 : :
491 : : // Determine nodes inside user-defined IC box
492 [ + + ][ + - ]: 1068 : for (auto& eq : g_cgpde) eq.IcBoxNodes( d->Coord(), m_boxnodes );
493 : :
494 : : // Compute volume of user-defined box IC
495 : 534 : d->boxvol( m_boxnodes );
496 : :
497 : : // Query time history field output labels from all PDEs integrated
498 : 534 : const auto& hist_points = g_inputdeck.get< tag::history, tag::point >();
499 [ + + ]: 534 : if (!hist_points.empty()) {
500 : 56 : std::vector< std::string > histnames;
501 [ + + ]: 56 : for (const auto& eq : g_cgpde) {
502 [ + - ]: 28 : auto n = eq.histNames();
503 [ + - ]: 28 : histnames.insert( end(histnames), begin(n), end(n) );
504 : : }
505 [ + - ]: 28 : d->histheader( std::move(histnames) );
506 : : }
507 : 534 : }
508 : : //! [setup]
509 : :
510 : : void
511 : 232012 : ALECG::volumetric( tk::Fields& u, const std::vector< tk::real >& v )
512 : : // *****************************************************************************
513 : : // Multiply solution with mesh volume
514 : : //! \param[in,out] u Solution vector
515 : : //! \param[in] v Volume to multiply with
516 : : // *****************************************************************************
517 : : {
518 [ - + ][ - - ]: 232012 : Assert( v.size() == u.nunk(), "Size mismatch" );
[ - - ][ - - ]
519 : :
520 [ + + ]: 49045317 : for (std::size_t i=0; i<u.nunk(); ++i)
521 [ + + ]: 180196974 : for (ncomp_t c=0; c<u.nprop(); ++c)
522 : 131383669 : u(i,c,0) *= v[i];
523 : 232012 : }
524 : :
525 : : void
526 : 231478 : ALECG::conserved( tk::Fields& u, const std::vector< tk::real >& v )
527 : : // *****************************************************************************
528 : : // Divide solution with mesh volume
529 : : //! \param[in,out] u Solution vector
530 : : //! \param[in] v Volume to divide with
531 : : // *****************************************************************************
532 : : {
533 [ - + ][ - - ]: 231478 : Assert( v.size() == u.nunk(), "Size mismatch" );
[ - - ][ - - ]
534 : :
535 [ + + ]: 48844107 : for (std::size_t i=0; i<u.nunk(); ++i)
536 [ + + ]: 179626638 : for (ncomp_t c=0; c<u.nprop(); ++c) {
537 : 131014009 : u(i,c,0) /= v[i];
538 : : }
539 : 231478 : }
540 : :
541 : : void
542 : 534 : ALECG::box( tk::real v )
543 : : // *****************************************************************************
544 : : // Receive total box IC volume and set conditions in box
545 : : //! \param[in] v Total volume within user-specified box
546 : : // *****************************************************************************
547 : : {
548 : 534 : auto d = Disc();
549 : :
550 : : // Store user-defined box IC volume
551 : 534 : d->Boxvol() = v;
552 : :
553 : : // Set initial conditions for all PDEs
554 [ + + ]: 1068 : for (auto& eq : g_cgpde)
555 [ + - ]: 534 : eq.initialize( d->Coord(), m_u, d->T(), d->Boxvol(), m_boxnodes );
556 : :
557 : : // Multiply conserved variables with mesh volume
558 : 534 : volumetric( m_u, Disc()->Vol() );
559 : :
560 : : // Initiate IC transfer (if coupled)
561 : 534 : Disc()->transfer( m_u );
562 : :
563 : : // Initialize nodal mesh volumes at previous time step stage
564 : 534 : d->Voln() = d->Vol();
565 : :
566 : : // Start computing the mesh mesh velocity for ALE
567 : 534 : meshvelstart();
568 : 534 : }
569 : :
570 : : void
571 : 47340 : ALECG::meshvelstart()
572 : : // *****************************************************************************
573 : : // Start computing the mesh mesh velocity for ALE
574 : : // *****************************************************************************
575 : : {
576 [ + - ]: 47340 : auto d = Disc();
577 : :
578 : : // Apply boundary conditions on numerical solution
579 [ + - ]: 47340 : BC();
580 : :
581 [ + - ]: 47340 : conserved( m_u, d->Vol() );
582 : :
583 : : // query fluid velocity across all systems integrated
584 : 94680 : tk::UnsMesh::Coords vel;
585 [ + + ][ + - ]: 94680 : for (const auto& eq : g_cgpde) eq.velocity( m_u, vel );
586 : : // query speed of sound in mesh nodes across all systems integrated
587 : 47340 : std::vector< tk::real > soundspeed;
588 [ + + ][ + - ]: 94680 : for (const auto& eq : g_cgpde) eq.soundspeed( m_u, soundspeed );
589 : :
590 [ + - ]: 47340 : volumetric( m_u, d->Vol() );
591 : :
592 : : // Start computing the mesh mesh velocity for ALE
593 [ + - ][ + - ]: 47340 : d->meshvelStart( vel, soundspeed, m_bnorm, rkcoef[m_stage] * d->Dt(),
594 [ + - ][ + - ]: 94680 : CkCallback(CkIndex_ALECG::meshveldone(), thisProxy[thisIndex]) );
[ + - ]
595 : 47340 : }
596 : :
597 : : void
598 : 47340 : ALECG::meshveldone()
599 : : // *****************************************************************************
600 : : // Done with computing the mesh velocity for ALE
601 : : // *****************************************************************************
602 : : {
603 : : // Assess and record mesh velocity linear solver conergence
604 : 47340 : Disc()->meshvelConv();
605 : :
606 : : // Continue
607 [ + + ]: 47340 : if (Disc()->Initial()) lhs(); else ale();
608 : 47340 : }
609 : :
610 : : //! [start]
611 : : void
612 : 534 : ALECG::start()
613 : : // *****************************************************************************
614 : : // Start time stepping
615 : : // *****************************************************************************
616 : : {
617 : : // Set flag that indicates that we are now during time stepping
618 : 534 : Disc()->Initial( 0 );
619 : : // Start timer measuring time stepping wall clock time
620 : 534 : Disc()->Timer().zero();
621 : : // Zero grind-timer
622 : 534 : Disc()->grindZero();
623 : : // Continue to first time step
624 : 534 : next();
625 : 534 : }
626 : : //! [start]
627 : :
628 : : //! [Compute lhs]
629 : : void
630 : 1614 : ALECG::lhs()
631 : : // *****************************************************************************
632 : : // Compute the left-hand side of transport equations
633 : : //! \details Also (re-)compute all data structures if the mesh changed.
634 : : // *****************************************************************************
635 : : {
636 : : // No need for LHS in ALECG
637 : :
638 : : // (Re-)compute boundary point-, and dual-face normals
639 : 1614 : norm();
640 : 1614 : }
641 : : //! [Compute lhs]
642 : :
643 : : //! [Merge normals and continue]
644 : : void
645 : 1614 : ALECG::mergelhs()
646 : : // *****************************************************************************
647 : : // The own and communication portion of the left-hand side is complete
648 : : // *****************************************************************************
649 : : {
650 : : // Combine own and communicated contributions of normals
651 : 1614 : normfinal();
652 : :
653 [ + + ]: 1614 : if (Disc()->Initial()) {
654 : : // Output initial conditions to file
655 [ + - ][ + - ]: 534 : writeFields( CkCallback(CkIndex_ALECG::start(), thisProxy[thisIndex]) );
[ + - ]
656 : : } else {
657 : 1080 : norm_complete();
658 : : }
659 : 1614 : }
660 : : //! [Merge normals and continue]
661 : :
662 : : void
663 : 1614 : ALECG::normfinal()
664 : : // *****************************************************************************
665 : : // Finish computing dual-face and boundary point normals
666 : : // *****************************************************************************
667 : : {
668 [ + - ]: 1614 : auto d = Disc();
669 : 1614 : const auto& lid = d->Lid();
670 : :
671 : : // Combine own and communicated contributions to boundary point normals
672 [ + + ]: 4451 : for (const auto& [s,norms] : m_bnormc) {
673 [ + - ]: 2837 : auto& bnorms = m_bnorm[s];
674 [ + + ]: 159488 : for (const auto& [p,n] : norms) {
675 [ + - ]: 156651 : auto& norm = bnorms[p];
676 : 156651 : norm[0] += n[0];
677 : 156651 : norm[1] += n[1];
678 : 156651 : norm[2] += n[2];
679 : 156651 : norm[3] += n[3];
680 : : }
681 : : }
682 : 1614 : tk::destroy( m_bnormc );
683 : :
684 : : // Divide summed point normals by the sum of inverse distance squared
685 [ + + ]: 4860 : for (auto& [s,norms] : m_bnorm)
686 [ + + ]: 777710 : for (auto& [p,n] : norms) {
687 : 774464 : n[0] /= n[3];
688 : 774464 : n[1] /= n[3];
689 : 774464 : n[2] /= n[3];
690 [ - + ][ - - ]: 774464 : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
[ - - ][ - - ]
691 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
692 : : "Non-unit normal" );
693 : : }
694 : :
695 : : // Replace global->local ids associated to boundary point normals
696 : 3228 : decltype(m_bnorm) bnorm;
697 [ + + ]: 4860 : for (auto& [s,norms] : m_bnorm) {
698 [ + - ]: 3246 : auto& bnorms = bnorm[s];
699 [ + + ]: 777710 : for (auto&& [g,n] : norms)
700 [ + - ][ + - ]: 774464 : bnorms[ tk::cref_find(lid,g) ] = std::move(n);
701 : : }
702 : 1614 : m_bnorm = std::move(bnorm);
703 : :
704 : : // Count contributions to chare-boundary edges
705 : : std::unordered_map< tk::UnsMesh::Edge, std::size_t,
706 : 3228 : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > edge_node_count;
707 [ + + ]: 9616 : for (const auto& [c,edges] : d->EdgeCommMap())
708 [ + + ]: 973100 : for (const auto& e : edges)
709 [ + - ]: 965098 : ++edge_node_count[e];
710 : :
711 : : // Combine and weigh communicated contributions to dual-face normals
712 [ + + ]: 899372 : for (auto& [e,n] : m_dfnormc) {
713 [ + - ]: 897758 : const auto& dfn = tk::cref_find( m_dfnorm, e );
714 : 897758 : n[0] += dfn[0];
715 : 897758 : n[1] += dfn[1];
716 : 897758 : n[2] += dfn[2];
717 [ + - ]: 897758 : auto count = static_cast< tk::real >( tk::cref_find( edge_node_count, e ) );
718 : 897758 : auto factor = 1.0/(count + 1.0);
719 [ + + ]: 3591032 : for (auto & x : n) x *= factor;
720 : : }
721 : :
722 : : // Generate list of unique edges
723 : 3228 : tk::UnsMesh::EdgeSet uedge;
724 [ + + ]: 1497960 : for (std::size_t p=0; p<m_u.nunk(); ++p)
725 [ + + ]: 17908402 : for (auto q : tk::Around(m_psup,p))
726 [ + - ]: 16412056 : uedge.insert( {p,q} );
727 : :
728 : : // Flatten edge list
729 [ + - ]: 1614 : m_edgenode.resize( uedge.size() * 2 );
730 : 1614 : std::size_t f = 0;
731 : 1614 : const auto& gid = d->Gid();
732 [ + + ]: 8207642 : for (auto&& [p,q] : uedge) {
733 [ + + ]: 8206028 : if (gid[p] > gid[q]) {
734 : 97633 : m_edgenode[f+0] = std::move(q);
735 : 97633 : m_edgenode[f+1] = std::move(p);
736 : : } else {
737 : 8108395 : m_edgenode[f+0] = std::move(p);
738 : 8108395 : m_edgenode[f+1] = std::move(q);
739 : : }
740 : 8206028 : f += 2;
741 : : }
742 : 1614 : tk::destroy(uedge);
743 : :
744 : : // Convert dual-face normals to streamable (and vectorizable) data structure
745 [ + - ]: 1614 : m_dfn.resize( m_edgenode.size() * 3 ); // 2 vectors per access
746 : : std::unordered_map< tk::UnsMesh::Edge, std::size_t,
747 : 3228 : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > eid;
748 [ + + ]: 8207642 : for (std::size_t e=0; e<m_edgenode.size()/2; ++e) {
749 : 8206028 : auto p = m_edgenode[e*2+0];
750 : 8206028 : auto q = m_edgenode[e*2+1];
751 [ + - ]: 8206028 : eid[{p,q}] = e;
752 : 8206028 : std::array< std::size_t, 2 > g{ gid[p], gid[q] };
753 [ + - ]: 8206028 : auto n = tk::cref_find( m_dfnorm, g );
754 : : // figure out if this is an edge on the parallel boundary
755 [ + - ]: 8206028 : auto nit = m_dfnormc.find( g );
756 [ + + ]: 8206028 : auto m = ( nit != m_dfnormc.end() ) ? nit->second : n;
757 : 8206028 : m_dfn[e*6+0] = n[0];
758 : 8206028 : m_dfn[e*6+1] = n[1];
759 : 8206028 : m_dfn[e*6+2] = n[2];
760 : 8206028 : m_dfn[e*6+3] = m[0];
761 : 8206028 : m_dfn[e*6+4] = m[1];
762 : 8206028 : m_dfn[e*6+5] = m[2];
763 : : }
764 : :
765 : 1614 : tk::destroy( m_dfnorm );
766 : 1614 : tk::destroy( m_dfnormc );
767 : :
768 : : // Flatten edge id data structure
769 [ + - ]: 1614 : m_edgeid.resize( m_psup.first.size() );
770 [ + + ]: 1497960 : for (std::size_t p=0,k=0; p<m_u.nunk(); ++p)
771 [ + + ]: 17908402 : for (auto q : tk::Around(m_psup,p))
772 [ + - ]: 16412056 : m_edgeid[k++] = tk::cref_find( eid, {p,q} );
773 : 1614 : }
774 : :
775 : : void
776 : 47340 : ALECG::BC()
777 : : // *****************************************************************************
778 : : // Apply boundary conditions
779 : : // \details The following BC enforcement changes the initial condition or
780 : : //! updated solution (dependending on when it is called) to ensure strong
781 : : //! imposition of the BCs. This is a matter of choice. Another alternative is
782 : : //! to only apply BCs when computing fluxes at boundary faces, thereby only
783 : : //! weakly enforcing the BCs. The former is conventionally used in continunous
784 : : //! Galerkin finite element methods (such as ALECG implements), whereas the
785 : : //! latter, in finite volume methods.
786 : : // *****************************************************************************
787 : : {
788 : 47340 : const auto& coord = Disc()->Coord();
789 : :
790 : 47340 : conserved( m_u, Disc()->Vol() );
791 : :
792 : : // Apply Dirichlet BCs
793 [ + + ]: 631884 : for (const auto& [b,bc] : m_dirbc)
794 [ + + ]: 3507264 : for (ncomp_t c=0; c<m_u.nprop(); ++c)
795 [ + - ][ + - ]: 2922720 : if (bc[c].first) m_u(b,c,0) = bc[c].second;
796 : :
797 : : // Apply symmetry BCs
798 [ + + ]: 94680 : for (const auto& eq : g_cgpde)
799 [ + - ]: 47340 : eq.symbc( m_u, coord, m_bnorm, m_symbcnodes );
800 : :
801 : : // Apply farfield BCs
802 [ + + ]: 94680 : for (const auto& eq : g_cgpde)
803 [ + - ]: 47340 : eq.farfieldbc( m_u, coord, m_bnorm, m_farfieldbcnodes );
804 : :
805 : : // Apply sponge conditions
806 [ + + ]: 94680 : for (const auto& eq : g_cgpde)
807 [ + - ]: 47340 : eq.sponge( m_u, coord, m_spongenodes );
808 : :
809 : : // Apply user defined time dependent BCs
810 [ + + ]: 94680 : for (const auto& eq : g_cgpde)
811 [ + - ][ + - ]: 47340 : eq.timedepbc( Disc()->T(), m_u, m_timedepbcnodes, m_timedepbcFn );
812 : :
813 : 47340 : volumetric( m_u, Disc()->Vol() );
814 : 47340 : }
815 : :
816 : : void
817 : 15602 : ALECG::next()
818 : : // *****************************************************************************
819 : : // Continue to next time step
820 : : // *****************************************************************************
821 : : {
822 : 15602 : dt();
823 : 15602 : }
824 : :
825 : : void
826 : 15602 : ALECG::dt()
827 : : // *****************************************************************************
828 : : // Compute time step size
829 : : // *****************************************************************************
830 : : {
831 : 15602 : tk::real mindt = std::numeric_limits< tk::real >::max();
832 : :
833 : 15602 : auto const_dt = g_inputdeck.get< tag::discr, tag::dt >();
834 : 15602 : auto def_const_dt = g_inputdeck_defaults.get< tag::discr, tag::dt >();
835 : 15602 : auto eps = std::numeric_limits< tk::real >::epsilon();
836 : :
837 [ + - ]: 15602 : auto d = Disc();
838 : :
839 : : // use constant dt if configured
840 [ + + ]: 15602 : if (std::abs(const_dt - def_const_dt) > eps) {
841 : :
842 : 7290 : mindt = const_dt;
843 : :
844 : : } else { // compute dt based on CFL
845 : :
846 : : //! [Find the minimum dt across all PDEs integrated]
847 [ + - ][ + - ]: 8312 : conserved( m_u, Disc()->Vol() );
848 [ + + ]: 8312 : if (g_inputdeck.get< tag::discr, tag::steady_state >()) {
849 : :
850 : : // compute new dt for each mesh point
851 [ + + ]: 400 : for (const auto& eq : g_cgpde)
852 [ + - ]: 200 : eq.dt( d->It(), d->Vol(), m_u, m_dtp );
853 : :
854 : : // find the smallest dt of all nodes on this chare
855 [ + - ]: 200 : mindt = *std::min_element( begin(m_dtp), end(m_dtp) );
856 : :
857 : : } else { // compute new dt for this chare
858 : :
859 : : // find the smallest dt of all equations on this chare
860 [ + + ]: 16224 : for (const auto& eq : g_cgpde) {
861 [ + - ]: 8112 : auto eqdt = eq.dt( d->Coord(), d->Inpoel(), d->T(), d->Dtn(), m_u,
862 : 8112 : d->Vol(), d->Voln() );
863 [ + - ]: 8112 : if (eqdt < mindt) mindt = eqdt;
864 : : }
865 : :
866 : : }
867 [ + - ][ + - ]: 8312 : volumetric( m_u, Disc()->Vol() );
868 : : //! [Find the minimum dt across all PDEs integrated]
869 : :
870 : : }
871 : :
872 : : //! [Advance]
873 : : // Actiavate SDAG waits for next time step stage
874 [ + - ][ + - ]: 15602 : thisProxy[ thisIndex ].wait4grad();
875 [ + - ][ + - ]: 15602 : thisProxy[ thisIndex ].wait4rhs();
876 : :
877 : : // Contribute to minimum dt across all chares and advance to next step
878 [ + - ]: 15602 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
879 [ + - ][ + - ]: 31204 : CkCallback(CkReductionTarget(ALECG,advance), thisProxy) );
880 : : //! [Advance]
881 : 15602 : }
882 : :
883 : : void
884 : 15602 : ALECG::advance( tk::real newdt )
885 : : // *****************************************************************************
886 : : // Advance equations to next time step
887 : : //! \param[in] newdt The smallest dt across the whole problem
888 : : // *****************************************************************************
889 : : {
890 : 15602 : auto d = Disc();
891 : :
892 : : // Set new time step size
893 [ + - ]: 15602 : if (m_stage == 0) d->setdt( newdt );
894 : :
895 : : // Compute gradients for next time step
896 : 15602 : chBndGrad();
897 : 15602 : }
898 : :
899 : : void
900 : 46806 : ALECG::chBndGrad()
901 : : // *****************************************************************************
902 : : // Compute nodal gradients at chare-boundary nodes. Gradients at internal nodes
903 : : // are calculated locally as needed and are not stored.
904 : : // *****************************************************************************
905 : : {
906 : 46806 : auto d = Disc();
907 : :
908 : : // Divide solution with mesh volume
909 : 46806 : conserved( m_u, Disc()->Vol() );
910 : : // Compute own portion of gradients for all equations
911 [ + + ]: 93612 : for (const auto& eq : g_cgpde)
912 [ + - ]: 46806 : eq.chBndGrad( d->Coord(), d->Inpoel(), m_bndel, d->Gid(), d->Bid(), m_u,
913 : 46806 : m_chBndGrad );
914 : : // Multiply solution with mesh volume
915 : 46806 : volumetric( m_u, Disc()->Vol() );
916 : :
917 : : // Communicate gradients to other chares on chare-boundary
918 [ + + ]: 46806 : if (d->NodeCommMap().empty()) // in serial we are done
919 : 1512 : comgrad_complete();
920 : : else // send gradients contributions to chare-boundary nodes to fellow chares
921 [ + + ]: 506070 : for (const auto& [c,n] : d->NodeCommMap()) {
922 [ + - ]: 460776 : std::vector< std::vector< tk::real > > g( n.size() );
923 : 460776 : std::size_t j = 0;
924 [ + + ][ + - ]: 7779774 : for (auto i : n) g[ j++ ] = m_chBndGrad[ tk::cref_find(d->Bid(),i) ];
[ + - ]
925 [ + - ][ + - ]: 460776 : thisProxy[c].comChBndGrad( std::vector<std::size_t>(begin(n),end(n)), g );
[ + - ]
926 : : }
927 : :
928 : 46806 : owngrad_complete();
929 : 46806 : }
930 : :
931 : : void
932 : 460776 : ALECG::comChBndGrad( const std::vector< std::size_t >& gid,
933 : : const std::vector< std::vector< tk::real > >& G )
934 : : // *****************************************************************************
935 : : // Receive contributions to nodal gradients on chare-boundaries
936 : : //! \param[in] gid Global mesh node IDs at which we receive grad contributions
937 : : //! \param[in] G Partial contributions of gradients to chare-boundary nodes
938 : : //! \details This function receives contributions to m_chBndGrad, which stores
939 : : //! nodal gradients at mesh chare-boundary nodes. While m_chBndGrad stores
940 : : //! own contributions, m_chBndGradc collects the neighbor chare
941 : : //! contributions during communication. This way work on m_chBndGrad and
942 : : //! m_chBndGradc is overlapped. The two are combined in rhs().
943 : : // *****************************************************************************
944 : : {
945 [ - + ][ - - ]: 460776 : Assert( G.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
946 : :
947 : : using tk::operator+=;
948 : :
949 [ + + ]: 7779774 : for (std::size_t i=0; i<gid.size(); ++i) m_chBndGradc[ gid[i] ] += G[i];
950 : :
951 [ + + ]: 460776 : if (++m_ngrad == Disc()->NodeCommMap().size()) {
952 : 45294 : m_ngrad = 0;
953 : 45294 : comgrad_complete();
954 : : }
955 : 460776 : }
956 : :
957 : : void
958 : 46806 : ALECG::rhs()
959 : : // *****************************************************************************
960 : : // Compute right-hand side of transport equations
961 : : // *****************************************************************************
962 : : {
963 : 46806 : auto d = Disc();
964 : :
965 : : // Combine own and communicated contributions to nodal gradients
966 [ + + ]: 4327527 : for (const auto& [gid,g] : m_chBndGradc) {
967 [ + - ]: 4280721 : auto bid = tk::cref_find( d->Bid(), gid );
968 [ + + ]: 35617056 : for (ncomp_t c=0; c<m_chBndGrad.nprop(); ++c)
969 [ + - ]: 31336335 : m_chBndGrad(bid,c,0) += g[c];
970 : : }
971 : :
972 : : // clear gradients receive buffer
973 : 46806 : tk::destroy(m_chBndGradc);
974 : :
975 : 46806 : const auto steady = g_inputdeck.get< tag::discr, tag::steady_state >();
976 : :
977 : : // Compute own portion of right-hand side for all equations
978 [ + + ]: 46806 : auto prev_rkcoef = m_stage == 0 ? 0.0 : rkcoef[m_stage-1];
979 [ + + ]: 46806 : if (steady)
980 [ + + ]: 67560 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] += prev_rkcoef * m_dtp[p];
981 : 46806 : conserved( m_u, Disc()->Vol() );
982 [ + + ]: 93612 : for (const auto& eq : g_cgpde) {
983 [ + - ]: 93612 : eq.rhs( d->T() + prev_rkcoef * d->Dt(), d->Coord(), d->Inpoel(),
984 : 46806 : m_triinpoel, d->Gid(), d->Bid(), d->Lid(), m_dfn, m_psup, m_esup,
985 : 46806 : m_symbctri, m_spongenodes, d->Vol(), m_edgenode, m_edgeid,
986 [ + - ]: 46806 : m_boxnodes, m_chBndGrad, m_u, d->meshvel(), m_tp, d->Boxvol(),
987 : 46806 : m_rhs );
988 : : }
989 : 46806 : volumetric( m_u, Disc()->Vol() );
990 [ + + ]: 46806 : if (steady)
991 [ + + ]: 67560 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] -= prev_rkcoef * m_dtp[p];
992 : :
993 : : // Query/update boundary-conditions-related data structures from user input
994 : 46806 : queryBnd();
995 : :
996 : : // Communicate rhs to other chares on chare-boundary
997 [ + + ]: 46806 : if (d->NodeCommMap().empty()) // in serial we are done
998 : 1512 : comrhs_complete();
999 : : else // send contributions of rhs to chare-boundary nodes to fellow chares
1000 [ + + ]: 506070 : for (const auto& [c,n] : d->NodeCommMap()) {
1001 [ + - ]: 460776 : std::vector< std::vector< tk::real > > r( n.size() );
1002 : 460776 : std::size_t j = 0;
1003 [ + + ][ + - ]: 7779774 : for (auto i : n) r[ j++ ] = m_rhs[ tk::cref_find(d->Lid(),i) ];
[ + - ]
1004 [ + - ][ + - ]: 460776 : thisProxy[c].comrhs( std::vector<std::size_t>(begin(n),end(n)), r );
[ + - ]
1005 : : }
1006 : :
1007 : 46806 : ownrhs_complete();
1008 : 46806 : }
1009 : :
1010 : : void
1011 : 460776 : ALECG::comrhs( const std::vector< std::size_t >& gid,
1012 : : const std::vector< std::vector< tk::real > >& R )
1013 : : // *****************************************************************************
1014 : : // Receive contributions to right-hand side vector on chare-boundaries
1015 : : //! \param[in] gid Global mesh node IDs at which we receive RHS contributions
1016 : : //! \param[in] R Partial contributions of RHS to chare-boundary nodes
1017 : : //! \details This function receives contributions to m_rhs, which stores the
1018 : : //! right hand side vector at mesh nodes. While m_rhs stores own
1019 : : //! contributions, m_rhsc collects the neighbor chare contributions during
1020 : : //! communication. This way work on m_rhs and m_rhsc is overlapped. The two
1021 : : //! are combined in solve().
1022 : : // *****************************************************************************
1023 : : {
1024 [ - + ][ - - ]: 460776 : Assert( R.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
1025 : :
1026 : : using tk::operator+=;
1027 : :
1028 [ + + ]: 7779774 : for (std::size_t i=0; i<gid.size(); ++i) m_rhsc[ gid[i] ] += R[i];
1029 : :
1030 : : // When we have heard from all chares we communicate with, this chare is done
1031 [ + + ]: 460776 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
1032 : 45294 : m_nrhs = 0;
1033 : 45294 : comrhs_complete();
1034 : : }
1035 : 460776 : }
1036 : :
1037 : : void
1038 : 46806 : ALECG::solve()
1039 : : // *****************************************************************************
1040 : : // Advance systems of equations
1041 : : // *****************************************************************************
1042 : : {
1043 : 46806 : auto d = Disc();
1044 : :
1045 : : // Combine own and communicated contributions to rhs
1046 [ + + ]: 4327527 : for (const auto& b : m_rhsc) {
1047 [ + - ]: 4280721 : auto lid = tk::cref_find( d->Lid(), b.first );
1048 [ + + ][ + - ]: 14726166 : for (ncomp_t c=0; c<m_rhs.nprop(); ++c) m_rhs(lid,c,0) += b.second[c];
1049 : : }
1050 : :
1051 : : // clear receive buffer
1052 : 46806 : tk::destroy(m_rhsc);
1053 : :
1054 : : // Update state at time n
1055 [ + + ]: 46806 : if (m_stage == 0) {
1056 : 15602 : m_un = m_u;
1057 [ + + ]: 15602 : if (g_inputdeck.get< tag::ale, tag::ale >()) d->UpdateCoordn();
1058 : : }
1059 : :
1060 : : // Solve the sytem
1061 [ + + ]: 46806 : if (g_inputdeck.get< tag::discr, tag::steady_state >()) {
1062 : :
1063 : : // Advance solution, converging to steady state
1064 [ + + ]: 67560 : for (std::size_t i=0; i<m_u.nunk(); ++i)
1065 [ + + ]: 401760 : for (ncomp_t c=0; c<m_u.nprop(); ++c)
1066 : 334800 : m_u(i,c,0) = m_un(i,c,0) + rkcoef[m_stage] * m_dtp[i] * m_rhs(i,c,0);
1067 : :
1068 : : } else {
1069 : :
1070 : 46206 : auto adt = rkcoef[m_stage] * d->Dt();
1071 : :
1072 : : // Advance unsteady solution
1073 [ + - ]: 46206 : m_u = m_un + adt * m_rhs;
1074 : :
1075 : : // Advance mesh if ALE is enabled
1076 [ + + ]: 46206 : if (g_inputdeck.get< tag::ale, tag::ale >()) {
1077 : 1080 : auto& coord = d->Coord();
1078 : 1080 : const auto& w = d->meshvel();
1079 [ + + ]: 3480 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
1080 [ + + ]: 1828620 : for (std::size_t i=0; i<coord[j].size(); ++i)
1081 [ + - ]: 1826220 : coord[j][i] = d->Coordn()[j][i] + adt * w(i,j,0);
1082 : : }
1083 : :
1084 : : }
1085 : :
1086 : 46806 : m_newmesh = 0; // recompute normals after ALE (if enabled)
1087 : : // Activate SDAG waits
1088 [ + - ]: 46806 : thisProxy[ thisIndex ].wait4norm();
1089 [ + - ]: 46806 : thisProxy[ thisIndex ].wait4mesh();
1090 : :
1091 : : //! [Continue after solve]
1092 : : // Recompute mesh volumes if ALE is enabled
1093 [ + + ]: 46806 : if (g_inputdeck.get< tag::ale, tag::ale >()) {
1094 : :
1095 [ + - ]: 1080 : transfer_complete();
1096 : : // Save nodal volumes at previous time step stage
1097 [ + - ]: 1080 : d->Voln() = d->Vol();
1098 : : // Prepare for recomputing the nodal volumes
1099 [ + - ]: 1080 : d->startvol();
1100 : 1080 : auto meshid = d->MeshId();
1101 [ + - ]: 1080 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1102 [ + - ][ + - ]: 2160 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1103 : :
1104 : : } else {
1105 : :
1106 : 45726 : norm_complete();
1107 : 45726 : resized();
1108 : :
1109 : : }
1110 : : //! [Continue after solve]
1111 : 46806 : }
1112 : :
1113 : : void
1114 : 46806 : ALECG::ale()
1115 : : // *****************************************************************************
1116 : : // Continue after computing the new mesh velocity for ALE
1117 : : // *****************************************************************************
1118 : : {
1119 : 46806 : auto d = Disc();
1120 : :
1121 [ + + ]: 46806 : if (m_stage < 2) {
1122 : :
1123 : : // Activate SDAG wait for next time step stage
1124 [ + - ]: 31204 : thisProxy[ thisIndex ].wait4grad();
1125 [ + - ]: 31204 : thisProxy[ thisIndex ].wait4rhs();
1126 : :
1127 : : // continue to mesh-to-mesh transfer (if coupled)
1128 : 31204 : transfer();
1129 : :
1130 : : } else {
1131 : :
1132 : : // Ensure new field output file if mesh moved if ALE is enabled
1133 [ + + ]: 15602 : if (g_inputdeck.get< tag::ale, tag::ale >()) {
1134 : 360 : d->Itf() = 0; // Zero field output iteration count if mesh moved
1135 : 360 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1136 : : }
1137 : :
1138 : : // Compute diagnostics, e.g., residuals
1139 : 15602 : conserved( m_u, Disc()->Vol() );
1140 : 15602 : conserved( m_un, Disc()->Voln() );
1141 : 31204 : auto diag_computed = m_diag.compute( *d, m_u, m_un, m_bnorm,
1142 : 15602 : m_symbcnodes, m_farfieldbcnodes );
1143 : 15602 : volumetric( m_u, Disc()->Vol() );
1144 : 15602 : volumetric( m_un, Disc()->Voln() );
1145 : : // Increase number of iterations and physical time
1146 : 15602 : d->next();
1147 : : // Advance physical time for local time stepping
1148 [ + + ]: 15602 : if (g_inputdeck.get< tag::discr, tag::steady_state >())
1149 [ + + ]: 22520 : for (std::size_t i=0; i<m_u.nunk(); ++i) m_tp[i] += m_dtp[i];
1150 : : // Continue to mesh refinement (if configured)
1151 [ + + ][ + - ]: 15602 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ + - ]
1152 : :
1153 : : }
1154 : 46806 : }
1155 : :
1156 : : //! [Refine]
1157 : : void
1158 : 15602 : ALECG::refine( const std::vector< tk::real >& l2res )
1159 : : // *****************************************************************************
1160 : : // Optionally refine/derefine mesh
1161 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1162 : : //! computed across the whole problem
1163 : : // *****************************************************************************
1164 : : {
1165 : 15602 : auto d = Disc();
1166 : :
1167 : 15602 : const auto steady = g_inputdeck.get< tag::discr, tag::steady_state >();
1168 : 15602 : const auto residual = g_inputdeck.get< tag::discr, tag::residual >();
1169 : 15602 : const auto rc = g_inputdeck.get< tag::discr, tag::rescomp >() - 1;
1170 : :
1171 [ + + ]: 15602 : if (steady) {
1172 : :
1173 : : // this is the last time step if max time of max number of time steps
1174 : : // reached or the residual has reached its convergence criterion
1175 [ + - ][ + + ]: 200 : if (d->finished() or l2res[rc] < residual) m_finished = 1;
[ + + ]
1176 : :
1177 : : } else {
1178 : :
1179 : : // this is the last time step if max time or max iterations reached
1180 [ + + ]: 15402 : if (d->finished()) m_finished = 1;
1181 : :
1182 : : }
1183 : :
1184 : 15602 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
1185 : 15602 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
1186 : :
1187 : : // Activate SDAG waits for re-computing the normals
1188 : 15602 : m_newmesh = 1; // recompute normals after AMR (if enabled)
1189 [ + - ]: 15602 : thisProxy[ thisIndex ].wait4norm();
1190 [ + - ]: 15602 : thisProxy[ thisIndex ].wait4mesh();
1191 : :
1192 : : // if t>0 refinement enabled and we hit the frequency
1193 [ - + ][ - - ]: 15602 : if (dtref && !(d->It() % dtfreq)) { // refine
[ - + ]
1194 : :
1195 : 0 : d->startvol();
1196 [ - - ]: 0 : d->Ref()->dtref( {}, m_bnode, {} );
1197 : 0 : d->refined() = 1;
1198 : :
1199 : : } else { // do not refine
1200 : :
1201 : 15602 : d->refined() = 0;
1202 : 15602 : norm_complete();
1203 : 15602 : resized();
1204 : :
1205 : : }
1206 : 15602 : }
1207 : : //! [Refine]
1208 : :
1209 : : //! [Resize]
1210 : : void
1211 : 0 : ALECG::resizePostAMR(
1212 : : const std::vector< std::size_t >& /*ginpoel*/,
1213 : : const tk::UnsMesh::Chunk& chunk,
1214 : : const tk::UnsMesh::Coords& coord,
1215 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1216 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1217 : : const std::set< std::size_t >& removedNodes,
1218 : : const tk::NodeCommMap& nodeCommMap,
1219 : : const std::map< int, std::vector< std::size_t > >& bface,
1220 : : const std::map< int, std::vector< std::size_t > >& bnode,
1221 : : const std::vector< std::size_t >& triinpoel )
1222 : : // *****************************************************************************
1223 : : // Receive new mesh from Refiner
1224 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1225 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1226 : : //! \param[in] coord New mesh node coordinates
1227 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1228 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1229 : : //! \param[in] removedNodes Newly removed mesh nodes (local ids)
1230 : : //! \param[in] nodeCommMap New node communication map
1231 : : //! \param[in] bface Boundary-faces mapped to side set ids
1232 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1233 : : //! \param[in] triinpoel Boundary-face connectivity
1234 : : // *****************************************************************************
1235 : : {
1236 [ - - ]: 0 : auto d = Disc();
1237 : :
1238 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1239 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1240 : :
1241 : : // Resize mesh data structures after mesh refinement
1242 [ - - ]: 0 : d->resizePostAMR( chunk, coord, nodeCommMap );
1243 : :
1244 : : // Remove newly removed nodes from solution vectors
1245 [ - - ]: 0 : m_u.rm(removedNodes);
1246 [ - - ]: 0 : m_un.rm(removedNodes);
1247 [ - - ]: 0 : m_rhs.rm(removedNodes);
1248 : :
1249 : : // Resize auxiliary solution vectors
1250 : 0 : auto npoin = coord[0].size();
1251 : 0 : auto nprop = m_u.nprop();
1252 [ - - ]: 0 : m_u.resize( npoin );
1253 [ - - ]: 0 : m_un.resize( npoin );
1254 [ - - ]: 0 : m_rhs.resize( npoin );
1255 [ - - ]: 0 : m_chBndGrad.resize( d->Bid().size() );
1256 : :
1257 : : // Update solution on new mesh
1258 [ - - ]: 0 : for (const auto& n : addedNodes)
1259 [ - - ]: 0 : for (std::size_t c=0; c<nprop; ++c)
1260 [ - - ][ - - ]: 0 : m_u(n.first,c,0) = (m_u(n.second[0],c,0) + m_u(n.second[1],c,0))/2.0;
[ - - ]
1261 : :
1262 : : // Update physical-boundary node-, face-, and element lists
1263 [ - - ]: 0 : m_bnode = bnode;
1264 [ - - ]: 0 : m_bface = bface;
1265 [ - - ]: 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1266 : :
1267 : 0 : auto meshid = d->MeshId();
1268 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1269 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1270 : 0 : }
1271 : : //! [Resize]
1272 : :
1273 : : void
1274 : 62408 : ALECG::resized()
1275 : : // *****************************************************************************
1276 : : // Resizing data sutrctures after mesh refinement has been completed
1277 : : // *****************************************************************************
1278 : : {
1279 : 62408 : resize_complete();
1280 : 62408 : }
1281 : :
1282 : : void
1283 : 46806 : ALECG::transfer()
1284 : : // *****************************************************************************
1285 : : // Transfer solution to other solver and mesh if coupled
1286 : : // *****************************************************************************
1287 : : {
1288 : : // Initiate solution transfer (if coupled)
1289 : : //Disc()->transfer( m_u,
1290 : : // CkCallback(CkIndex_ALECG::stage(), thisProxy[thisIndex]) );
1291 [ + - ]: 46806 : thisProxy[thisIndex].stage();
1292 : 46806 : }
1293 : :
1294 : : //! [stage]
1295 : : void
1296 : 46806 : ALECG::stage()
1297 : : // *****************************************************************************
1298 : : // Evaluate whether to continue with next time step stage
1299 : : // *****************************************************************************
1300 : : {
1301 : : // Increment Runge-Kutta stage counter
1302 : 46806 : ++m_stage;
1303 : :
1304 : : // if not all Runge-Kutta stages complete, continue to next time stage,
1305 : : // otherwise output field data to file(s)
1306 [ + + ]: 46806 : if (m_stage < 3) chBndGrad(); else out();
1307 : 46806 : }
1308 : : //! [stage]
1309 : :
1310 : : void
1311 : 1994 : ALECG::writeFields( CkCallback c )
1312 : : // *****************************************************************************
1313 : : // Output mesh-based fields to file
1314 : : //! \param[in] c Function to continue with after the write
1315 : : // *****************************************************************************
1316 : : {
1317 [ + + ]: 1994 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) {
1318 : :
1319 : 608 : c.send();
1320 : :
1321 : : } else {
1322 : :
1323 [ + - ]: 1386 : auto d = Disc();
1324 : 1386 : const auto& coord = d->Coord();
1325 : :
1326 : : // Query fields names requested by user
1327 [ + - ]: 2772 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
1328 : : // Collect field output from numerical solution requested by user
1329 [ + - ][ + - ]: 1386 : conserved( m_u, Disc()->Vol() );
1330 [ + - ]: 2772 : auto nodefields = numericFieldOutput( m_u, tk::Centering::NODE );
1331 [ + - ][ + - ]: 1386 : volumetric( m_u, Disc()->Vol() );
1332 : :
1333 : : //! Lambda to put in a field for output if not empty
1334 : 300 : auto add_node_field = [&]( const auto& name, const auto& field ){
1335 [ + - ]: 300 : if (not field.empty()) {
1336 [ + - ][ + - ]: 300 : nodefieldnames.push_back( name );
1337 : 300 : nodefields.push_back( field );
1338 : : }
1339 : 1686 : };
1340 : :
1341 : : // Output mesh velocity if ALE is enabled
1342 [ + + ]: 1386 : if (g_inputdeck.get< tag::ale, tag::ale >()) {
1343 [ + - ]: 75 : const auto& w = d->meshvel();
1344 [ + - ][ + - ]: 75 : add_node_field( "x-mesh-velocity", w.extract(0,0) );
1345 [ + - ][ + - ]: 75 : add_node_field( "y-mesh-velocity", w.extract(1,0) );
1346 [ + - ][ + - ]: 75 : add_node_field( "z-mesh-velocity", w.extract(2,0) );
1347 [ + - ]: 75 : add_node_field( "volume", d->Vol() );
1348 : : }
1349 : :
1350 : : // Collect field output names for analytical solutions
1351 [ + + ]: 2772 : for (const auto& eq : g_cgpde)
1352 [ + - ]: 1386 : analyticFieldNames( eq, tk::Centering::NODE, nodefieldnames );
1353 : :
1354 : : // Collect field output from analytical solutions (if exist)
1355 [ + + ]: 2772 : for (const auto& eq : g_cgpde)
1356 [ + - ]: 2772 : analyticFieldOutput( eq, tk::Centering::NODE, coord[0], coord[1],
1357 : 1386 : coord[2], d->T(), nodefields );
1358 : :
1359 : : // Query and collect block and surface field names from PDEs integrated
1360 : 2772 : std::vector< std::string > nodesurfnames;
1361 [ + + ]: 2772 : for (const auto& eq : g_cgpde) {
1362 [ + - ]: 1386 : auto s = eq.surfNames();
1363 [ + - ]: 1386 : nodesurfnames.insert( end(nodesurfnames), begin(s), end(s) );
1364 : : }
1365 : :
1366 : : // Collect node block and surface field solution
1367 : 1386 : std::vector< std::vector< tk::real > > nodesurfs;
1368 [ + - ][ + - ]: 1386 : conserved( m_u, Disc()->Vol() );
1369 [ + + ]: 2772 : for (const auto& eq : g_cgpde) {
1370 [ + - ][ + - ]: 1386 : auto s = eq.surfOutput( tk::bfacenodes(m_bface,m_triinpoel), m_u );
1371 [ + - ]: 1386 : nodesurfs.insert( end(nodesurfs), begin(s), end(s) );
1372 : : }
1373 [ + - ][ + - ]: 1386 : volumetric( m_u, Disc()->Vol() );
1374 : :
1375 [ - + ][ - - ]: 1386 : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
1376 : :
1377 : : // Send mesh and fields data (solution dump) for output to file
1378 [ + - ][ + - ]: 2772 : d->write( d->Inpoel(), coord, m_bface, tk::remap(m_bnode,d->Lid()),
1379 : 1386 : m_triinpoel, {}, nodefieldnames, nodesurfnames, {}, nodefields,
1380 : : nodesurfs, c );
1381 : :
1382 : : }
1383 : 1994 : }
1384 : :
1385 : : void
1386 : 15602 : ALECG::out()
1387 : : // *****************************************************************************
1388 : : // Output mesh field data
1389 : : // *****************************************************************************
1390 : : {
1391 : 15602 : auto d = Disc();
1392 : :
1393 : : // Output time history
1394 [ + + ][ + + ]: 15602 : if (d->histiter() or d->histtime() or d->histrange()) {
[ + + ][ + + ]
1395 : 1796 : std::vector< std::vector< tk::real > > hist;
1396 [ + - ][ + - ]: 898 : conserved( m_u, Disc()->Vol() );
1397 [ + + ]: 1796 : for (const auto& eq : g_cgpde) {
1398 [ + - ]: 898 : auto h = eq.histOutput( d->Hist(), d->Inpoel(), m_u );
1399 [ + - ]: 898 : hist.insert( end(hist), begin(h), end(h) );
1400 : : }
1401 [ + - ][ + - ]: 898 : volumetric( m_u, Disc()->Vol() );
1402 [ + - ]: 898 : d->history( std::move(hist) );
1403 : : }
1404 : :
1405 : : // Output field data
1406 [ + + ][ + + ]: 15602 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished)
[ + - ][ + + ]
[ + + ]
1407 [ + - ][ + - ]: 1460 : writeFields( CkCallback(CkIndex_ALECG::step(), thisProxy[thisIndex]) );
[ + - ]
1408 : : else
1409 : 14142 : step();
1410 : 15602 : }
1411 : :
1412 : : void
1413 : 15068 : ALECG::evalLB( int nrestart )
1414 : : // *****************************************************************************
1415 : : // Evaluate whether to do load balancing
1416 : : //! \param[in] nrestart Number of times restarted
1417 : : // *****************************************************************************
1418 : : {
1419 : 15068 : auto d = Disc();
1420 : :
1421 : : // Detect if just returned from a checkpoint and if so, zero timers and
1422 : : // finished flag
1423 [ - + ]: 15068 : if (d->restarted( nrestart )) m_finished = 0;
1424 : :
1425 : 15068 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1426 : 15068 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
1427 : :
1428 : : // Load balancing if user frequency is reached or after the second time-step
1429 [ + + ][ + + ]: 15068 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
[ + + ]
1430 : :
1431 : 10306 : AtSync();
1432 [ - + ]: 10306 : if (nonblocking) next();
1433 : :
1434 : : } else {
1435 : :
1436 : 4762 : next();
1437 : :
1438 : : }
1439 : 15068 : }
1440 : :
1441 : : void
1442 : 15068 : ALECG::evalRestart()
1443 : : // *****************************************************************************
1444 : : // Evaluate whether to save checkpoint/restart
1445 : : // *****************************************************************************
1446 : : {
1447 : 15068 : auto d = Disc();
1448 : :
1449 : 15068 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
1450 : 15068 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1451 : :
1452 [ + + ][ - + ]: 15068 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
1453 : :
1454 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1455 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1456 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
1457 : :
1458 : : } else {
1459 : :
1460 : 15068 : evalLB( /* nrestart = */ -1 );
1461 : :
1462 : : }
1463 : 15068 : }
1464 : :
1465 : : void
1466 : 15602 : ALECG::step()
1467 : : // *****************************************************************************
1468 : : // Evaluate whether to continue with next time step
1469 : : // *****************************************************************************
1470 : : {
1471 : 15602 : auto d = Disc();
1472 : :
1473 : : // Output one-liner status report to screen
1474 : 15602 : d->status();
1475 : : // Reset Runge-Kutta stage counter
1476 : 15602 : m_stage = 0;
1477 : :
1478 [ + + ]: 15602 : if (not m_finished) {
1479 : :
1480 : 15068 : evalRestart();
1481 : :
1482 : : } else {
1483 : :
1484 : 534 : auto meshid = d->MeshId();
1485 [ + - ]: 534 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1486 [ + - ][ + - ]: 1068 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1487 : :
1488 : : }
1489 : 15602 : }
1490 : :
1491 : : #include "NoWarning/alecg.def.h"
|
