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