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