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