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