Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ALE.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 Definitions file for distributed ALE mesh motion
9 : : \details Definitions file for asynchronous distributed
10 : : arbitrary Lagrangian-Eulerian (ALE) mesh motion using Charm++.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include "ALE.hpp"
15 : : #include "DerivedData.hpp"
16 : : #include "Vector.hpp"
17 : : #include "Inciter/InputDeck/InputDeck.hpp"
18 : :
19 : : namespace inciter {
20 : :
21 : : extern ctr::InputDeck g_inputdeck;
22 : :
23 : : } // inciter::
24 : :
25 : : using inciter::ALE;
26 : :
27 : 31 : ALE::ALE( const tk::CProxy_ConjugateGradients& conjugategradientsproxy,
28 : : const std::array< std::vector< tk::real >, 3 >& coord,
29 : : const std::vector< std::size_t >& inpoel,
30 : : const std::vector< std::size_t >& gid,
31 : : const std::unordered_map< std::size_t, std::size_t >& lid,
32 : 31 : const tk::NodeCommMap& nodeCommMap ) :
33 : : m_conjugategradients( conjugategradientsproxy ),
34 : : m_done(),
35 : : m_nvort( 0 ),
36 : : m_ndiv( 0 ),
37 : : m_npot( 0 ),
38 : : m_nwf( 0 ),
39 : : m_nodeCommMap( nodeCommMap ),
40 : : m_lid( lid ),
41 : : m_coord0( coord ),
42 : : m_inpoel( inpoel ),
43 : : m_vol0(),
44 : : m_vol(),
45 : : m_it( 0 ),
46 : : m_t( 0.0 ),
47 : : m_w( gid.size(), 3 ),
48 : : m_wf( gid.size(), 3 ),
49 : : m_wfc(),
50 : : m_veldiv(),
51 : : m_veldivc(),
52 : : m_gradpot(),
53 : : m_gradpotc(),
54 : : m_vorticity(),
55 : : m_vorticityc(),
56 : : m_meshveldirbcnodes(),
57 : : m_meshvelsymbcnodes(),
58 [ + - ][ + - ]: 124 : m_move( moveCfg() )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
59 : : // *****************************************************************************
60 : : // Constructor
61 : : //! \param[in] conjugategradientsproxy Distributed Conjugrate Gradients linear
62 : : //! solver proxy (bound to ALE proxy)
63 : : //! \param[in] coord Mesh node coordinates
64 : : //! \param[in] inpoel Mesh element connectivity
65 : : //! \param[in] gid Local->global node id map
66 : : //! \param[in] lid Global->locbal node id map
67 : : //! \param[in] nodeCommMap Node communication map
68 : : // *****************************************************************************
69 : : {
70 : 31 : auto smoother = g_inputdeck.get< tag::ale, tag::smoother >();
71 : :
72 [ + + ]: 31 : if (smoother == ctr::MeshVelocitySmootherType::LAPLACE)
73 : 17 : m_conjugategradients[ thisIndex ]. // solve for mesh velocity
74 [ + - ][ + - ]: 34 : insert( Laplacian( 3, coord ), gid, m_lid, m_nodeCommMap );
[ + - ][ - + ]
[ - - ]
75 [ + + ]: 14 : else if (smoother == ctr::MeshVelocitySmootherType::HELMHOLTZ)
76 : 4 : m_conjugategradients[ thisIndex ]. // solve for scalar potential
77 [ + - ][ + - ]: 8 : insert( Laplacian( 1, coord ), gid, m_lid, m_nodeCommMap );
[ + - ][ - + ]
[ - - ]
78 : :
79 : : // Zero ALE mesh velocity
80 : : m_w.fill( 0.0 );
81 : :
82 : : // Activate SDAG wait for initially computing prerequisites for ALE
83 [ + - ][ + - ]: 31 : thisProxy[ thisIndex ].wait4vel();
84 [ + - ][ + - ]: 31 : thisProxy[ thisIndex ].wait4pot();
85 [ + - ][ + - ]: 31 : thisProxy[ thisIndex ].wait4for();
86 : 31 : }
87 : :
88 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
89 : 21 : ALE::Laplacian( std::size_t ncomp,
90 : : const std::array< std::vector< tk::real >, 3 >& coord ) const
91 : : // *****************************************************************************
92 : : // Generate {A,x,b} for Laplacian mesh velocity smoother
93 : : //! \param[in] ncomp Number of scalar components
94 : : //! \param[in] coord Mesh node coordinates
95 : : //! \return {A,x,b} with a Laplacian, unknown, and rhs initialized with zeros
96 : : //! \see Waltz, et al. "A three-dimensional finite element arbitrary
97 : : //! Lagrangian-Eulerian method for shock hydrodynamics on unstructured
98 : : //! grids", Computers& Fluids, 2013, and Bakosi, et al. "Improved ALE mesh
99 : : //! velocities for complex flows, International Journal for Numerical Methods
100 : : //! in Fluids, 2017.
101 : : // *****************************************************************************
102 : : {
103 [ + - ][ + - ]: 63 : tk::CSR A( ncomp, tk::genPsup(m_inpoel,4,tk::genEsup(m_inpoel,4)) );
104 : :
105 : : const auto& X = coord[0];
106 : : const auto& Y = coord[1];
107 : : const auto& Z = coord[2];
108 : :
109 : : // fill matrix with Laplacian
110 [ + + ]: 93293 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
111 : : // access node IDs
112 : : const std::array< std::size_t, 4 >
113 : 93272 : N{{ m_inpoel[e*4+0], m_inpoel[e*4+1], m_inpoel[e*4+2], m_inpoel[e*4+3] }};
114 : : // compute element Jacobi determinant
115 : : const std::array< tk::real, 3 >
116 : 93272 : ba{{ X[N[1]]-X[N[0]], Y[N[1]]-Y[N[0]], Z[N[1]]-Z[N[0]] }},
117 : 93272 : ca{{ X[N[2]]-X[N[0]], Y[N[2]]-Y[N[0]], Z[N[2]]-Z[N[0]] }},
118 : 93272 : da{{ X[N[3]]-X[N[0]], Y[N[3]]-Y[N[0]], Z[N[3]]-Z[N[0]] }};
119 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
120 : : Assert( J > 0, "Element Jacobian non-positive" );
121 : :
122 : : // shape function derivatives, nnode*ndim [4][3]
123 : : std::array< std::array< tk::real, 3 >, 4 > grad;
124 : : grad[1] = tk::crossdiv( ca, da, J );
125 : 93272 : grad[2] = tk::crossdiv( da, ba, J );
126 : 93272 : grad[3] = tk::crossdiv( ba, ca, J );
127 [ + + ]: 373088 : for (std::size_t i=0; i<3; ++i)
128 : 279816 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
129 : :
130 [ + + ]: 466360 : for (std::size_t a=0; a<4; ++a)
131 [ + + ]: 1492352 : for (std::size_t k=0; k<3; ++k)
132 [ + + ]: 5596320 : for (std::size_t b=0; b<4; ++b)
133 [ + + ]: 17838144 : for (std::size_t i=0; i<ncomp; ++i)
134 [ + - ]: 13361088 : A(N[a],N[b],i) -= J/6 * grad[a][k] * grad[b][k];
135 : : }
136 : :
137 [ + - ]: 21 : auto npoin = coord[0].size();
138 [ + - ][ + - ]: 21 : std::vector< tk::real > b(npoin*ncomp,0.0), x(npoin*ncomp,0.0);
[ - - ]
139 : :
140 : 21 : return { std::move(A), std::move(x), std::move(b) };
141 : : }
142 : :
143 : : decltype(ALE::m_move)
144 : 31 : ALE::moveCfg()
145 : : // *****************************************************************************
146 : : // Initialize user-defined functions for ALE moving sides
147 : : //! \details This function fills in only part of the data structure
148 : : //! returned, containing the user-defined functions in discrete form that will
149 : : //! be sampled in time. The node lists will be initialized later.
150 : : // *****************************************************************************
151 : : {
152 : : decltype(m_move) cfg;
153 : :
154 [ + + ]: 39 : for (const auto& m : g_inputdeck.get< tag::ale, tag::move >()) {
155 : : const auto& fn = m.get< tag::fn >();
156 : : Assert( fn.size() % 4 == 0, "Incomplete user-defined function" );
157 [ + - ]: 8 : cfg.emplace_back();
158 : : // store user-defined function type
159 : 8 : std::get<0>(cfg.back()) = m.get< tag::fntype >();
160 : : // store user-defined function discrete data
161 [ + + ]: 1204 : for (std::size_t i=0; i<fn.size()/4; ++i)
162 : : std::get<1>(cfg.back()).
163 [ + - ]: 1196 : push_back( {{ fn[i*4+0], fn[i*4+1], fn[i*4+2], fn[i*4+3] }} );
164 : : }
165 : :
166 : 31 : return cfg;
167 : : }
168 : :
169 : : void
170 : 1111 : ALE::meshvelBnd(
171 : : const std::map< int, std::vector< std::size_t > >& bface,
172 : : const std::map< int, std::vector< std::size_t > >& bnode,
173 : : const std::vector< std::size_t >& triinpoel )
174 : : // *****************************************************************************
175 : : // Query mesh velocity boundary conditions node lists and node list at which
176 : : // ALE moves boundaries
177 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
178 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
179 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set
180 : : // *****************************************************************************
181 : : {
182 : : // Prepare unique set of mesh velocity Dirichlet BC nodes
183 : 1111 : tk::destroy( m_meshveldirbcnodes );
184 : : std::unordered_map<int, std::unordered_set< std::size_t >> meshveldirbcnodes;
185 [ + + ]: 3333 : for (const auto& s : g_inputdeck.template get< tag::ale, tag::dirichlet >()) {
186 : 2222 : auto k = bface.find(static_cast<int>(s));
187 [ + + ]: 2222 : if (k != end(bface)) {
188 [ + - ]: 1054 : auto& n = meshveldirbcnodes[ k->first ]; // associate set id
189 [ + + ]: 65844 : for (auto f : k->second) { // face ids on side set
190 [ + - ]: 64790 : n.insert( triinpoel[f*3+0] );
191 [ + - ]: 64790 : n.insert( triinpoel[f*3+1] );
192 [ + - ]: 64790 : n.insert( triinpoel[f*3+2] );
193 : : }
194 : : }
195 : : }
196 [ + + ]: 2165 : for (const auto& [s,nodes] : meshveldirbcnodes)
197 : : m_meshveldirbcnodes.insert( begin(nodes), end(nodes) );
198 : :
199 : : // Prepare unique set of mesh velocity symmetry BC nodes. Note that somewhat
200 : : // counter-intuitively, we interrogate the boundary nodes instead of boundary
201 : : // faces here. This is because if we query the boundary faces, then we will
202 : : // get the mathematically correctly defined finite discrete surfaces
203 : : // (triangles) where mesh velocity symmetry BCs are configured by the user.
204 : : // However, in parallel, decomposing the domain and the boundary in various
205 : : // ways can produce situations on the boundary where boundary nodes are part
206 : : // of the given side set for mesh velocity symmetry BCs but not a full
207 : : // triangle face because, not all 3 nodes lie on the boundary. Thus
208 : : // interrogating the boundary nodes will be a superset and will include those
209 : : // nodes that are part of imposing symmetry BCs on nodes of faces that are
210 : : // only partial due to domain decomposition.
211 : 1111 : tk::destroy( m_meshvelsymbcnodes );
212 : : std::unordered_map<int, std::unordered_set< std::size_t >> meshvelsymbcnodes;
213 [ + + ]: 3343 : for (const auto& s : g_inputdeck.template get< tag::ale, tag::symmetry >()) {
214 : 2232 : auto k = bnode.find(static_cast<int>(s));
215 [ + + ]: 2232 : if (k != end(bnode)) {
216 [ + - ]: 1674 : auto& n = meshvelsymbcnodes[ k->first ]; // associate set id
217 [ + + ]: 918065 : for (auto g : k->second) { // node ids on side set
218 : : n.insert( tk::cref_find(m_lid,g) ); // store local ids
219 : : }
220 : : }
221 : : }
222 [ + + ]: 2785 : for (const auto& [s,nodes] : meshvelsymbcnodes)
223 : : m_meshvelsymbcnodes.insert( begin(nodes), end(nodes) );
224 : :
225 : : // Prepare unique sets of boundary nodes at which ALE moves the boundary
226 : : // based on user-defined functions.
227 : : std::unordered_map< int, std::unordered_set< std::size_t > > movenodes;
228 : : std::size_t i = 0;
229 [ + + ]: 1359 : for (const auto& m : g_inputdeck.get< tag::ale, tag::move >()) {
230 [ + + ]: 496 : for (const auto& s : m.get< tag::sideset >()) {
231 : 248 : auto k = bnode.find(static_cast<int>(s));
232 [ + + ]: 248 : if (k != end(bnode)) {
233 [ + - ]: 186 : auto& n = movenodes[ k->first ]; // associate set id
234 [ + + ]: 9145 : for (auto g : k->second) { // node ids on side set
235 : : n.insert( tk::cref_find(m_lid,g) ); // store local ids
236 : : }
237 : : }
238 : : }
239 : : // store all nodes from multiple side sets moved by this usrdef fn
240 : 248 : auto& n = std::get<2>(m_move[i]);
241 : : n.clear();
242 [ + + ]: 434 : for (const auto& [s,nodes] : movenodes) n.insert(begin(nodes), end(nodes));
243 : : // increment move ... end configuration block counter
244 : 248 : ++i;
245 : : }
246 : 1111 : }
247 : :
248 : : bool
249 : 67301 : ALE::move( std::size_t i ) const
250 : : // *****************************************************************************
251 : : // Find Dirichlet BCs on mesh velocity with prescribed movement
252 : : //! \param[in] i Local node id to check
253 : : //! \return True of node falls on a boundary that is prescribed to move
254 : : // *****************************************************************************
255 : : {
256 [ + + ]: 67301 : for (const auto& m : m_move)
257 [ - + ]: 2728 : if (std::get<2>(m).find(i) != end(std::get<2>(m)))
258 : : return true;
259 : :
260 : : return false;
261 : : }
262 : :
263 : : bool
264 : 801 : ALE::converged() const
265 : : // *****************************************************************************
266 : : // Query the solution of the Conjugrate Gradients linear solver
267 : : //! \return Solution to the Conjugate Gradients linear solve
268 : : // *****************************************************************************
269 : : {
270 [ - + ]: 1602 : return m_conjugategradients[ thisIndex ].ckLocal()->converged();
271 : : }
272 : :
273 : : void
274 : 1111 : ALE::start(
275 : : const tk::UnsMesh::Coords vel,
276 : : const std::vector< tk::real >& soundspeed,
277 : : CkCallback done,
278 : : const std::array< std::vector< tk::real >, 3 >& coord,
279 : : const tk::UnsMesh::Coords coordn,
280 : : const std::vector< tk::real >& vol0,
281 : : const std::vector< tk::real >& vol,
282 : : const std::unordered_map< int,
283 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm,
284 : : std::size_t initial,
285 : : std::size_t it,
286 : : tk::real t,
287 : : tk::real adt )
288 : : // *****************************************************************************
289 : : // Start computing new mesh velocity for ALE mesh motion
290 : : //! \param[in] vel Fluid velocity at mesh nodes
291 : : //! \param[in] soundspeed Speed of sound at mesh nodes
292 : : //! \param[in] done Function to continue with when mesh velocity has been
293 : : //! computed
294 : : //! \param[in] coord Mesh node coordinates
295 : : //! \param[in] coordn Mesh node coordinates at the previous time step
296 : : //! \param[in] vol0 Nodal mesh volumes at t=t0
297 : : //! \param[in] vol Nodal mesh volumes
298 : : //! \param[in] bnorm Face normals in boundary points associated to side sets
299 : : //! \param[in] initial Nonzero during the first time step stage, zero otherwise
300 : : //! \param[in] it Iteration count
301 : : //! \param[in] t Physics time
302 : : //! \param[in] adt alpha*dt of the RK time step
303 : : // *****************************************************************************
304 : : {
305 : 1111 : m_done = done;
306 : :
307 : : m_coord = coord;
308 : 1111 : m_soundspeed = soundspeed;
309 : 1111 : m_vol0 = vol0;
310 : 1111 : m_vol = vol;
311 : : m_bnorm = bnorm;
312 : 1111 : m_it = it;
313 : 1111 : m_t = t;
314 : 1111 : m_adt = adt;
315 : :
316 : : // assign mesh velocity
317 : 1111 : auto meshveltype = g_inputdeck.get< tag::ale, tag::mesh_velocity >();
318 [ + + ]: 1111 : if (meshveltype == ctr::MeshVelocityType::SINE) {
319 : :
320 : : // prescribe mesh velocity with a sine function during setup
321 [ + + ]: 279 : if (initial)
322 [ + + ]: 1812 : for (std::size_t i=0; i<m_w.nunk(); ++i)
323 : 1803 : m_w(i,0) = std::pow( std::sin(coord[0][i]*M_PI), 2.0 );
324 : :
325 [ + + ]: 832 : } else if (meshveltype == ctr::MeshVelocityType::FLUID) {
326 : :
327 : : // equate mesh velocity with fluid velocity
328 [ + + ]: 1726 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
329 [ + + ]: 1108793 : for (std::size_t i=0; i<vel[j].size(); ++i)
330 : 1107651 : m_w(i,j) = vel[j][i];
331 : :
332 [ + - ]: 248 : } else if (meshveltype == ctr::MeshVelocityType::USER_DEFINED) {
333 : :
334 : : // assign mesh velocity to sidesets from user-defined functions
335 [ + + ]: 496 : for (const auto& m : m_move)
336 [ + + ]: 248 : if (std::get<0>(m) == tk::ctr::UserTableType::VELOCITY) {
337 : 124 : auto meshvel = tk::sample<3>( t, std::get<1>(m) );
338 [ + + ]: 1922 : for (auto i : std::get<2>(m))
339 [ + + ]: 5394 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
340 : 3596 : m_w(i,j) = meshvel[j];
341 [ + - ]: 124 : } else if (std::get<0>(m) == tk::ctr::UserTableType::POSITION) {
342 : : auto eps = std::numeric_limits< tk::real >::epsilon();
343 [ + + ]: 124 : if (adt > eps) { // dt == 0 during setup
344 : 120 : auto pos = tk::sample<3>( t+adt, std::get<1>(m) );
345 [ + + ]: 1440 : for (auto i : std::get<2>(m))
346 [ + + ]: 3960 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
347 : 2640 : m_w(i,j) = (m_coord0[j][i] + pos[j] - coordn[j][i]) / adt;
348 : : }
349 : : }
350 : :
351 : : }
352 : :
353 : : // start computing the fluid vorticity
354 [ + + ]: 2222 : m_vorticity = tk::curl( coord, m_inpoel, vel );
355 : : // communicate vorticity sums to other chares on chare-boundary
356 [ + + ]: 1111 : if (m_nodeCommMap.empty()) {
357 : 154 : comvort_complete();
358 : : } else {
359 [ + + ]: 3927 : for (const auto& [c,n] : m_nodeCommMap) {
360 [ + - ]: 2970 : std::vector< std::array< tk::real, 3 > > v( n.size() );
361 : : std::size_t j = 0;
362 [ + + ]: 300090 : for (auto i : n) {
363 : 297120 : auto lid = tk::cref_find( m_lid, i );
364 : 297120 : v[j][0] = m_vorticity[0][lid];
365 : 297120 : v[j][1] = m_vorticity[1][lid];
366 : 297120 : v[j][2] = m_vorticity[2][lid];
367 : 297120 : ++j;
368 : : }
369 [ + - ][ + - ]: 5940 : thisProxy[c].comvort( std::vector<std::size_t>(begin(n),end(n)), v );
[ + - ][ - + ]
[ + - ][ - - ]
370 : : }
371 : : }
372 : 1111 : ownvort_complete();
373 : :
374 : : // start computing the fluid velocity divergence
375 [ + + ]: 2222 : m_veldiv = tk::div( coord, m_inpoel, vel );
376 : : // communicate vorticity sums to other chares on chare-boundary
377 [ + + ]: 1111 : if (m_nodeCommMap.empty()) {
378 : 154 : comdiv_complete();
379 : : } else {
380 [ + + ]: 3927 : for (const auto& [c,n] : m_nodeCommMap) {
381 [ + - ]: 2970 : std::vector< tk::real > v( n.size() );
382 : : std::size_t j = 0;
383 [ + + ]: 300090 : for (auto i : n) v[j++] = m_veldiv[ tk::cref_find( m_lid, i ) ];
384 [ + - ][ + - ]: 5940 : thisProxy[c].comdiv( std::vector<std::size_t>(begin(n),end(n)), v );
[ + - ][ - + ]
[ + - ][ - - ]
385 : : }
386 : : }
387 : 1111 : owndiv_complete();
388 : 1111 : }
389 : :
390 : : void
391 : 2970 : ALE::comvort( const std::vector< std::size_t >& gid,
392 : : const std::vector< std::array< tk::real, 3 > >& v )
393 : : // *****************************************************************************
394 : : // Receive contributions to vorticity on chare-boundaries
395 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
396 : : //! \param[in] v Partial contributions to chare-boundary nodes
397 : : // *****************************************************************************
398 : : {
399 : : Assert( v.size() == gid.size(), "Size mismatch" );
400 : : using tk::operator+=;
401 [ + + ]: 300090 : for (std::size_t i=0; i<gid.size(); ++i) m_vorticityc[ gid[i] ] += v[i];
402 : :
403 : : // When we have heard from all chares we communicate with, this chare is done
404 [ + + ]: 2970 : if (++m_nvort == m_nodeCommMap.size()) {
405 : 957 : m_nvort = 0;
406 : 957 : comvort_complete();
407 : : }
408 : 2970 : }
409 : :
410 : : void
411 : 2970 : ALE::comdiv( const std::vector< std::size_t >& gid,
412 : : const std::vector< tk::real >& v )
413 : : // *****************************************************************************
414 : : // Receive contributions to velocity divergence on chare-boundaries
415 : : //! \param[in] gid Global mesh node IDs at which we receive RHS contributions
416 : : //! \param[in] v Partial contributions to chare-boundary nodes
417 : : // *****************************************************************************
418 : : {
419 : : Assert( v.size() == gid.size(), "Size mismatch" );
420 : : using tk::operator+=;
421 [ + + ]: 300090 : for (std::size_t i=0; i<gid.size(); ++i) m_veldivc[ gid[i] ] += v[i];
422 : :
423 : : // When we have heard from all chares we communicate with, this chare is done
424 [ + + ]: 2970 : if (++m_ndiv == m_nodeCommMap.size()) {
425 : 957 : m_ndiv = 0;
426 : 957 : comdiv_complete();
427 : : }
428 : 2970 : }
429 : :
430 : : void
431 : 1111 : ALE::mergevel()
432 : : // *****************************************************************************
433 : : // Finalize computing fluid vorticity and velocity divergence for ALE
434 : : // *****************************************************************************
435 : : {
436 : : // combine own and communicated contributions to vorticity
437 [ + + ]: 259140 : for (const auto& [g,v] : m_vorticityc) {
438 : 258029 : auto lid = tk::cref_find( m_lid, g );
439 : 258029 : m_vorticity[0][lid] += v[0];
440 : 258029 : m_vorticity[1][lid] += v[1];
441 : 258029 : m_vorticity[2][lid] += v[2];
442 : : }
443 : : // clear fluid vorticity receive buffer
444 : 1111 : tk::destroy(m_vorticityc);
445 : :
446 : : // finish computing vorticity dividing the weak sum by the nodal volumes
447 [ + + ]: 4444 : for (std::size_t j=0; j<3; ++j)
448 [ + + ]: 3975612 : for (std::size_t p=0; p<m_vorticity[j].size(); ++p)
449 : 3972279 : m_vorticity[j][p] /= m_vol[p];
450 : :
451 : : // compute vorticity magnitude
452 [ + + ]: 1325204 : for (std::size_t p=0; p<m_vorticity[0].size(); ++p)
453 : 1324093 : m_vorticity[0][p] =
454 : 1324093 : tk::length( m_vorticity[0][p], m_vorticity[1][p], m_vorticity[2][p] );
455 : :
456 : : // get rid of the y and z vorticity components, since we just overwrote the
457 : : // x component with the magnitude
458 : : tk::destroy( m_vorticity[1] );
459 : : tk::destroy( m_vorticity[2] );
460 : :
461 : : // compute max vorticity magnitude
462 : : auto maxv =
463 [ + - ]: 2222 : *std::max_element( m_vorticity[0].cbegin(), m_vorticity[0].cend() );
464 : :
465 : : // combine own and communicated contributions to velocidy divergence
466 [ + + ]: 259140 : for (const auto& [g,v] : m_veldivc)
467 : 258029 : m_veldiv[ tk::cref_find( m_lid, g ) ] += v;
468 : : // clear velocity divergence receive buffer
469 : 1111 : tk::destroy(m_veldivc);
470 : :
471 : : // finish computing velocity divergence dividing weak sum by the nodal volumes
472 [ + + ]: 1325204 : for (std::size_t p=0; p<m_veldiv.size(); ++p) m_veldiv[p] /= m_vol[p];
473 : :
474 : : // Compute max vorticity magnitude across all chares
475 [ + - ]: 1111 : contribute( sizeof(tk::real), &maxv, CkReduction::max_double,
476 : 1111 : CkCallback(CkReductionTarget(ALE,meshvelbc), thisProxy) );
477 : 1111 : }
478 : :
479 : : void
480 : 1111 : ALE::meshvelbc( tk::real maxv )
481 : : // *****************************************************************************
482 : : // Apply mesh velocity smoother boundary conditions for ALE mesh motion
483 : : //! \param[in] maxv The largest vorticity magnitude across the whole problem
484 : : // *****************************************************************************
485 : : {
486 : : std::size_t ignorebc = false;
487 : :
488 : : // smooth mesh velocity if needed
489 : 1111 : auto smoother = g_inputdeck.get< tag::ale, tag::smoother >();
490 : :
491 [ + + ]: 1111 : if (smoother == ctr::MeshVelocitySmootherType::LAPLACE) {
492 : :
493 : : // scale mesh velocity with a function of the fluid vorticity
494 [ + + ]: 677 : if (maxv > 1.0e-8) {
495 : 652 : auto mult = g_inputdeck.get< tag::ale, tag::vortmult >();
496 [ + + ]: 1780 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
497 [ + + ]: 1593997 : for (std::size_t p=0; p<m_vorticity[0].size(); ++p)
498 [ + + ]: 3185645 : m_w(p,j) *= std::max( 0.0, 1.0 - mult*m_vorticity[0][p]/maxv );
499 : : }
500 : :
501 : : // Set mesh velocity smoother linear solve boundary conditions
502 : : std::unordered_map< std::size_t,
503 : : std::vector< std::pair< bool, tk::real > > > wbc;
504 : :
505 : : // Dirichlet BCs on mesh velocity with prescribed movement
506 [ + + ]: 925 : for (const auto& m : m_move)
507 [ + + ]: 248 : if (std::get<0>(m) == tk::ctr::UserTableType::VELOCITY) {
508 : 124 : auto meshvel = tk::sample<3>( m_t, std::get<1>(m) );
509 [ + + ]: 1922 : for (auto i : std::get<2>(m))
510 [ + + ]: 5394 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
511 : 3596 : m_w(i,j) = meshvel[j];
512 [ + - ]: 124 : } else if (std::get<0>(m) == tk::ctr::UserTableType::POSITION) {
513 : : auto eps = std::numeric_limits< tk::real >::epsilon();
514 [ + + ]: 124 : if (m_adt > eps) {
515 : 120 : ignorebc = m_it > 0;
516 [ + + ]: 1440 : for (auto i : std::get<2>(m))
517 [ + - ]: 1320 : if (m_meshveldirbcnodes.find(i) != end(m_meshveldirbcnodes))
518 [ + - ]: 2640 : wbc[i] = {{ {false,0}, {false,0}, {false,0} }};
519 : : }
520 : : }
521 : :
522 : : // Dirichlet BCs where user specified mesh velocity BCs
523 [ + + ]: 30375 : for (auto i : m_meshveldirbcnodes)
524 [ + - ][ + + ]: 58032 : if (not move(i)) wbc[i] = {{ {true,0}, {true,0}, {true,0} }};
[ + - ]
525 : :
526 : : // initialize mesh velocity smoother linear solver
527 : 677 : m_conjugategradients[ thisIndex ].ckLocal()->
528 [ + - ][ + - ]: 2708 : init( m_w.flat(), {}, wbc, ignorebc,
[ + - ][ - + ]
[ - - ]
529 [ + - ][ + - ]: 2031 : CkCallback(CkIndex_ALE::applied(nullptr), thisProxy[thisIndex]) );
[ - + ][ - + ]
[ - - ][ - - ]
530 : :
531 [ + + ]: 434 : } else if (smoother == ctr::MeshVelocitySmootherType::HELMHOLTZ) {
532 : :
533 : : // Set scalar potential linear solve boundary conditions
534 : : std::unordered_map< std::size_t,
535 : : std::vector< std::pair< bool, tk::real > > > pbc;
536 : :
537 : : // Dirichlet BCs where user specified mesh velocity BCs
538 [ + + ][ + - ]: 8029 : for (auto i : m_meshveldirbcnodes) pbc[i] = {{ {true,0} }};
539 : :
540 : : // prepare velocity divergence as weak sum required for Helmholtz solve
541 [ + - ]: 124 : decltype(m_veldiv) wveldiv = m_veldiv;
542 [ + + ]: 10013 : for (std::size_t p=0; p<wveldiv.size(); ++p) wveldiv[p] *= m_vol[p];
543 : :
544 : : // initialize Helmholtz decomposition linear solver
545 : 124 : m_conjugategradients[ thisIndex ].ckLocal()->
546 [ + - ][ + - ]: 496 : init( {}, wveldiv, pbc, ignorebc,
[ + - ][ - - ]
547 [ + - ][ + - ]: 372 : CkCallback(CkIndex_ALE::applied(nullptr), thisProxy[thisIndex]) );
[ - + ][ - + ]
[ - - ][ - - ]
548 : :
549 : : } else {
550 : :
551 : : // continue
552 : 310 : applied();
553 : :
554 : : }
555 : 1111 : }
556 : :
557 : : void
558 : 1111 : ALE::applied( [[maybe_unused]] CkDataMsg* msg )
559 : : // *****************************************************************************
560 : : // Solve mesh velocity linear solve for ALE mesh motion
561 : : // *****************************************************************************
562 : : {
563 : : //if (msg != nullptr) {
564 : : // auto *norm = static_cast< tk::real * >( msg->getData() );
565 : : // std::cout << "applied: " << *norm << '\n';
566 : : //}
567 : :
568 : 1111 : auto smoother = g_inputdeck.get< tag::ale, tag::smoother >();
569 : :
570 [ + + ]: 1111 : if (smoother == ctr::MeshVelocitySmootherType::LAPLACE) {
571 : :
572 [ + - ]: 2031 : m_conjugategradients[ thisIndex ].ckLocal()->solve(
573 : : g_inputdeck.get< tag::ale, tag::maxit >(),
574 : : g_inputdeck.get< tag::ale, tag::tolerance >(),
575 [ + - ][ + - ]: 2031 : CkCallback(CkIndex_ALE::solved(nullptr), thisProxy[thisIndex]) );
[ - + ][ - + ]
[ - - ][ - - ]
576 : :
577 [ + + ]: 434 : } else if (smoother == ctr::MeshVelocitySmootherType::HELMHOLTZ) {
578 : :
579 [ + - ]: 372 : m_conjugategradients[ thisIndex ].ckLocal()->solve(
580 : : g_inputdeck.get< tag::ale, tag::maxit >(),
581 : : g_inputdeck.get< tag::ale, tag::tolerance >(),
582 [ + - ][ + - ]: 372 : CkCallback(CkIndex_ALE::helmholtz(nullptr), thisProxy[thisIndex]) );
[ - + ][ - + ]
[ - - ][ - - ]
583 : :
584 : : } else {
585 : :
586 : 310 : solved();
587 : :
588 : : }
589 : 1111 : }
590 : :
591 : : void
592 : 124 : ALE::helmholtz( [[maybe_unused]] CkDataMsg* msg )
593 : : // *****************************************************************************
594 : : // Compute the gradient of the scalar potential for ALE
595 : : // *****************************************************************************
596 : : {
597 : : //if (msg != nullptr) {
598 : : // auto *norm = static_cast< tk::real * >( msg->getData() );
599 : : // std::cout << "solved: " << *norm << '\n';
600 : : //}
601 : :
602 : : // compute gradient of scalar potential for ALE (own portion)
603 [ + - ][ + - ]: 248 : m_gradpot = tk::grad( m_coord, m_inpoel,
604 [ - + ][ - + ]: 372 : m_conjugategradients[ thisIndex ].ckLocal()->solution() );
[ - - ]
605 : :
606 : : // communicate scalar potential sums to other chares on chare-boundary
607 [ - + ]: 124 : if (m_nodeCommMap.empty())
608 : 0 : compot_complete();
609 : : else
610 [ + + ]: 496 : for (const auto& [c,n] : m_nodeCommMap) {
611 [ + - ]: 372 : std::vector< std::array< tk::real, 3 > > v( n.size() );
612 : : std::size_t j = 0;
613 [ + + ]: 6386 : for (auto i : n) {
614 : 6014 : auto lid = tk::cref_find( m_lid, i );
615 : 6014 : v[j][0] = m_gradpot[0][lid];
616 : 6014 : v[j][1] = m_gradpot[1][lid];
617 : 6014 : v[j][2] = m_gradpot[2][lid];
618 : 6014 : ++j;
619 : : }
620 [ + - ][ + - ]: 744 : thisProxy[c].compot( std::vector<std::size_t>(begin(n),end(n)), v );
[ + - ][ - + ]
[ + - ][ - - ]
621 : : }
622 : 124 : ownpot_complete();
623 : 124 : }
624 : :
625 : : void
626 : 372 : ALE::compot( const std::vector< std::size_t >& gid,
627 : : const std::vector< std::array< tk::real, 3 > >& v )
628 : : // *****************************************************************************
629 : : // Receive contributions to scalar potential gradient on chare-boundaries
630 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
631 : : //! \param[in] v Partial contributions to chare-boundary nodes
632 : : // *****************************************************************************
633 : : {
634 : : Assert( v.size() == gid.size(), "Size mismatch" );
635 : : using tk::operator+=;
636 [ + + ]: 6386 : for (std::size_t i=0; i<gid.size(); ++i) m_gradpotc[ gid[i] ] += v[i];
637 : :
638 : : // When we have heard from all chares we communicate with, this chare is done
639 [ + + ]: 372 : if (++m_npot == m_nodeCommMap.size()) {
640 : 124 : m_npot = 0;
641 : 124 : compot_complete();
642 : : }
643 : 372 : }
644 : :
645 : : void
646 : 124 : ALE::gradpot()
647 : : // *****************************************************************************
648 : : // Finalize computing gradient of the scalar potential for ALE
649 : : // *****************************************************************************
650 : : {
651 : : // combine own and communicated contributions to scalar potential gradient
652 [ + + ]: 4712 : for (const auto& [g,v] : m_gradpotc) {
653 : 4588 : auto lid = tk::cref_find( m_lid, g );
654 : 4588 : m_gradpot[0][lid] += v[0];
655 : 4588 : m_gradpot[1][lid] += v[1];
656 : 4588 : m_gradpot[2][lid] += v[2];
657 : : }
658 : : // clear receive buffer
659 : 124 : tk::destroy(m_gradpotc);
660 : :
661 : : // finish computing the gradient dividing weak sum by the nodal volumes
662 [ + + ]: 496 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
663 [ + + ]: 30039 : for (std::size_t p=0; p<m_gradpot[j].size(); ++p)
664 : 29667 : m_gradpot[j][p] /= m_vol[p];
665 : :
666 : 124 : solved();
667 : 124 : }
668 : :
669 : : void
670 : 1111 : ALE::solved( [[maybe_unused]] CkDataMsg* msg )
671 : : // *****************************************************************************
672 : : // Mesh smoother linear solver converged
673 : : // *****************************************************************************
674 : : {
675 : : //if (msg != nullptr) {
676 : : // auto *normres = static_cast< tk::real * >( msg->getData() );
677 : : // std::cout << "solved: " << *normres << '\n';
678 : : //}
679 : :
680 : 1111 : auto smoother = g_inputdeck.get< tag::ale, tag::smoother >();
681 : :
682 [ + + ]: 1111 : if (smoother == ctr::MeshVelocitySmootherType::LAPLACE) {
683 : :
684 : : // Read out linear solution
685 : 1354 : auto w = m_conjugategradients[ thisIndex ].ckLocal()->solution();
686 : :
687 : : // Assign mesh velocity from linear solution skipping dimensions that are
688 : : // not allowed to move
689 [ + + ]: 1850 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
690 [ + + ]: 1653928 : for (std::size_t i=0; i<m_w.nunk(); ++i)
691 : 1652755 : m_w(i,j) = w[i*m_w.nprop()+j];
692 : :
693 [ + + ]: 434 : } else if (smoother == ctr::MeshVelocitySmootherType::HELMHOLTZ) {
694 : :
695 : 124 : auto a1 = g_inputdeck.get< tag::ale, tag::vortmult >();
696 [ + + ]: 496 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
697 [ + + ]: 30039 : for (std::size_t p=0; p<m_w.nunk(); ++p)
698 : 29667 : m_w(p,j) += a1 * (m_gradpot[j][p] - m_w(p,j));
699 : :
700 : : }
701 : :
702 : : // continue to applying a mesh force to the mesh velocity
703 : 1111 : startforce();
704 : 1111 : }
705 : :
706 : : void
707 : 1111 : ALE::startforce()
708 : : // *****************************************************************************
709 : : // Compute mesh force for the ALE mesh velocity
710 : : //! \details Compute mesh forces. See Sec.4 in Bakosi, Waltz, Morgan, Improved
711 : : //! ALE mesh velocities for complex flows, International Journal for Numerical
712 : : //! Methods in Fluids, 2017.
713 : : // *****************************************************************************
714 : : {
715 : : const auto& x = m_coord[0];
716 : : const auto& y = m_coord[1];
717 : : const auto& z = m_coord[2];
718 : :
719 : : // compute pseudo-pressure gradient for mesh force
720 : : const auto& f = g_inputdeck.get< tag::ale, tag::meshforce >();
721 : : m_wf.fill( 0.0 );
722 [ + + ]: 5107021 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
723 : : // access node IDs
724 : : const std::array< std::size_t, 4 >
725 : 5105910 : N{{m_inpoel[e*4+0], m_inpoel[e*4+1], m_inpoel[e*4+2], m_inpoel[e*4+3]}};
726 : : // compute element Jacobi determinant
727 : : const std::array< tk::real, 3 >
728 : 5105910 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
729 : 5105910 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
730 : 5105910 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
731 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
732 : 5105910 : auto J24 = J/24.0;
733 : : Assert( J > 0, "Element Jacobian non-positive" );
734 : :
735 : : // shape function derivatives, nnode*ndim [4][3]
736 : : std::array< std::array< tk::real, 3 >, 4 > grad;
737 : 5105910 : grad[1] = tk::crossdiv( ca, da, J );
738 : 5105910 : grad[2] = tk::crossdiv( da, ba, J );
739 : 5105910 : grad[3] = tk::crossdiv( ba, ca, J );
740 [ + + ]: 20423640 : for (std::size_t i=0; i<3; ++i)
741 : 15317730 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
742 : :
743 : : // max sound speed across nodes of element
744 : 5105910 : tk::real c = 1.0e-3; // floor on sound speed
745 [ + + ][ + + ]: 25529550 : for (std::size_t a=0; a<4; ++a) c = std::max( c, m_soundspeed[N[a]] );
746 : :
747 : : // mesh force in nodes
748 : 5105910 : auto V = J/6.0;
749 : 5105910 : auto L = std::cbrt(V); // element length scale, L=V^(1/3)
750 : 5105910 : auto dv = m_vol0[e] / V;
751 : 5105910 : std::array< tk::real, 4 > q{0,0,0,0};
752 [ + + ]: 25529550 : for (std::size_t a=0; a<4; ++a) {
753 : 20423640 : auto du = L * m_veldiv[N[a]];
754 : 20423640 : q[a] = - du*(f[0]*c + f[1]*std::abs(du) + f[2]*du*du)
755 : 20423640 : + f[3]*c*c*std::abs(dv-1.0);
756 : : }
757 : :
758 : : // scatter add pseudo-pressure gradient
759 [ + + ]: 25529550 : for (std::size_t a=0; a<4; ++a)
760 [ + + ]: 102118200 : for (std::size_t b=0; b<4; ++b)
761 [ + + ]: 326778240 : for (std::size_t i=0; i<3; ++i)
762 : 245083680 : m_wf(N[a],i) += J24 * grad[b][i] * q[b];
763 : : }
764 : :
765 : : // communicate mesh force sums to other chares on chare-boundary
766 [ + + ]: 1111 : if (m_nodeCommMap.empty()) {
767 : 154 : comfor_complete();
768 : : } else {
769 [ + + ]: 3927 : for (const auto& [c,n] : m_nodeCommMap) {
770 [ + - ]: 2970 : std::vector< std::array< tk::real, 3 > > w( n.size() );
771 : : std::size_t j = 0;
772 [ + + ]: 300090 : for (auto i : n) {
773 : 297120 : auto lid = tk::cref_find( m_lid, i );
774 : 297120 : w[j][0] = m_wf(lid,0);
775 : 297120 : w[j][1] = m_wf(lid,1);
776 : 297120 : w[j][2] = m_wf(lid,2);
777 : 297120 : ++j;
778 : : }
779 [ + - ][ + - ]: 5940 : thisProxy[c].comfor(std::vector<std::size_t>(begin(n),end(n)), w);
[ + - ][ - + ]
[ + - ][ - - ]
780 : : }
781 : : }
782 : 1111 : ownfor_complete();
783 : 1111 : }
784 : :
785 : : void
786 : 2970 : ALE::comfor( const std::vector< std::size_t >& gid,
787 : : const std::vector< std::array< tk::real, 3 > >& w )
788 : : // *****************************************************************************
789 : : // Receive contributions to ALE mesh force on chare-boundaries
790 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
791 : : //! \param[in] w Partial contributions to chare-boundary nodes
792 : : // *****************************************************************************
793 : : {
794 : : Assert( w.size() == gid.size(), "Size mismatch" );
795 : : using tk::operator+=;
796 [ + + ]: 300090 : for (std::size_t i=0; i<gid.size(); ++i) m_wfc[ gid[i] ] += w[i];
797 : :
798 : : // When we have heard from all chares we communicate with, this chare is done
799 [ + + ]: 2970 : if (++m_nwf == m_nodeCommMap.size()) {
800 : 957 : m_nwf = 0;
801 : 957 : comfor_complete();
802 : : }
803 : 2970 : }
804 : :
805 : : void
806 : 1111 : ALE::meshforce()
807 : : // *****************************************************************************
808 : : // Apply mesh force
809 : : // *****************************************************************************
810 : : {
811 : : // combine own and communicated contributions to mesh force
812 [ + + ]: 259140 : for (const auto& [g,w] : m_wfc) {
813 : 258029 : auto lid = tk::cref_find( m_lid, g );
814 : 258029 : m_wf(lid,0) += w[0];
815 : 258029 : m_wf(lid,1) += w[1];
816 : 258029 : m_wf(lid,2) += w[2];
817 : : }
818 : : // clear receive buffer
819 : 1111 : tk::destroy(m_wfc);
820 : :
821 : : // finish computing the mesh force by dviding weak sum by the nodal volumes
822 [ + + ]: 4444 : for (std::size_t j=0; j<3; ++j)
823 [ + + ]: 3975612 : for (std::size_t p=0; p<m_wf.nunk(); ++p)
824 : 3972279 : m_wf(p,j) /= m_vol[p];
825 : :
826 : : // advance mesh velocity in time due to pseudo-pressure gradient mesh force
827 [ + + ]: 3586 : for (auto j : g_inputdeck.get< tag::ale, tag::mesh_motion >())
828 [ + + ]: 1874803 : for (std::size_t i=0; i<m_w.nunk(); ++i)
829 : : // This is likely incorrect. It should be m_w = m_w0 + ...
830 : 1872328 : m_w(i,j) += m_adt * m_wf(i,j);
831 : :
832 : : // Enforce mesh velocity Dirichlet BCs where user specfied but did not
833 : : // prescribe a move
834 [ + + ]: 38714 : for (auto i : m_meshveldirbcnodes)
835 [ + + ]: 37603 : if (not move(i)) m_w(i,0) = m_w(i,1) = m_w(i,2) = 0.0;
836 : :
837 : : // On meshvel symmetry BCs remove normal component of mesh velocity
838 : : const auto& sbc = g_inputdeck.get< tag::ale, tag::symmetry >();
839 [ + + ]: 297843 : for (auto p : m_meshvelsymbcnodes) {
840 [ + + ]: 2077124 : for (const auto& s : sbc) {
841 : 1780392 : auto j = m_bnorm.find(static_cast<int>(s));
842 [ + + ]: 1780392 : if (j != end(m_bnorm)) {
843 : : auto i = j->second.find(p);
844 [ + + ]: 1380660 : if (i != end(j->second)) {
845 : : std::array< tk::real, 3 >
846 : 320550 : n{ i->second[0], i->second[1], i->second[2] },
847 : 320550 : v{ m_w(p,0), m_w(p,1), m_w(p,2) };
848 : : auto v_dot_n = tk::dot( v, n );
849 : : // symbc: remove normal component of mesh velocity
850 : 320550 : m_w(p,0) -= v_dot_n * n[0];
851 : 320550 : m_w(p,1) -= v_dot_n * n[1];
852 : 320550 : m_w(p,2) -= v_dot_n * n[2];
853 : : }
854 : : }
855 : : }
856 : : }
857 : :
858 : : // Activate SDAG wait for re-computing prerequisites for ALE
859 [ + - ]: 1111 : thisProxy[ thisIndex ].wait4vel();
860 [ + - ]: 1111 : thisProxy[ thisIndex ].wait4pot();
861 [ + - ]: 1111 : thisProxy[ thisIndex ].wait4for();
862 : :
863 : 1111 : m_done.send();
864 : 1111 : }
865 : :
866 : : #include "NoWarning/ale.def.h"
|
