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