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