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