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