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