Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/DG.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 DG advances a system of PDEs with the discontinuous Galerkin scheme
9 : : \details DG advances a system of partial differential equations (PDEs) using
10 : : discontinuous Galerkin (DG) finite element (FE) spatial discretization (on
11 : : tetrahedron elements) combined with Runge-Kutta (RK) time stepping.
12 : : \see The documentation in DG.h.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <algorithm>
17 : : #include <numeric>
18 : : #include <sstream>
19 : :
20 : : #include "DG.hpp"
21 : : #include "Discretization.hpp"
22 : : #include "DGPDE.hpp"
23 : : #include "DiagReducer.hpp"
24 : : #include "DerivedData.hpp"
25 : : #include "ElemDiagnostics.hpp"
26 : : #include "Inciter/InputDeck/InputDeck.hpp"
27 : : #include "Refiner.hpp"
28 : : #include "Limiter.hpp"
29 : : #include "Reorder.hpp"
30 : : #include "Vector.hpp"
31 : : #include "Around.hpp"
32 : : #include "Integrate/Basis.hpp"
33 : : #include "FieldOutput.hpp"
34 : : #include "ChareStateCollector.hpp"
35 : : #include "PDE/MultiMat/MultiMatIndexing.hpp"
36 : :
37 : : namespace inciter {
38 : :
39 : : extern ctr::InputDeck g_inputdeck;
40 : : extern std::vector< DGPDE > g_dgpde;
41 : :
42 : : //! Runge-Kutta coefficients
43 : : static const std::array< std::array< tk::real, 3 >, 2 >
44 : : rkcoef{{ {{ 0.0, 3.0/4.0, 1.0/3.0 }}, {{ 1.0, 1.0/4.0, 2.0/3.0 }} }};
45 : :
46 : : //! Implicit-Explicit Runge-Kutta Coefficients
47 : : static const tk::real rk_gamma = (2.0-std::sqrt(2.0))/2.0;
48 : : static const tk::real rk_delta = -2.0*std::sqrt(2.0)/3.0;
49 : : static const tk::real c2 =
50 : : (27.0 + std::pow(2187.0-1458.0*std::sqrt(2.0),1.0/3.0)
51 : : + 9.0*std::pow(3.0+2.0*std::sqrt(2.0),1.0/3.0))/54.0;
52 : : static const tk::real c3 = c2/(6.0*std::pow(c2,2.0)-3.0*c2+1.0);
53 : : static const tk::real b2 = (3.0*c2-1.0)/(6.0*std::pow(c2,2.0));
54 : : static const tk::real b3 =
55 : : (6.0*std::pow(c2,2.0)-3.0*c2+1.0)/(6.0*std::pow(c2,2.0));
56 : : static const tk::real a22_impl = c2;
57 : : static const tk::real a21_expl = c2;
58 : : static const tk::real a32_expl = c3;
59 : : static const tk::real a33_impl =
60 : : (1.0/6.0-b2*std::pow(c2,2.0)-b3*c2*c3)/(b3*(c3-c2));
61 : : static const tk::real a32_impl = a33_impl-c3;
62 : : static const std::array< std::array< tk::real, 3 >, 2 >
63 : : expl_rkcoef{{ {{ 0.0, 0.0, b2 }},
64 : : {{ a21_expl, a32_expl, b3 }} }};
65 : : static const std::array< std::array< tk::real, 3 >, 2>
66 : : impl_rkcoef{{ {{ 0.0, a32_impl, b2 }},
67 : : {{ a22_impl, a33_impl, b3}} }};
68 : :
69 : : } // inciter::
70 : :
71 : : extern tk::CProxy_ChareStateCollector stateProxy;
72 : :
73 : : using inciter::DG;
74 : :
75 : 671 : DG::DG( const CProxy_Discretization& disc,
76 : : const CProxy_Ghosts& ghostsproxy,
77 : : const std::map< int, std::vector< std::size_t > >& bface,
78 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
79 : 671 : const std::vector< std::size_t >& triinpoel ) :
80 : : m_disc( disc ),
81 : : m_ghosts( ghostsproxy ),
82 : : m_ndof_NodalExtrm( 3 ), // for the first order derivatives in 3 directions
83 : : m_nsol( 0 ),
84 : : m_ninitsol( 0 ),
85 : : m_nlim( 0 ),
86 : : m_nnod( 0 ),
87 : : m_nrefine( 0 ),
88 : : m_nsmooth( 0 ),
89 : : m_nreco( 0 ),
90 : : m_nnodalExtrema( 0 ),
91 [ + - ]: 671 : m_nstiffeq( g_dgpde[Disc()->MeshId()].nstiffeq() ),
92 [ + - ]: 671 : m_nnonstiffeq( g_dgpde[Disc()->MeshId()].nnonstiffeq() ),
93 [ + - ]: 671 : m_u( Disc()->Inpoel().size()/4,
94 : 671 : g_inputdeck.get< tag::rdof >()*
95 : 671 : g_inputdeck.get< tag::ncomp >() ),
96 : : m_un( m_u.nunk(), m_u.nprop() ),
97 : 671 : m_p( m_u.nunk(), g_inputdeck.get< tag::rdof >()*
98 [ + - ][ + - ]: 671 : g_dgpde[Disc()->MeshId()].nprim() ),
99 : : m_lhs( m_u.nunk(),
100 : 671 : g_inputdeck.get< tag::ndof >()*
101 : 671 : g_inputdeck.get< tag::ncomp >() ),
102 : : m_rhs( m_u.nunk(), m_lhs.nprop() ),
103 : : m_rhsprev( m_u.nunk(), m_lhs.nprop() ),
104 : 671 : m_stiffrhs( m_u.nunk(), g_inputdeck.get< tag::ndof >()*
105 [ + - ][ + - ]: 671 : g_dgpde[Disc()->MeshId()].nstiffeq() ),
106 : 671 : m_stiffrhsprev( m_u.nunk(), g_inputdeck.get< tag::ndof >()*
107 [ + - ][ + - ]: 671 : g_dgpde[Disc()->MeshId()].nstiffeq() ),
108 [ + - ]: 671 : m_stiffEqIdx( g_dgpde[Disc()->MeshId()].nstiffeq() ),
109 [ + - ]: 671 : m_nonStiffEqIdx( g_dgpde[Disc()->MeshId()].nnonstiffeq() ),
110 : : m_mtInv(
111 : : tk::invMassMatTaylorRefEl(g_inputdeck.get< tag::rdof >()) ),
112 : : m_uNodalExtrm(),
113 : : m_pNodalExtrm(),
114 : : m_uNodalExtrmc(),
115 : : m_pNodalExtrmc(),
116 [ + - ]: 671 : m_npoin( Disc()->Coord()[0].size() ),
117 : : m_diag(),
118 : : m_nstage( 3 ),
119 : : m_stage( 0 ),
120 : : m_ndof(),
121 : : m_interface(),
122 : : m_numEqDof(),
123 : : m_uc(),
124 : : m_pc(),
125 : : m_ndofc(),
126 : : m_interfacec(),
127 : : m_initial( 1 ),
128 : : m_uElemfields( m_u.nunk(),
129 : : g_inputdeck.get< tag::ncomp >() ),
130 : : m_pElemfields( m_u.nunk(),
131 : 671 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
132 : : m_uNodefields( m_npoin,
133 : : g_inputdeck.get< tag::ncomp >() ),
134 : : m_pNodefields( m_npoin,
135 : 671 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
136 : : m_uNodefieldsc(),
137 : : m_pNodefieldsc(),
138 : : m_outmesh(),
139 : : m_boxelems(),
140 [ + - ][ + - ]: 7381 : m_shockmarker(m_u.nunk(), 1)
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
141 : : // *****************************************************************************
142 : : // Constructor
143 : : //! \param[in] disc Discretization proxy
144 : : //! \param[in] bface Boundary-faces mapped to side set ids
145 : : //! \param[in] triinpoel Boundary-face connectivity
146 : : // *****************************************************************************
147 : : {
148 [ + + ]: 671 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
149 [ + + ]: 633 : g_inputdeck.get< tag::cmd, tag::quiescence >())
150 [ + - ][ + - ]: 892 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
151 : : "DG" );
152 : :
153 : : // assign number of dofs for each equation in all pde systems
154 [ + - ][ + - ]: 671 : g_dgpde[Disc()->MeshId()].numEquationDofs(m_numEqDof);
155 : :
156 : : // Allocate storage for the vector of nodal extrema
157 [ + - ][ + - ]: 671 : m_uNodalExtrm.resize( Disc()->Bid().size(),
158 : 0 : std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
159 [ + - ]: 671 : g_inputdeck.get< tag::ncomp >() ) );
160 [ + - ][ + - ]: 671 : m_pNodalExtrm.resize( Disc()->Bid().size(),
161 : 0 : std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
162 [ + - ]: 671 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ) );
163 : :
164 : : // Initialization for the buffer vector of nodal extrema
165 [ + - ]: 671 : resizeNodalExtremac();
166 : :
167 : 671 : usesAtSync = true; // enable migration at AtSync
168 : :
169 : : // Enable SDAG wait for initially building the solution vector and limiting
170 [ + - ]: 671 : if (m_initial) {
171 [ + - ][ + - ]: 671 : thisProxy[ thisIndex ].wait4sol();
172 [ + - ][ + - ]: 671 : thisProxy[ thisIndex ].wait4refine();
173 [ + - ][ + - ]: 671 : thisProxy[ thisIndex ].wait4smooth();
174 [ + - ][ + - ]: 671 : thisProxy[ thisIndex ].wait4lim();
175 [ + - ][ + - ]: 671 : thisProxy[ thisIndex ].wait4nod();
176 [ + - ][ + - ]: 671 : thisProxy[ thisIndex ].wait4reco();
177 [ + - ][ + - ]: 1342 : thisProxy[ thisIndex ].wait4nodalExtrema();
178 : : }
179 : :
180 [ + - ][ + - ]: 1342 : m_ghosts[thisIndex].insert(m_disc, bface, triinpoel, m_u.nunk(),
181 [ + - ][ + - ]: 2013 : CkCallback(CkIndex_DG::resizeSolVectors(), thisProxy[thisIndex]));
[ - + ][ - + ]
[ - - ][ - - ]
182 : :
183 : : // global-sync to call doneInserting on m_ghosts
184 [ + - ]: 671 : auto meshid = Disc()->MeshId();
185 [ + - ]: 671 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
186 [ + - ][ - - ]: 671 : CkCallback(CkReductionTarget(Transporter,doneInsertingGhosts),
187 [ + - ]: 671 : Disc()->Tr()) );
188 : 671 : }
189 : :
190 : : void
191 : 529 : DG::registerReducers()
192 : : // *****************************************************************************
193 : : // Configure Charm++ reduction types
194 : : //! \details Since this is a [initnode] routine, the runtime system executes the
195 : : //! routine exactly once on every logical node early on in the Charm++ init
196 : : //! sequence. Must be static as it is called without an object. See also:
197 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
198 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
199 : : // *****************************************************************************
200 : : {
201 : 529 : ElemDiagnostics::registerReducers();
202 : 529 : }
203 : :
204 : : void
205 : 12084 : DG::ResumeFromSync()
206 : : // *****************************************************************************
207 : : // Return from migration
208 : : //! \details This is called when load balancing (LB) completes. The presence of
209 : : //! this function does not affect whether or not we block on LB.
210 : : // *****************************************************************************
211 : : {
212 [ - + ][ - - ]: 12084 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
213 : :
214 [ + - ]: 12084 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
215 : 12084 : }
216 : :
217 : : void
218 : 671 : DG::resizeSolVectors()
219 : : // *****************************************************************************
220 : : // Resize solution vectors after extension due to Ghosts and continue with setup
221 : : // *****************************************************************************
222 : : {
223 : : // Resize solution vectors, lhs and rhs by the number of ghost tets
224 : 671 : m_u.resize( myGhosts()->m_nunk );
225 : 671 : m_un.resize( myGhosts()->m_nunk );
226 : 671 : m_p.resize( myGhosts()->m_nunk );
227 : 671 : m_lhs.resize( myGhosts()->m_nunk );
228 : 671 : m_rhs.resize( myGhosts()->m_nunk );
229 : 671 : m_rhsprev.resize( myGhosts()->m_nunk );
230 : 671 : m_stiffrhs.resize( myGhosts()->m_nunk );
231 : 671 : m_stiffrhsprev.resize( myGhosts()->m_nunk );
232 : :
233 : : // Size communication buffer for solution and number of degrees of freedom
234 [ + + ]: 2684 : for (auto& n : m_ndofc) n.resize( myGhosts()->m_bid.size() );
235 [ + + ]: 2684 : for (auto& u : m_uc) u.resize( myGhosts()->m_bid.size() );
236 [ + + ]: 2684 : for (auto& p : m_pc) p.resize( myGhosts()->m_bid.size() );
237 [ + + ]: 1342 : for (auto& i : m_interfacec) i.resize( myGhosts()->m_bid.size() );
238 : :
239 : : // Initialize number of degrees of freedom in mesh elements
240 : 671 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
241 [ + + ]: 671 : if( pref )
242 : : {
243 : 134 : const auto ndofmax = g_inputdeck.get< tag::pref, tag::ndofmax >();
244 : 134 : m_ndof.resize( myGhosts()->m_nunk, ndofmax );
245 : : }
246 : : else
247 : : {
248 : 537 : const auto ndof = g_inputdeck.get< tag::ndof >();
249 : 537 : m_ndof.resize( myGhosts()->m_nunk, ndof );
250 : : }
251 : 671 : m_interface.resize( myGhosts()->m_nunk, 0 );
252 : :
253 : : // Ensure that we also have all the geometry and connectivity data
254 : : // (including those of ghosts)
255 : : Assert( myGhosts()->m_geoElem.nunk() == m_u.nunk(),
256 : : "GeoElem unknowns size mismatch" );
257 : :
258 : : // Signal the runtime system that all workers have received their adjacency
259 : 671 : std::vector< std::size_t > meshdata{ myGhosts()->m_initial, Disc()->MeshId() };
260 : 671 : contribute( meshdata, CkReduction::sum_ulong,
261 [ + - ][ + - ]: 2013 : CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
[ + - ][ - - ]
262 : 671 : }
263 : :
264 : : void
265 : 671 : DG::setup()
266 : : // *****************************************************************************
267 : : // Set initial conditions, generate lhs, output mesh
268 : : // *****************************************************************************
269 : : {
270 [ + + ]: 671 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
271 [ + + ]: 633 : g_inputdeck.get< tag::cmd, tag::quiescence >())
272 [ + - ][ + - ]: 892 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
273 : : "setup" );
274 : :
275 : 671 : auto d = Disc();
276 : :
277 : : // Basic error checking on sizes of element geometry data and connectivity
278 : : Assert( myGhosts()->m_geoElem.nunk() == m_lhs.nunk(),
279 : : "Size mismatch in DG::setup()" );
280 : :
281 : : // Compute left-hand side of discrete PDEs
282 : 671 : lhs();
283 : :
284 : : // Determine elements inside user-defined IC box
285 : 671 : g_dgpde[d->MeshId()].IcBoxElems( myGhosts()->m_geoElem,
286 : 671 : myGhosts()->m_fd.Esuel().size()/4, m_boxelems );
287 : :
288 : : // Compute volume of user-defined box IC
289 [ + - ][ - + ]: 671 : d->boxvol( {}, {}, 0 ); // punt for now
290 : :
291 : : // Query time history field output labels from all PDEs integrated
292 : : const auto& hist_points = g_inputdeck.get< tag::history_output, tag::point >();
293 [ - + ]: 671 : if (!hist_points.empty()) {
294 : 0 : std::vector< std::string > histnames;
295 [ - - ]: 0 : auto n = g_dgpde[d->MeshId()].histNames();
296 [ - - ]: 0 : histnames.insert( end(histnames), begin(n), end(n) );
297 [ - - ]: 0 : d->histheader( std::move(histnames) );
298 : : }
299 : :
300 : : // If working with IMEX-RK, Store stiff equations into m_stiffEqIdx
301 [ - + ]: 671 : if (g_inputdeck.get< tag::imex_runge_kutta >())
302 : : {
303 : 0 : g_dgpde[Disc()->MeshId()].setStiffEqIdx(m_stiffEqIdx);
304 : 0 : g_dgpde[Disc()->MeshId()].setNonStiffEqIdx(m_nonStiffEqIdx);
305 : : }
306 : 671 : }
307 : :
308 : : void
309 : 671 : DG::box( tk::real v, const std::vector< tk::real >& )
310 : : // *****************************************************************************
311 : : // Receive total box IC volume and set conditions in box
312 : : //! \param[in] v Total volume within user-specified box
313 : : // *****************************************************************************
314 : : {
315 : 671 : auto d = Disc();
316 : :
317 : : // Store user-defined box IC volume
318 : 671 : d->Boxvol() = v;
319 : :
320 : : // Set initial conditions for all PDEs
321 : 1342 : g_dgpde[d->MeshId()].initialize( m_lhs, myGhosts()->m_inpoel,
322 : 671 : myGhosts()->m_coord, m_boxelems, d->ElemBlockId(), m_u, d->T(),
323 : 671 : myGhosts()->m_fd.Esuel().size()/4 );
324 : 671 : g_dgpde[d->MeshId()].updatePrimitives( m_u, m_lhs, myGhosts()->m_geoElem, m_p,
325 : 671 : myGhosts()->m_fd.Esuel().size()/4, m_ndof );
326 : :
327 : : m_un = m_u;
328 : :
329 : : // Output initial conditions to file (regardless of whether it was requested)
330 [ + - ][ + - ]: 2013 : startFieldOutput( CkCallback(CkIndex_DG::start(), thisProxy[thisIndex]) );
[ - + ][ - - ]
331 : 671 : }
332 : :
333 : : void
334 : 671 : DG::start()
335 : : // *****************************************************************************
336 : : // Start time stepping
337 : : // *****************************************************************************
338 : : {
339 : : // Free memory storing output mesh
340 : 671 : m_outmesh.destroy();
341 : :
342 : : // Start timer measuring time stepping wall clock time
343 : 671 : Disc()->Timer().zero();
344 : : // Zero grind-timer
345 : 671 : Disc()->grindZero();
346 : : // Start time stepping by computing the size of the next time step)
347 : 671 : next();
348 : 671 : }
349 : :
350 : : void
351 : 13426 : DG::startFieldOutput( CkCallback c )
352 : : // *****************************************************************************
353 : : // Start preparing fields for output to file
354 : : //! \param[in] c Function to continue with after the write
355 : : // *****************************************************************************
356 : : {
357 : : // No field output in benchmark mode or if field output frequency not hit
358 [ + + ][ + + ]: 13426 : if (g_inputdeck.get< tag::cmd, tag::benchmark >() || !fieldOutput()) {
359 : :
360 : 12101 : c.send();
361 : :
362 : : } else {
363 : :
364 : : // Optionally refine mesh for field output
365 : 1325 : auto d = Disc();
366 : :
367 [ + + ]: 1325 : if (refinedOutput()) {
368 : :
369 : 33 : const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
370 [ + - ][ + - ]: 66 : d->Ref()->outref( myGhosts()->m_fd.Bface(), {}, tr, c );
[ + - ][ + - ]
[ - - ]
371 : :
372 : : } else {
373 : :
374 : : // cut off ghosts from mesh connectivity and coordinates
375 : 1292 : const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
376 [ + - ][ + - ]: 2584 : extractFieldOutput( {}, d->Chunk(), d->Coord(), {}, {},
[ + + ][ - - ]
377 [ + - ]: 1292 : d->NodeCommMap(), myGhosts()->m_fd.Bface(), {}, tr, c );
378 : :
379 : : }
380 : :
381 : : }
382 : 13426 : }
383 : :
384 : : void
385 : 38265 : DG::next()
386 : : // *****************************************************************************
387 : : // Advance equations to next time step
388 : : // *****************************************************************************
389 : : {
390 : 38265 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
391 : :
392 : 38265 : auto d = Disc();
393 : :
394 [ + + ][ + + ]: 38265 : if (pref && m_stage == 0 && d->T() > 0)
[ + + ]
395 : 3272 : g_dgpde[d->MeshId()].eval_ndof( myGhosts()->m_nunk, myGhosts()->m_coord,
396 : 1636 : myGhosts()->m_inpoel,
397 : 1636 : myGhosts()->m_fd, m_u, m_p,
398 : : g_inputdeck.get< tag::pref, tag::indicator >(),
399 : : g_inputdeck.get< tag::ndof >(),
400 : : g_inputdeck.get< tag::pref, tag::ndofmax >(),
401 : : g_inputdeck.get< tag::pref, tag::tolref >(),
402 : 1636 : m_ndof );
403 : :
404 : : // communicate solution ghost data (if any)
405 [ + + ]: 38265 : if (myGhosts()->m_sendGhost.empty())
406 : 3390 : comsol_complete();
407 : : else
408 [ + + ]: 454800 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
409 [ + - ]: 385050 : std::vector< std::size_t > tetid( ghostdata.size() );
410 [ + - ][ + - ]: 770100 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
411 [ + - ][ + - ]: 770100 : prim( ghostdata.size() );
412 [ + - ][ + - ]: 385050 : std::vector< std::size_t > interface( ghostdata.size() );
413 [ + - ][ - - ]: 385050 : std::vector< std::size_t > ndof( ghostdata.size() );
414 : : std::size_t j = 0;
415 [ + + ]: 6854730 : for(const auto& i : ghostdata) {
416 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
417 : : "Sending solution ghost data" );
418 [ + - ]: 6469680 : tetid[j] = i;
419 [ + - ][ - + ]: 6469680 : u[j] = m_u[i];
420 [ + - ][ - + ]: 6469680 : prim[j] = m_p[i];
421 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0) {
422 : 395110 : ndof[j] = m_ndof[i];
423 : 395110 : interface[j] = m_interface[i];
424 : : }
425 : 6469680 : ++j;
426 : : }
427 [ + - ][ + - ]: 770100 : thisProxy[ cid ].comsol( thisIndex, m_stage, tetid, u, prim, interface, ndof );
[ + - ][ - - ]
428 : : }
429 : :
430 : 38265 : ownsol_complete();
431 : 38265 : }
432 : :
433 : : void
434 : 385050 : DG::comsol( int fromch,
435 : : std::size_t fromstage,
436 : : const std::vector< std::size_t >& tetid,
437 : : const std::vector< std::vector< tk::real > >& u,
438 : : const std::vector< std::vector< tk::real > >& prim,
439 : : const std::vector< std::size_t >& interface,
440 : : const std::vector< std::size_t >& ndof )
441 : : // *****************************************************************************
442 : : // Receive chare-boundary solution ghost data from neighboring chares
443 : : //! \param[in] fromch Sender chare id
444 : : //! \param[in] fromstage Sender chare time step stage
445 : : //! \param[in] tetid Ghost tet ids we receive solution data for
446 : : //! \param[in] u Solution ghost data
447 : : //! \param[in] prim Primitive variables in ghost cells
448 : : //! \param[in] interface Interface marker in ghost cells
449 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
450 : : //! \details This function receives contributions to the unlimited solution
451 : : //! from fellow chares.
452 : : // *****************************************************************************
453 : : {
454 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comsol()" );
455 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comsol()" );
456 : :
457 : 385050 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
458 : :
459 : 385050 : if (pref && fromstage == 0) {
460 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsol()" );
461 : : Assert( interface.size() == tetid.size(), "Size mismatch in DG::comsol()" );
462 : : }
463 : :
464 : : // Find local-to-ghost tet id map for sender chare
465 : 385050 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
466 : :
467 [ + + ]: 6854730 : for (std::size_t i=0; i<tetid.size(); ++i) {
468 : 6469680 : auto j = tk::cref_find( n, tetid[i] );
469 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
470 : : "Receiving solution non-ghost data" );
471 : 6469680 : auto b = tk::cref_find( myGhosts()->m_bid, j );
472 : : Assert( b < m_uc[0].size(), "Indexing out of bounds" );
473 : : Assert( b < m_pc[0].size(), "Indexing out of bounds" );
474 : 6469680 : m_uc[0][b] = u[i];
475 : 6469680 : m_pc[0][b] = prim[i];
476 [ + + ]: 6469680 : if (pref && fromstage == 0) {
477 : : Assert( b < m_ndofc[0].size(), "Indexing out of bounds" );
478 : 395110 : m_ndofc[0][b] = ndof[i];
479 : : Assert( b < m_interfacec[0].size(), "Indexing out of bounds" );
480 : 395110 : m_interfacec[0][b] = interface[i];
481 : : }
482 : : }
483 : :
484 : : // if we have received all solution ghost contributions from neighboring
485 : : // chares (chares we communicate along chare-boundary faces with), and
486 : : // contributed our solution to these neighbors, proceed to reconstructions
487 [ + + ]: 385050 : if (++m_nsol == myGhosts()->m_sendGhost.size()) {
488 : 34875 : m_nsol = 0;
489 : 34875 : comsol_complete();
490 : : }
491 : 385050 : }
492 : :
493 : : void
494 : 1325 : DG::extractFieldOutput(
495 : : const std::vector< std::size_t >& /*ginpoel*/,
496 : : const tk::UnsMesh::Chunk& chunk,
497 : : const tk::UnsMesh::Coords& coord,
498 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
499 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
500 : : const tk::NodeCommMap& nodeCommMap,
501 : : const std::map< int, std::vector< std::size_t > >& bface,
502 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
503 : : const std::vector< std::size_t >& triinpoel,
504 : : CkCallback c )
505 : : // *****************************************************************************
506 : : // Extract field output going to file
507 : : //! \param[in] chunk Field-output mesh chunk (connectivity and global<->local
508 : : //! id maps)
509 : : //! \param[in] coord Field-output mesh node coordinates
510 : : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
511 : : //! \param[in] nodeCommMap Field-output mesh node communication map
512 : : //! \param[in] bface Field-output meshndary-faces mapped to side set ids
513 : : //! \param[in] triinpoel Field-output mesh boundary-face connectivity
514 : : //! \param[in] c Function to continue with after the write
515 : : // *****************************************************************************
516 : : {
517 : : m_outmesh.chunk = chunk;
518 : : m_outmesh.coord = coord;
519 : 1325 : m_outmesh.triinpoel = triinpoel;
520 : : m_outmesh.bface = bface;
521 : : m_outmesh.nodeCommMap = nodeCommMap;
522 : :
523 : : const auto& inpoel = std::get< 0 >( chunk );
524 : :
525 : : // Evaluate element solution on incoming mesh
526 : 1325 : evalSolution( *Disc(), inpoel, coord, addedTets, m_ndof, m_u, m_p,
527 : 1325 : m_uElemfields, m_pElemfields, m_uNodefields, m_pNodefields );
528 : :
529 : : // Send node fields contributions to neighbor chares
530 [ + + ]: 1325 : if (nodeCommMap.empty())
531 : 161 : comnodeout_complete();
532 : : else {
533 : : const auto& lid = std::get< 2 >( chunk );
534 : 2328 : auto esup = tk::genEsup( inpoel, 4 );
535 [ + + ]: 11544 : for(const auto& [ch,nodes] : nodeCommMap) {
536 : : // Pack node field data in chare boundary nodes
537 : : std::vector< std::vector< tk::real > >
538 [ + - ][ + - ]: 31140 : lu( m_uNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
[ + - ]
539 : : std::vector< std::vector< tk::real > >
540 [ + - ][ + - ]: 31140 : lp( m_pNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
541 [ + + ]: 69172 : for (std::size_t f=0; f<m_uNodefields.nprop(); ++f) {
542 : : std::size_t j = 0;
543 [ + + ]: 603358 : for (auto g : nodes)
544 : 544566 : lu[f][j++] = m_uNodefields(tk::cref_find(lid,g),f);
545 : : }
546 [ + + ]: 21008 : for (std::size_t f=0; f<m_pNodefields.nprop(); ++f) {
547 : : std::size_t j = 0;
548 [ + + ]: 120108 : for (auto g : nodes)
549 : 109480 : lp[f][j++] = m_pNodefields(tk::cref_find(lid,g),f);
550 : : }
551 : : // Pack (partial) number of elements surrounding chare boundary nodes
552 [ + - ]: 10380 : std::vector< std::size_t > nesup( nodes.size() );
553 : : std::size_t j = 0;
554 [ + + ]: 107480 : for (auto g : nodes) {
555 : 97100 : auto i = tk::cref_find( lid, g );
556 : 97100 : nesup[j++] = esup.second[i+1] - esup.second[i];
557 : : }
558 [ + - ][ + - ]: 31140 : thisProxy[ch].comnodeout(
[ + - ][ - - ]
559 [ + - ][ - + ]: 20760 : std::vector<std::size_t>(begin(nodes),end(nodes)), nesup, lu, lp );
560 : : }
561 : : }
562 : :
563 [ + - ]: 1325 : ownnod_complete( c, addedTets );
564 : 1325 : }
565 : :
566 : : void
567 : 671 : DG::lhs()
568 : : // *****************************************************************************
569 : : // Compute left-hand side of discrete transport equations
570 : : // *****************************************************************************
571 : : {
572 : 671 : g_dgpde[Disc()->MeshId()].lhs( myGhosts()->m_geoElem, m_lhs );
573 : :
574 [ - + ]: 671 : if (!m_initial) stage();
575 : 671 : }
576 : :
577 : 38265 : void DG::refine()
578 : : // *****************************************************************************
579 : : // Add the protective layer for ndof refinement
580 : : // *****************************************************************************
581 : : {
582 : 38265 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
583 : :
584 : : // Combine own and communicated contributions of unreconstructed solution and
585 : : // degrees of freedom in cells (if p-adaptive)
586 [ + + ]: 6546210 : for (const auto& b : myGhosts()->m_bid) {
587 : : Assert( m_uc[0][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
588 : : Assert( m_pc[0][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
589 [ + + ]: 173340045 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
590 : 166870365 : m_u(b.first,c) = m_uc[0][b.second][c];
591 : : }
592 [ + + ]: 34324605 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
593 : 27854925 : m_p(b.first,c) = m_pc[0][b.second][c];
594 : : }
595 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0) {
596 : 395110 : m_ndof[ b.first ] = m_ndofc[0][ b.second ];
597 : 395110 : m_interface[ b.first ] = m_interfacec[0][ b.second ];
598 : : }
599 : : }
600 : :
601 [ + + ][ + + ]: 38265 : if (pref && m_stage==0) refine_ndof();
602 : :
603 [ + + ]: 38265 : if (myGhosts()->m_sendGhost.empty())
604 : 3390 : comrefine_complete();
605 : : else
606 [ + + ]: 454800 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
607 [ + - ]: 385050 : std::vector< std::size_t > tetid( ghostdata.size() );
608 [ + - ][ + - ]: 770100 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
609 [ + - ][ + - ]: 770100 : prim( ghostdata.size() );
610 [ + - ]: 385050 : std::vector< std::size_t > ndof( ghostdata.size() );
611 : : std::size_t j = 0;
612 [ + + ]: 6854730 : for(const auto& i : ghostdata) {
613 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending refined ndof "
614 : : "data" );
615 [ + + ]: 6469680 : tetid[j] = i;
616 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0) ndof[j] = m_ndof[i];
617 : 6469680 : ++j;
618 : : }
619 [ + - ][ + - ]: 770100 : thisProxy[ cid ].comrefine( thisIndex, tetid, ndof );
[ + - ][ - - ]
620 : : }
621 : :
622 : 38265 : ownrefine_complete();
623 : 38265 : }
624 : :
625 : : void
626 : 385050 : DG::comrefine( int fromch,
627 : : const std::vector< std::size_t >& tetid,
628 : : const std::vector< std::size_t >& ndof )
629 : : // *****************************************************************************
630 : : // Receive chare-boundary ghost data from neighboring chares
631 : : //! \param[in] fromch Sender chare id
632 : : //! \param[in] tetid Ghost tet ids we receive solution data for
633 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
634 : : //! \details This function receives contributions to the refined ndof data
635 : : //! from fellow chares.
636 : : // *****************************************************************************
637 : : {
638 : 385050 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
639 : :
640 : : if (pref && m_stage == 0)
641 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comrefine()" );
642 : :
643 : : // Find local-to-ghost tet id map for sender chare
644 : 385050 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
645 : :
646 [ + + ]: 6854730 : for (std::size_t i=0; i<tetid.size(); ++i) {
647 : 6469680 : auto j = tk::cref_find( n, tetid[i] );
648 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
649 : : "Receiving solution non-ghost data" );
650 : 6469680 : auto b = tk::cref_find( myGhosts()->m_bid, j );
651 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0) {
652 : : Assert( b < m_ndofc[1].size(), "Indexing out of bounds" );
653 : 395110 : m_ndofc[1][b] = ndof[i];
654 : : }
655 : : }
656 : :
657 : : // if we have received all solution ghost contributions from neighboring
658 : : // chares (chares we communicate along chare-boundary faces with), and
659 : : // contributed our solution to these neighbors, proceed to limiting
660 [ + + ]: 385050 : if (++m_nrefine == myGhosts()->m_sendGhost.size()) {
661 : 34875 : m_nrefine = 0;
662 : 34875 : comrefine_complete();
663 : : }
664 : 385050 : }
665 : :
666 : : void
667 : 38265 : DG::smooth()
668 : : // *****************************************************************************
669 : : // Smooth the refined ndof distribution
670 : : // *****************************************************************************
671 : : {
672 : 38265 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
673 : :
674 [ + + ]: 6546210 : for (const auto& b : myGhosts()->m_bid) {
675 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0)
676 : 395110 : m_ndof[ b.first ] = m_ndofc[1][ b.second ];
677 : : }
678 : :
679 [ + + ][ + + ]: 38265 : if (pref && m_stage==0) smooth_ndof();
680 : :
681 [ + + ]: 38265 : if (myGhosts()->m_sendGhost.empty())
682 : 3390 : comsmooth_complete();
683 : : else
684 [ + + ]: 454800 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
685 : 385050 : std::vector< std::size_t > tetid( ghostdata.size() );
686 : : std::vector< std::size_t > ndof;
687 : : std::size_t j = 0;
688 [ + + ]: 6854730 : for(const auto& i : ghostdata) {
689 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending ndof data" );
690 [ + + ]: 6469680 : tetid[j] = i;
691 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
[ + - ]
692 : 6469680 : ++j;
693 : : }
694 [ + - ][ + - ]: 770100 : thisProxy[ cid ].comsmooth( thisIndex, tetid, ndof );
[ + + ][ - - ]
695 : : }
696 : :
697 : 38265 : ownsmooth_complete();
698 : 38265 : }
699 : :
700 : : void
701 : 385050 : DG::comsmooth( int fromch,
702 : : const std::vector< std::size_t >& tetid,
703 : : const std::vector< std::size_t >& ndof )
704 : : // *****************************************************************************
705 : : // Receive chare-boundary ghost data from neighboring chares
706 : : //! \param[in] fromch Sender chare id
707 : : //! \param[in] tetid Ghost tet ids we receive solution data for
708 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
709 : : //! \details This function receives contributions to the smoothed ndof data
710 : : //! from fellow chares.
711 : : // *****************************************************************************
712 : : {
713 : 385050 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
714 : :
715 : : if (pref && m_stage == 0)
716 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsmooth()" );
717 : :
718 : 385050 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
719 : :
720 [ + + ]: 6854730 : for (std::size_t i=0; i<tetid.size(); ++i) {
721 : 6469680 : auto j = tk::cref_find( n, tetid[i] );
722 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4, "Receiving ndof data" );
723 : 6469680 : auto b = tk::cref_find( myGhosts()->m_bid, j );
724 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0) {
725 : : Assert( b < m_ndofc[2].size(), "Indexing out of bounds" );
726 : 395110 : m_ndofc[2][b] = ndof[i];
727 : : }
728 : : }
729 : :
730 [ + + ]: 385050 : if (++m_nsmooth == myGhosts()->m_sendGhost.size()) {
731 : 34875 : m_nsmooth = 0;
732 : 34875 : comsmooth_complete();
733 : : }
734 : 385050 : }
735 : :
736 : : void
737 : 38265 : DG::reco()
738 : : // *****************************************************************************
739 : : // Compute reconstructions
740 : : // *****************************************************************************
741 : : {
742 : 38265 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
743 : 38265 : const auto rdof = g_inputdeck.get< tag::rdof >();
744 : :
745 : : // Combine own and communicated contributions of unreconstructed solution and
746 : : // degrees of freedom in cells (if p-adaptive)
747 [ + + ]: 6546210 : for (const auto& b : myGhosts()->m_bid) {
748 [ + + ][ + + ]: 6469680 : if (pref && m_stage == 0) {
749 : 395110 : m_ndof[ b.first ] = m_ndofc[2][ b.second ];
750 : : }
751 : : }
752 : :
753 : 38265 : auto d = Disc();
754 [ + + ][ + + ]: 38265 : if (pref && m_stage==0) {
755 : 1770 : g_dgpde[d->MeshId()].resetAdapSol( myGhosts()->m_fd, m_u, m_p, m_ndof );
756 : : }
757 : :
758 [ + + ]: 38265 : if (rdof > 1)
759 : : // Reconstruct second-order solution and primitive quantities
760 : 49680 : g_dgpde[d->MeshId()].reconstruct( d->T(), myGhosts()->m_geoFace,
761 : 24840 : myGhosts()->m_geoElem,
762 : 24840 : myGhosts()->m_fd, myGhosts()->m_esup, myGhosts()->m_inpoel,
763 : 24840 : myGhosts()->m_coord, m_u, m_p, pref, m_ndof );
764 : :
765 : : // Send reconstructed solution to neighboring chares
766 [ + + ]: 38265 : if (myGhosts()->m_sendGhost.empty())
767 : 3390 : comreco_complete();
768 : : else
769 [ + + ]: 454800 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
770 [ + - ]: 385050 : std::vector< std::size_t > tetid( ghostdata.size() );
771 [ + - ][ + - ]: 770100 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
772 [ + - ]: 385050 : prim( ghostdata.size() );
773 : : std::size_t j = 0;
774 [ + + ]: 6854730 : for(const auto& i : ghostdata) {
775 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending reconstructed ghost "
776 : : "data" );
777 [ + - ]: 6469680 : tetid[j] = i;
778 [ + - ][ - + ]: 6469680 : u[j] = m_u[i];
779 [ + - ][ - + ]: 6469680 : prim[j] = m_p[i];
780 : 6469680 : ++j;
781 : : }
782 [ + - ][ + - ]: 770100 : thisProxy[ cid ].comreco( thisIndex, tetid, u, prim );
783 : : }
784 : :
785 : 38265 : ownreco_complete();
786 : 38265 : }
787 : :
788 : : void
789 : 385050 : DG::comreco( int fromch,
790 : : const std::vector< std::size_t >& tetid,
791 : : const std::vector< std::vector< tk::real > >& u,
792 : : const std::vector< std::vector< tk::real > >& prim )
793 : : // *****************************************************************************
794 : : // Receive chare-boundary reconstructed ghost data from neighboring chares
795 : : //! \param[in] fromch Sender chare id
796 : : //! \param[in] tetid Ghost tet ids we receive solution data for
797 : : //! \param[in] u Reconstructed high-order solution
798 : : //! \param[in] prim Limited high-order primitive quantities
799 : : //! \details This function receives contributions to the reconstructed solution
800 : : //! from fellow chares.
801 : : // *****************************************************************************
802 : : {
803 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comreco()" );
804 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comreco()" );
805 : :
806 : : // Find local-to-ghost tet id map for sender chare
807 : 385050 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
808 : :
809 [ + + ]: 6854730 : for (std::size_t i=0; i<tetid.size(); ++i) {
810 : 6469680 : auto j = tk::cref_find( n, tetid[i] );
811 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
812 : : "Receiving solution non-ghost data" );
813 : 6469680 : auto b = tk::cref_find( myGhosts()->m_bid, j );
814 : : Assert( b < m_uc[1].size(), "Indexing out of bounds" );
815 : : Assert( b < m_pc[1].size(), "Indexing out of bounds" );
816 : 6469680 : m_uc[1][b] = u[i];
817 : 6469680 : m_pc[1][b] = prim[i];
818 : : }
819 : :
820 : : // if we have received all solution ghost contributions from neighboring
821 : : // chares (chares we communicate along chare-boundary faces with), and
822 : : // contributed our solution to these neighbors, proceed to limiting
823 [ + + ]: 385050 : if (++m_nreco == myGhosts()->m_sendGhost.size()) {
824 : 34875 : m_nreco = 0;
825 : 34875 : comreco_complete();
826 : : }
827 : 385050 : }
828 : :
829 : : void
830 : 38265 : DG::nodalExtrema()
831 : : // *****************************************************************************
832 : : // Compute nodal extrema at chare-boundary nodes. Extrema at internal nodes
833 : : // are calculated in limiter function.
834 : : // *****************************************************************************
835 : : {
836 : 38265 : auto d = Disc();
837 : 38265 : auto gid = d->Gid();
838 : : auto bid = d->Bid();
839 : 38265 : const auto rdof = g_inputdeck.get< tag::rdof >();
840 : 38265 : const auto ncomp = m_u.nprop() / rdof;
841 : 38265 : const auto nprim = m_p.nprop() / rdof;
842 : :
843 : : // Combine own and communicated contributions of unlimited solution, and
844 : : // if a p-adaptive algorithm is used, degrees of freedom in cells
845 [ + - ][ + + ]: 6546210 : for (const auto& [boundary, localtet] : myGhosts()->m_bid) {
846 : : Assert( m_uc[1][localtet].size() == m_u.nprop(), "ncomp size mismatch" );
847 : : Assert( m_pc[1][localtet].size() == m_p.nprop(), "ncomp size mismatch" );
848 [ + + ]: 173340045 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
849 : 166870365 : m_u(boundary,c) = m_uc[1][localtet][c];
850 : : }
851 [ + + ]: 34324605 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
852 : 27854925 : m_p(boundary,c) = m_pc[1][localtet][c];
853 : : }
854 : : }
855 : :
856 : : // Initialize nodal extrema vector
857 : : auto large = std::numeric_limits< tk::real >::max();
858 [ + + ]: 1317765 : for(std::size_t i = 0; i<bid.size(); i++)
859 : : {
860 [ + + ]: 9304560 : for (std::size_t c=0; c<ncomp; ++c)
861 : : {
862 [ + + ]: 32100240 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
863 : : {
864 : 24075180 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
865 : 24075180 : auto min_mark = max_mark + 1;
866 : 24075180 : m_uNodalExtrm[i][max_mark] = -large;
867 : 24075180 : m_uNodalExtrm[i][min_mark] = large;
868 : : }
869 : : }
870 [ + + ]: 3460350 : for (std::size_t c=0; c<nprim; ++c)
871 : : {
872 [ + + ]: 8723400 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
873 : : {
874 : 6542550 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
875 : 6542550 : auto min_mark = max_mark + 1;
876 : 6542550 : m_pNodalExtrm[i][max_mark] = -large;
877 : 6542550 : m_pNodalExtrm[i][min_mark] = large;
878 : : }
879 : : }
880 : : }
881 : :
882 : : // Evaluate the max/min value for the chare-boundary nodes
883 [ + + ]: 38265 : if(rdof > 4) {
884 [ + - ]: 14220 : evalNodalExtrmRefEl(ncomp, nprim, m_ndof_NodalExtrm, d->bndel(),
885 [ + - ][ + - ]: 7110 : myGhosts()->m_inpoel, gid, bid, m_u, m_p, m_uNodalExtrm, m_pNodalExtrm);
886 : : }
887 : :
888 : : // Communicate extrema at nodes to other chares on chare-boundary
889 [ + + ]: 38265 : if (d->NodeCommMap().empty()) // in serial we are done
890 [ + - ]: 3390 : comnodalExtrema_complete();
891 : : else // send nodal extrema to chare-boundary nodes to fellow chares
892 : : {
893 [ + + ]: 419925 : for (const auto& [c,n] : d->NodeCommMap()) {
894 [ + - ][ + - ]: 770100 : std::vector< std::vector< tk::real > > g1( n.size() ), g2( n.size() );
895 : : std::size_t j = 0;
896 [ + + ]: 3269040 : for (auto i : n)
897 : : {
898 : 2883990 : auto p = tk::cref_find(d->Bid(),i);
899 [ + - ]: 2883990 : g1[ j ] = m_uNodalExtrm[ p ];
900 [ + - ]: 2883990 : g2[ j++ ] = m_pNodalExtrm[ p ];
901 : : }
902 [ + - ][ + - ]: 770100 : thisProxy[c].comnodalExtrema( std::vector<std::size_t>(begin(n),end(n)),
[ + - ][ - + ]
903 : : g1, g2 );
904 : : }
905 : : }
906 [ + - ]: 38265 : ownnodalExtrema_complete();
907 : 38265 : }
908 : :
909 : : void
910 : 385050 : DG::comnodalExtrema( const std::vector< std::size_t >& gid,
911 : : const std::vector< std::vector< tk::real > >& G1,
912 : : const std::vector< std::vector< tk::real > >& G2 )
913 : : // *****************************************************************************
914 : : // Receive contributions to nodal extrema on chare-boundaries
915 : : //! \param[in] gid Global mesh node IDs at which we receive grad contributions
916 : : //! \param[in] G1 Partial contributions of extrema for conservative variables to
917 : : //! chare-boundary nodes
918 : : //! \param[in] G2 Partial contributions of extrema for primitive variables to
919 : : //! chare-boundary nodes
920 : : //! \details This function receives contributions to m_uNodalExtrm/m_pNodalExtrm
921 : : //! , which stores nodal extrems at mesh chare-boundary nodes. While
922 : : //! m_uNodalExtrm/m_pNodalExtrm stores own contributions, m_uNodalExtrmc
923 : : //! /m_pNodalExtrmc collects the neighbor chare contributions during
924 : : //! communication.
925 : : // *****************************************************************************
926 : : {
927 : : Assert( G1.size() == gid.size(), "Size mismatch" );
928 : : Assert( G2.size() == gid.size(), "Size mismatch" );
929 : :
930 : 385050 : const auto rdof = g_inputdeck.get< tag::rdof >();
931 : 385050 : const auto ncomp = m_u.nprop() / rdof;
932 : 385050 : const auto nprim = m_p.nprop() / rdof;
933 : :
934 [ + + ]: 3269040 : for (std::size_t i=0; i<gid.size(); ++i)
935 : : {
936 : : auto& u = m_uNodalExtrmc[gid[i]];
937 : 2883990 : auto& p = m_pNodalExtrmc[gid[i]];
938 [ + + ]: 20524560 : for (std::size_t c=0; c<ncomp; ++c)
939 : : {
940 [ + + ]: 70562280 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
941 : : {
942 : 52921710 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
943 : 52921710 : auto min_mark = max_mark + 1;
944 [ + + ][ + + ]: 53469278 : u[max_mark] = std::max( G1[i][max_mark], u[max_mark] );
945 : 52921710 : u[min_mark] = std::min( G1[i][min_mark], u[min_mark] );
946 : : }
947 : : }
948 [ + + ]: 7291440 : for (std::size_t c=0; c<nprim; ++c)
949 : : {
950 [ + + ]: 17629800 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
951 : : {
952 : 13222350 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
953 : 13222350 : auto min_mark = max_mark + 1;
954 [ - + ][ - + ]: 13222350 : p[max_mark] = std::max( G2[i][max_mark], p[max_mark] );
955 : 13222350 : p[min_mark] = std::min( G2[i][min_mark], p[min_mark] );
956 : : }
957 : : }
958 : : }
959 : :
960 [ + + ]: 385050 : if (++m_nnodalExtrema == Disc()->NodeCommMap().size())
961 : : {
962 : 34875 : m_nnodalExtrema = 0;
963 : 34875 : comnodalExtrema_complete();
964 : : }
965 : 385050 : }
966 : :
967 : 38936 : void DG::resizeNodalExtremac()
968 : : // *****************************************************************************
969 : : // Resize the buffer vector of nodal extrema
970 : : // *****************************************************************************
971 : : {
972 : 38936 : const auto rdof = g_inputdeck.get< tag::rdof >();
973 : 38936 : const auto ncomp = m_u.nprop() / rdof;
974 : 38936 : const auto nprim = m_p.nprop() / rdof;
975 : :
976 : 38936 : auto large = std::numeric_limits< tk::real >::max();
977 [ + - ][ + + ]: 469946 : for (const auto& [c,n] : Disc()->NodeCommMap())
978 : : {
979 [ + + ][ + - ]: 3322548 : for (auto i : n) {
980 : : auto& u = m_uNodalExtrmc[i];
981 : : auto& p = m_pNodalExtrmc[i];
982 [ + - ]: 2930474 : u.resize( 2*m_ndof_NodalExtrm*ncomp, large );
983 [ + - ]: 2930474 : p.resize( 2*m_ndof_NodalExtrm*nprim, large );
984 : :
985 : : // Initialize the minimum nodal extrema
986 [ + + ]: 11721896 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
987 : : {
988 [ + + ]: 62506722 : for(std::size_t k = 0; k < ncomp; k++)
989 : 53715300 : u[2*k*m_ndof_NodalExtrm+2*idof] = -large;
990 [ + + ]: 22160034 : for(std::size_t k = 0; k < nprim; k++)
991 : 13368612 : p[2*k*m_ndof_NodalExtrm+2*idof] = -large;
992 : : }
993 : : }
994 : : }
995 : 38936 : }
996 : :
997 : 7110 : void DG::evalNodalExtrmRefEl(
998 : : const std::size_t ncomp,
999 : : const std::size_t nprim,
1000 : : const std::size_t ndof_NodalExtrm,
1001 : : const std::vector< std::size_t >& bndel,
1002 : : const std::vector< std::size_t >& inpoel,
1003 : : const std::vector< std::size_t >& gid,
1004 : : const std::unordered_map< std::size_t, std::size_t >& bid,
1005 : : const tk::Fields& U,
1006 : : const tk::Fields& P,
1007 : : std::vector< std::vector<tk::real> >& uNodalExtrm,
1008 : : std::vector< std::vector<tk::real> >& pNodalExtrm )
1009 : : // *****************************************************************************
1010 : : // Compute the nodal extrema of ref el derivatives for chare-boundary nodes
1011 : : //! \param[in] ncomp Number of conservative variables
1012 : : //! \param[in] nprim Number of primitive variables
1013 : : //! \param[in] ndof_NodalExtrm Degree of freedom for nodal extrema
1014 : : //! \param[in] bndel List of elements contributing to chare-boundary nodes
1015 : : //! \param[in] inpoel Element-node connectivity for element e
1016 : : //! \param[in] gid Local->global node id map
1017 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
1018 : : //! global node ids (key)
1019 : : //! \param[in] U Vector of conservative variables
1020 : : //! \param[in] P Vector of primitive variables
1021 : : //! \param[in,out] uNodalExtrm Chare-boundary nodal extrema for conservative
1022 : : //! variables
1023 : : //! \param[in,out] pNodalExtrm Chare-boundary nodal extrema for primitive
1024 : : //! variables
1025 : : // *****************************************************************************
1026 : : {
1027 : 7110 : const auto rdof = g_inputdeck.get< tag::rdof >();
1028 : :
1029 [ + + ]: 653190 : for (auto e : bndel)
1030 : : {
1031 : : // access node IDs
1032 : : const std::vector<std::size_t> N
1033 [ + - ]: 646080 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
1034 : :
1035 : : // Loop over nodes of element e
1036 [ + + ]: 3230400 : for(std::size_t ip=0; ip<4; ++ip)
1037 : : {
1038 : 2584320 : auto i = bid.find( gid[N[ip]] );
1039 [ + + ]: 2584320 : if (i != end(bid)) // If ip is the chare boundary point
1040 : : {
1041 : : // If DG(P2) is applied, find the nodal extrema of the gradients of
1042 : : // conservative/primitive variables in the reference element
1043 : :
1044 : : // Vector used to store the first order derivatives for both
1045 : : // conservative and primitive variables
1046 [ + - ][ - - ]: 1754730 : std::vector< std::array< tk::real, 3 > > gradc(ncomp, {0.0, 0.0, 0.0});
1047 [ + - ][ - - ]: 1754730 : std::vector< std::array< tk::real, 3 > > gradp(ncomp, {0.0, 0.0, 0.0});
1048 : :
1049 : : // Derivatives of the Dubiner basis
1050 : 1754730 : std::array< tk::real, 3 > center {{0.25, 0.25, 0.25}};
1051 [ + - ]: 1754730 : auto dBdxi = tk::eval_dBdxi(rdof, center);
1052 : :
1053 : : // Evaluate the first order derivative
1054 [ + + ]: 10528380 : for(std::size_t icomp = 0; icomp < ncomp; icomp++)
1055 : : {
1056 : 8773650 : auto mark = icomp * rdof;
1057 [ + + ]: 35094600 : for(std::size_t idir = 0; idir < 3; idir++)
1058 : : {
1059 : 26320950 : gradc[icomp][idir] = 0;
1060 [ + + ]: 263209500 : for(std::size_t idof = 1; idof < rdof; idof++)
1061 : 236888550 : gradc[icomp][idir] += U(e, mark+idof) * dBdxi[idir][idof];
1062 : : }
1063 : : }
1064 [ - + ]: 1754730 : for(std::size_t icomp = 0; icomp < nprim; icomp++)
1065 : : {
1066 : 0 : auto mark = icomp * rdof;
1067 [ - - ]: 0 : for(std::size_t idir = 0; idir < 3; idir++)
1068 : : {
1069 : 0 : gradp[icomp][idir] = 0;
1070 [ - - ]: 0 : for(std::size_t idof = 1; idof < rdof; idof++)
1071 : 0 : gradp[icomp][idir] += P(e, mark+idof) * dBdxi[idir][idof];
1072 : : }
1073 : : }
1074 : :
1075 : : // Store the extrema for the gradients
1076 [ + + ]: 10528380 : for (std::size_t c=0; c<ncomp; ++c)
1077 : : {
1078 [ + + ]: 35094600 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
1079 : : {
1080 : 26320950 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1081 : 26320950 : auto min_mark = max_mark + 1;
1082 [ + + ]: 26320950 : auto& ex = uNodalExtrm[i->second];
1083 [ + + ][ + + ]: 33003900 : ex[max_mark] = std::max(ex[max_mark], gradc[c][idof]);
1084 : 26320950 : ex[min_mark] = std::min(ex[min_mark], gradc[c][idof]);
1085 : : }
1086 : : }
1087 [ - + ]: 1754730 : for (std::size_t c=0; c<nprim; ++c)
1088 : : {
1089 [ - - ]: 0 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
1090 : : {
1091 : 0 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1092 : 0 : auto min_mark = max_mark + 1;
1093 [ - - ]: 0 : auto& ex = pNodalExtrm[i->second];
1094 [ - - ][ - - ]: 0 : ex[max_mark] = std::max(ex[max_mark], gradp[c][idof]);
1095 : 0 : ex[min_mark] = std::min(ex[min_mark], gradp[c][idof]);
1096 : : }
1097 : : }
1098 : : }
1099 : : }
1100 : : }
1101 : 7110 : }
1102 : :
1103 : : void
1104 : 38265 : DG::lim()
1105 : : // *****************************************************************************
1106 : : // Compute limiter function
1107 : : // *****************************************************************************
1108 : : {
1109 : 38265 : auto d = Disc();
1110 : 38265 : const auto rdof = g_inputdeck.get< tag::rdof >();
1111 : 38265 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1112 : 38265 : const auto ncomp = m_u.nprop() / rdof;
1113 : 38265 : const auto nprim = m_p.nprop() / rdof;
1114 : :
1115 : : // Combine own and communicated contributions to nodal extrema
1116 [ + + ]: 1317765 : for (const auto& [gid,g] : m_uNodalExtrmc) {
1117 : 1279500 : auto bid = tk::cref_find( d->Bid(), gid );
1118 [ + + ]: 9304560 : for (ncomp_t c=0; c<ncomp; ++c)
1119 : : {
1120 [ + + ]: 32100240 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1121 : : {
1122 : 24075180 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1123 : 24075180 : auto min_mark = max_mark + 1;
1124 : 24075180 : m_uNodalExtrm[bid][max_mark] =
1125 [ + + ][ + + ]: 24802452 : std::max(g[max_mark], m_uNodalExtrm[bid][max_mark]);
1126 : 24075180 : m_uNodalExtrm[bid][min_mark] =
1127 : 24075180 : std::min(g[min_mark], m_uNodalExtrm[bid][min_mark]);
1128 : : }
1129 : : }
1130 : : }
1131 [ + + ]: 1317765 : for (const auto& [gid,g] : m_pNodalExtrmc) {
1132 : 1279500 : auto bid = tk::cref_find( d->Bid(), gid );
1133 [ + + ]: 3460350 : for (ncomp_t c=0; c<nprim; ++c)
1134 : : {
1135 [ + + ]: 8723400 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1136 : : {
1137 : 6542550 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1138 : 6542550 : auto min_mark = max_mark + 1;
1139 : 6542550 : m_pNodalExtrm[bid][max_mark] =
1140 [ - + ][ - + ]: 6542550 : std::max(g[max_mark], m_pNodalExtrm[bid][max_mark]);
1141 : 6542550 : m_pNodalExtrm[bid][min_mark] =
1142 : 6542550 : std::min(g[min_mark], m_pNodalExtrm[bid][min_mark]);
1143 : : }
1144 : : }
1145 : : }
1146 : :
1147 : : // clear gradients receive buffer
1148 : 38265 : tk::destroy(m_uNodalExtrmc);
1149 : 38265 : tk::destroy(m_pNodalExtrmc);
1150 : :
1151 [ + + ]: 38265 : if (rdof > 1) {
1152 : 49680 : g_dgpde[d->MeshId()].limit( d->T(), pref, myGhosts()->m_geoFace,
1153 : 24840 : myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_esup,
1154 : 24840 : myGhosts()->m_inpoel, myGhosts()->m_coord, m_ndof, d->Gid(),
1155 : 24840 : d->Bid(), m_uNodalExtrm, m_pNodalExtrm, m_mtInv, m_u, m_p,
1156 : 24840 : m_shockmarker );
1157 : :
1158 [ + + ]: 24840 : if (g_inputdeck.get< tag::limsol_projection >())
1159 : 41580 : g_dgpde[d->MeshId()].CPL(m_p, myGhosts()->m_geoElem,
1160 : 20790 : myGhosts()->m_inpoel, myGhosts()->m_coord, m_u,
1161 : 20790 : myGhosts()->m_fd.Esuel().size()/4);
1162 : : }
1163 : :
1164 : : // Send limited solution to neighboring chares
1165 [ + + ]: 38265 : if (myGhosts()->m_sendGhost.empty())
1166 : 3390 : comlim_complete();
1167 : : else
1168 [ + + ]: 454800 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
1169 [ + - ]: 385050 : std::vector< std::size_t > tetid( ghostdata.size() );
1170 [ + - ][ + - ]: 770100 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
1171 [ + - ]: 385050 : prim( ghostdata.size() );
1172 : : std::vector< std::size_t > ndof;
1173 : : std::size_t j = 0;
1174 [ + + ]: 6854730 : for(const auto& i : ghostdata) {
1175 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
1176 : : "Sending limiter ghost data" );
1177 [ + - ]: 6469680 : tetid[j] = i;
1178 [ + - ][ - + ]: 6469680 : u[j] = m_u[i];
1179 [ + - ][ - + ]: 6469680 : prim[j] = m_p[i];
1180 : 6469680 : ++j;
1181 : : }
1182 [ + - ][ + - ]: 770100 : thisProxy[ cid ].comlim( thisIndex, tetid, u, prim );
1183 : : }
1184 : :
1185 : 38265 : ownlim_complete();
1186 : 38265 : }
1187 : :
1188 : : void
1189 : 1770 : DG::refine_ndof()
1190 : : // *****************************************************************************
1191 : : // p-refine all elements that are adjacent to p-refined elements
1192 : : //! \details This function p-refines all the neighbors of an element that has
1193 : : //! been p-refined as a result of an error indicator.
1194 : : // *****************************************************************************
1195 : : {
1196 : 1770 : auto d = Disc();
1197 : : const auto& coord = d->Coord();
1198 : 1770 : const auto& inpoel = d->Inpoel();
1199 : 1770 : const auto npoin = coord[0].size();
1200 : 1770 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
1201 : 1770 : std::vector<std::size_t> node_ndof(npoin, 1);
1202 : :
1203 : : // Mark the max ndof for each node and store in node_ndof
1204 [ + + ]: 188540 : for(std::size_t ip=0; ip<npoin; ip++)
1205 : : {
1206 [ + - ]: 186770 : const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
1207 [ + + ]: 2641810 : for(auto er : pesup)
1208 [ + + ]: 2536953 : node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
1209 : : }
1210 : :
1211 [ + + ]: 414490 : for(std::size_t e = 0; e < nelem; e++)
1212 : : {
1213 : : // Find if any node of this element has p1/p2 ndofs
1214 : : std::size_t counter_p2(0);
1215 : : std::size_t counter_p1(0);
1216 [ + + ]: 2063600 : for(std::size_t inode = 0; inode < 4; inode++)
1217 : : {
1218 [ + + ]: 1650880 : auto node = inpoel[4*e+inode];
1219 [ + + ]: 1650880 : if(node_ndof[node] == 10)
1220 : 320808 : counter_p2++;
1221 [ + + ]: 1330072 : else if (node_ndof[node] == 4)
1222 : 180140 : counter_p1++;
1223 : : }
1224 : :
1225 : : // If there is at least one node with p1/p2 ndofs, all of the elements
1226 : : // around this node are refined to p1/p2.
1227 [ + + ][ + + ]: 412720 : if(counter_p2 > 0 && m_ndof[e] < 10)
1228 : : {
1229 [ + + ]: 16717 : if(m_ndof[e] == 4)
1230 : 15693 : m_ndof[e] = 10;
1231 [ + + ]: 16717 : if(m_ndof[e] == 1)
1232 : 1024 : m_ndof[e] = 4;
1233 : : }
1234 [ + + ][ + + ]: 396003 : else if(counter_p1 > 0 && m_ndof[e] < 4)
1235 : 13452 : m_ndof[e] = 4;
1236 : : }
1237 : 1770 : }
1238 : :
1239 : 1770 : void DG::smooth_ndof()
1240 : : // *****************************************************************************
1241 : : // Smooth the refined ndof distribution to avoid zigzag refinement
1242 : : // *****************************************************************************
1243 : : {
1244 : 1770 : auto d = Disc();
1245 : 1770 : const auto& inpoel = d->Inpoel();
1246 : : const auto& coord = d->Coord();
1247 : 1770 : const auto npoin = coord[0].size();
1248 : 1770 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
1249 : 1770 : std::vector<std::size_t> node_ndof(npoin, 1);
1250 : :
1251 : : // Mark the max ndof for each node and store in node_ndof
1252 [ + + ]: 188540 : for(std::size_t ip=0; ip<npoin; ip++)
1253 : : {
1254 [ + - ]: 186770 : const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
1255 [ + + ]: 2641810 : for(auto er : pesup)
1256 [ + + ]: 2542909 : node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
1257 : : }
1258 : :
1259 [ + + ]: 414490 : for(std::size_t e = 0; e < nelem; e++)
1260 : : {
1261 : : // Find if any node of this element has p1/p2 ndofs
1262 : : std::size_t counter_p2(0);
1263 : : std::size_t counter_p1(0);
1264 [ + + ]: 2063600 : for(std::size_t inode = 0; inode < 4; inode++)
1265 : : {
1266 [ + + ]: 1650880 : auto node = inpoel[4*e+inode];
1267 [ + + ]: 1650880 : if(node_ndof[node] == 10)
1268 : 382503 : counter_p2++;
1269 [ + + ]: 1268377 : else if (node_ndof[node] == 4)
1270 : 182564 : counter_p1++;
1271 : : }
1272 : :
1273 : : // If all the nodes in the element are p1/p2, this element is refined to
1274 : : // p1/p2.
1275 [ + + ][ + + ]: 412720 : if(counter_p2 == 4 && m_ndof[e] == 4)
1276 : 1469 : m_ndof[e] = 10;
1277 [ + + ][ + + ]: 411251 : else if(counter_p1 == 4 && m_ndof[e] == 1)
1278 : 1553 : m_ndof[e] = 4;
1279 : : }
1280 : 1770 : }
1281 : :
1282 : : void
1283 : 385050 : DG::comlim( int fromch,
1284 : : const std::vector< std::size_t >& tetid,
1285 : : const std::vector< std::vector< tk::real > >& u,
1286 : : const std::vector< std::vector< tk::real > >& prim )
1287 : : // *****************************************************************************
1288 : : // Receive chare-boundary limiter ghost data from neighboring chares
1289 : : //! \param[in] fromch Sender chare id
1290 : : //! \param[in] tetid Ghost tet ids we receive solution data for
1291 : : //! \param[in] u Limited high-order solution
1292 : : //! \param[in] prim Limited high-order primitive quantities
1293 : : //! \details This function receives contributions to the limited solution from
1294 : : //! fellow chares.
1295 : : // *****************************************************************************
1296 : : {
1297 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comlim()" );
1298 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comlim()" );
1299 : :
1300 : : // Find local-to-ghost tet id map for sender chare
1301 : 385050 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
1302 : :
1303 [ + + ]: 6854730 : for (std::size_t i=0; i<tetid.size(); ++i) {
1304 : 6469680 : auto j = tk::cref_find( n, tetid[i] );
1305 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
1306 : : "Receiving solution non-ghost data" );
1307 : 6469680 : auto b = tk::cref_find( myGhosts()->m_bid, j );
1308 : : Assert( b < m_uc[2].size(), "Indexing out of bounds" );
1309 : : Assert( b < m_pc[2].size(), "Indexing out of bounds" );
1310 : 6469680 : m_uc[2][b] = u[i];
1311 : 6469680 : m_pc[2][b] = prim[i];
1312 : : }
1313 : :
1314 : : // if we have received all solution ghost contributions from neighboring
1315 : : // chares (chares we communicate along chare-boundary faces with), and
1316 : : // contributed our solution to these neighbors, proceed to limiting
1317 [ + + ]: 385050 : if (++m_nlim == myGhosts()->m_sendGhost.size()) {
1318 : 34875 : m_nlim = 0;
1319 : 34875 : comlim_complete();
1320 : : }
1321 : 385050 : }
1322 : :
1323 : : void
1324 : 38265 : DG::dt()
1325 : : // *****************************************************************************
1326 : : // Compute time step size
1327 : : // *****************************************************************************
1328 : : {
1329 : 38265 : auto d = Disc();
1330 : :
1331 : : // Combine own and communicated contributions of limited solution and degrees
1332 : : // of freedom in cells (if p-adaptive)
1333 [ + + ]: 6546210 : for (const auto& b : myGhosts()->m_bid) {
1334 : : Assert( m_uc[2][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
1335 : : Assert( m_pc[2][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
1336 [ + + ]: 173340045 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1337 : 166870365 : m_u(b.first,c) = m_uc[2][b.second][c];
1338 : : }
1339 [ + + ]: 34324605 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
1340 : 27854925 : m_p(b.first,c) = m_pc[2][b.second][c];
1341 : : }
1342 : : }
1343 : :
1344 : 38265 : auto mindt = std::numeric_limits< tk::real >::max();
1345 : :
1346 [ + + ]: 38265 : if (m_stage == 0)
1347 : : {
1348 [ + + ]: 12755 : auto const_dt = g_inputdeck.get< tag::dt >();
1349 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1350 : :
1351 : : // use constant dt if configured
1352 [ + + ]: 12755 : if (std::abs(const_dt) > eps) {
1353 : :
1354 : 9925 : mindt = const_dt;
1355 : :
1356 : : } else { // compute dt based on CFL
1357 : :
1358 : : // find the minimum dt across all PDEs integrated
1359 : : auto eqdt =
1360 : 5660 : g_dgpde[d->MeshId()].dt( myGhosts()->m_coord, myGhosts()->m_inpoel,
1361 : 2830 : myGhosts()->m_fd,
1362 : 2830 : myGhosts()->m_geoFace, myGhosts()->m_geoElem, m_ndof, m_u, m_p,
1363 : 2830 : myGhosts()->m_fd.Esuel().size()/4 );
1364 [ + - ]: 2830 : if (eqdt < mindt) mindt = eqdt;
1365 : :
1366 : 2830 : mindt *= g_inputdeck.get< tag::cfl >();
1367 : : }
1368 : : }
1369 : : else
1370 : : {
1371 : 25510 : mindt = d->Dt();
1372 : : }
1373 : :
1374 : : // Resize the buffer vector of nodal extrema
1375 : 38265 : resizeNodalExtremac();
1376 : :
1377 : : // Contribute to minimum dt across all chares then advance to next step
1378 [ + - ]: 38265 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
1379 : 38265 : CkCallback(CkReductionTarget(DG,solve), thisProxy) );
1380 : 38265 : }
1381 : :
1382 : : void
1383 : 38265 : DG::solve( tk::real newdt )
1384 : : // *****************************************************************************
1385 : : // Compute right-hand side of discrete transport equations
1386 : : //! \param[in] newdt Size of this new time step
1387 : : // *****************************************************************************
1388 : : {
1389 : : // Enable SDAG wait for building the solution vector during the next stage
1390 [ + - ]: 38265 : thisProxy[ thisIndex ].wait4sol();
1391 [ + - ]: 38265 : thisProxy[ thisIndex ].wait4refine();
1392 [ + - ]: 38265 : thisProxy[ thisIndex ].wait4smooth();
1393 [ + - ]: 38265 : thisProxy[ thisIndex ].wait4reco();
1394 [ + - ]: 38265 : thisProxy[ thisIndex ].wait4nodalExtrema();
1395 [ + - ]: 38265 : thisProxy[ thisIndex ].wait4lim();
1396 [ + - ]: 38265 : thisProxy[ thisIndex ].wait4nod();
1397 : :
1398 : 38265 : auto d = Disc();
1399 : 38265 : const auto rdof = g_inputdeck.get< tag::rdof >();
1400 : 38265 : const auto ndof = g_inputdeck.get< tag::ndof >();
1401 : 38265 : const auto neq = m_u.nprop()/rdof;
1402 : 38265 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1403 : :
1404 : : // Set new time step size
1405 [ + + ]: 38265 : if (m_stage == 0) d->setdt( newdt );
1406 : :
1407 : : // Update Un
1408 [ + + ]: 38265 : if (m_stage == 0) m_un = m_u;
1409 : :
1410 : : // Explicit or IMEX
1411 : 38265 : const auto imex_runge_kutta = g_inputdeck.get< tag::imex_runge_kutta >();
1412 : :
1413 : : // physical time at time-stage for computing exact source terms
1414 : 38265 : tk::real physT(d->T());
1415 [ + + ]: 38265 : if (m_stage == 1) {
1416 : 12755 : physT += d->Dt();
1417 : : }
1418 [ + + ]: 25510 : else if (m_stage == 2) {
1419 : 12755 : physT += 0.5*d->Dt();
1420 : : }
1421 : :
1422 [ - + ]: 38265 : if (imex_runge_kutta) {
1423 [ - - ]: 0 : if (m_stage == 0)
1424 : : {
1425 : : // Save previous rhs
1426 : : m_rhsprev = m_rhs;
1427 : : // Initialize m_stiffrhs to zero
1428 : : m_stiffrhs.fill(0.0);
1429 : : m_stiffrhsprev.fill(0.0);
1430 : : }
1431 : : }
1432 : :
1433 : 76530 : g_dgpde[d->MeshId()].rhs( physT, pref, myGhosts()->m_geoFace,
1434 : 38265 : myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_inpoel, m_boxelems,
1435 : 38265 : myGhosts()->m_coord, m_u, m_p, m_ndof, d->Dt(), m_rhs );
1436 : :
1437 [ + - ]: 38265 : if (!imex_runge_kutta) {
1438 : : // Explicit time-stepping using RK3 to discretize time-derivative
1439 [ + + ]: 16351845 : for(std::size_t e=0; e<myGhosts()->m_nunk; ++e)
1440 [ + + ]: 111412515 : for(std::size_t c=0; c<neq; ++c)
1441 : : {
1442 [ + + ]: 373799400 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
1443 : : {
1444 [ + + ]: 278700465 : if(k < m_ndof[e]) {
1445 : 193948005 : auto rmark = c*rdof+k;
1446 : 193948005 : auto mark = c*ndof+k;
1447 [ + + ]: 193948005 : m_u(e, rmark) = rkcoef[0][m_stage] * m_un(e, rmark)
1448 : 193948005 : + rkcoef[1][m_stage] * ( m_u(e, rmark)
1449 : 193948005 : + d->Dt() * m_rhs(e, mark)/m_lhs(e, mark) );
1450 [ + + ]: 193948005 : if(fabs(m_u(e, rmark)) < 1e-16)
1451 : 31378053 : m_u(e, rmark) = 0;
1452 : : }
1453 : : }
1454 : : }
1455 : : }
1456 : : else {
1457 : : // Implicit-Explicit time-stepping using RK3 to discretize time-derivative
1458 : 0 : DG::imex_integrate();
1459 : : }
1460 : :
1461 [ + + ]: 16351845 : for(std::size_t e=0; e<myGhosts()->m_nunk; ++e)
1462 [ + + ]: 111412515 : for(std::size_t c=0; c<neq; ++c)
1463 : : {
1464 : : // zero out unused/reconstructed dofs of equations using reduced dofs
1465 : : // (see DGMultiMat::numEquationDofs())
1466 [ + + ]: 95098935 : if (m_numEqDof[c] < rdof) {
1467 [ + + ]: 96926640 : for (std::size_t k=m_numEqDof[c]; k<rdof; ++k)
1468 : : {
1469 : 72694980 : auto rmark = c*rdof+k;
1470 : 72694980 : m_u(e, rmark) = 0.0;
1471 : : }
1472 : : }
1473 : : }
1474 : :
1475 : : // Update primitives based on the evolved solution
1476 : 38265 : g_dgpde[d->MeshId()].updateInterfaceCells( m_u,
1477 : 38265 : myGhosts()->m_fd.Esuel().size()/4, m_ndof, m_interface );
1478 : 38265 : g_dgpde[d->MeshId()].updatePrimitives( m_u, m_lhs, myGhosts()->m_geoElem, m_p,
1479 : 38265 : myGhosts()->m_fd.Esuel().size()/4, m_ndof );
1480 [ + - ]: 38265 : if (!g_inputdeck.get< tag::accuracy_test >()) {
1481 : 38265 : g_dgpde[d->MeshId()].cleanTraceMaterial( physT, myGhosts()->m_geoElem, m_u,
1482 : 38265 : m_p, myGhosts()->m_fd.Esuel().size()/4 );
1483 : : }
1484 : :
1485 [ + + ]: 38265 : if (m_stage < m_nstage-1) {
1486 : :
1487 : : // continue with next time step stage
1488 : 25510 : stage();
1489 : :
1490 : : } else {
1491 : :
1492 : : // Increase number of iterations and physical time
1493 : 12755 : d->next();
1494 : :
1495 : : // Compute diagnostics, e.g., residuals
1496 : 12755 : auto diag_computed = m_diag.compute( *d,
1497 : 12755 : m_u.nunk()-myGhosts()->m_fd.Esuel().size()/4, myGhosts()->m_geoElem,
1498 : 12755 : m_ndof, m_u, m_un );
1499 : :
1500 : : // Continue to mesh refinement (if configured)
1501 [ + + ][ + - ]: 15117 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 0.0 ) );
[ + - ]
1502 : :
1503 : : }
1504 : 38265 : }
1505 : :
1506 : : void
1507 : 12755 : DG::refine( [[maybe_unused]] const std::vector< tk::real >& l2res )
1508 : : // *****************************************************************************
1509 : : // Optionally refine/derefine mesh
1510 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1511 : : //! computed across the whole problem
1512 : : // *****************************************************************************
1513 : : {
1514 : 12755 : auto d = Disc();
1515 : :
1516 : 12755 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
1517 : 12755 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
1518 : :
1519 : : // if t>0 refinement enabled and we hit the dtref frequency
1520 [ - + ][ - - ]: 12755 : if (dtref && !(d->It() % dtfreq)) { // refine
1521 : :
1522 : 0 : d->startvol();
1523 [ - - ][ - - ]: 0 : d->Ref()->dtref( myGhosts()->m_fd.Bface(), {},
[ - - ][ - - ]
1524 : 0 : tk::remap(myGhosts()->m_fd.Triinpoel(),d->Gid()) );
1525 : 0 : d->refined() = 1;
1526 : :
1527 : : } else { // do not refine
1528 : :
1529 : 12755 : d->refined() = 0;
1530 : 12755 : stage();
1531 : :
1532 : : }
1533 : 12755 : }
1534 : :
1535 : : void
1536 : 0 : DG::resizePostAMR(
1537 : : const std::vector< std::size_t >& /*ginpoel*/,
1538 : : const tk::UnsMesh::Chunk& chunk,
1539 : : const tk::UnsMesh::Coords& coord,
1540 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
1541 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
1542 : : const std::set< std::size_t >& removedNodes,
1543 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
1544 : : const tk::NodeCommMap& nodeCommMap,
1545 : : const std::map< int, std::vector< std::size_t > >& bface,
1546 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
1547 : : const std::vector< std::size_t >& triinpoel,
1548 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid )
1549 : : // *****************************************************************************
1550 : : // Receive new mesh from Refiner
1551 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1552 : : //! \param[in] coord New mesh node coordinates
1553 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1554 : : //! \param[in] removedNodes Newly removed mesh node local ids
1555 : : //! \param[in] amrNodeMap Node id map after amr (local ids)
1556 : : //! \param[in] nodeCommMap New node communication map
1557 : : //! \param[in] bface Boundary-faces mapped to side set ids
1558 : : //! \param[in] triinpoel Boundary-face connectivity
1559 : : //! \param[in] elemblockid Local tet ids associated with mesh block ids
1560 : : // *****************************************************************************
1561 : : {
1562 : 0 : auto d = Disc();
1563 : :
1564 : : // Set flag that indicates that we are during time stepping
1565 : 0 : m_initial = 0;
1566 : 0 : myGhosts()->m_initial = 0;
1567 : :
1568 : : // Zero field output iteration count between two mesh refinement steps
1569 : 0 : d->Itf() = 0;
1570 : :
1571 : : // Increase number of iterations with mesh refinement
1572 : 0 : ++d->Itr();
1573 : :
1574 : : // Save old number of elements
1575 : 0 : [[maybe_unused]] auto old_nelem = myGhosts()->m_inpoel.size()/4;
1576 : :
1577 : : // Resize mesh data structures
1578 : 0 : d->resizePostAMR( chunk, coord, amrNodeMap, nodeCommMap, removedNodes,
1579 : : elemblockid );
1580 : :
1581 : : // Update state
1582 : 0 : myGhosts()->m_inpoel = d->Inpoel();
1583 : 0 : myGhosts()->m_coord = d->Coord();
1584 : 0 : auto nelem = myGhosts()->m_inpoel.size()/4;
1585 : : m_p.resize( nelem );
1586 : : m_u.resize( nelem );
1587 : : m_un.resize( nelem );
1588 : : m_lhs.resize( nelem );
1589 : : m_rhs.resize( nelem );
1590 : : m_rhsprev.resize( nelem );
1591 : : m_stiffrhs.resize( nelem );
1592 : : m_stiffrhsprev.resize( nelem );
1593 [ - - ][ - - ]: 0 : m_uNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
1594 [ - - ]: 0 : m_ndof_NodalExtrm*g_inputdeck.get< tag::ncomp >() ) );
1595 [ - - ][ - - ]: 0 : m_pNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
1596 [ - - ]: 0 : m_ndof_NodalExtrm*m_p.nprop()/g_inputdeck.get< tag::rdof >()));
1597 : :
1598 : : // Resize the buffer vector of nodal extrema
1599 : 0 : resizeNodalExtremac();
1600 : :
1601 [ - - ][ - - ]: 0 : myGhosts()->m_fd = FaceData( myGhosts()->m_inpoel, bface,
[ - - ][ - - ]
[ - - ]
1602 : 0 : tk::remap(triinpoel,d->Lid()) );
1603 : :
1604 [ - - ]: 0 : myGhosts()->m_geoFace =
1605 : 0 : tk::Fields( tk::genGeoFaceTri( myGhosts()->m_fd.Nipfac(),
1606 : 0 : myGhosts()->m_fd.Inpofa(), coord ) );
1607 [ - - ]: 0 : myGhosts()->m_geoElem = tk::Fields( tk::genGeoElemTet( myGhosts()->m_inpoel,
1608 : 0 : coord ) );
1609 : :
1610 : 0 : myGhosts()->m_nfac = myGhosts()->m_fd.Inpofa().size()/3;
1611 : 0 : myGhosts()->m_nunk = nelem;
1612 : 0 : m_npoin = coord[0].size();
1613 : 0 : myGhosts()->m_bndFace.clear();
1614 : 0 : myGhosts()->m_exptGhost.clear();
1615 : 0 : myGhosts()->m_sendGhost.clear();
1616 : 0 : myGhosts()->m_ghost.clear();
1617 : 0 : myGhosts()->m_esup.clear();
1618 : :
1619 : : // Update solution on new mesh, P0 (cell center value) only for now
1620 : : m_un = m_u;
1621 : : auto pn = m_p;
1622 : 0 : auto unprop = m_u.nprop();
1623 : : auto pnprop = m_p.nprop();
1624 [ - - ]: 0 : for (const auto& [child,parent] : addedTets) {
1625 : : Assert( child < nelem, "Indexing out of new solution vector" );
1626 : : Assert( parent < old_nelem, "Indexing out of old solution vector" );
1627 [ - - ]: 0 : for (std::size_t i=0; i<unprop; ++i) m_u(child,i) = m_un(parent,i);
1628 [ - - ]: 0 : for (std::size_t i=0; i<pnprop; ++i) m_p(child,i) = pn(parent,i);
1629 : : }
1630 : : m_un = m_u;
1631 : :
1632 : : // Resize communication buffers
1633 [ - - ][ - - ]: 0 : m_ghosts[thisIndex].resizeComm();
[ - - ][ - - ]
1634 : 0 : }
1635 : :
1636 : : bool
1637 : 10214 : DG::fieldOutput() const
1638 : : // *****************************************************************************
1639 : : // Decide wether to output field data
1640 : : //! \return True if field data is output in this step
1641 : : // *****************************************************************************
1642 : : {
1643 : 10214 : auto d = Disc();
1644 : :
1645 : : // Output field data
1646 [ + + ][ + - ]: 10214 : return d->fielditer() or d->fieldtime() or d->fieldrange() or d->finished();
[ + - ][ + + ]
1647 : : }
1648 : :
1649 : : bool
1650 : 1325 : DG::refinedOutput() const
1651 : : // *****************************************************************************
1652 : : // Decide if we write field output using a refined mesh
1653 : : //! \return True if field output will use a refined mesh
1654 : : // *****************************************************************************
1655 : : {
1656 [ + + ]: 1325 : return g_inputdeck.get< tag::field_output, tag::refined >() &&
1657 [ - + ]: 33 : g_inputdeck.get< tag::scheme >() != ctr::SchemeType::DGP0;
1658 : : }
1659 : :
1660 : : void
1661 : 1325 : DG::writeFields(
1662 : : CkCallback c,
1663 : : const std::unordered_map< std::size_t, std::size_t >& addedTets )
1664 : : // *****************************************************************************
1665 : : // Output mesh field data
1666 : : //! \param[in] c Function to continue with after the write
1667 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1668 : : // *****************************************************************************
1669 : : {
1670 : 1325 : auto d = Disc();
1671 : :
1672 : : const auto& inpoel = std::get< 0 >( m_outmesh.chunk );
1673 : 2650 : auto esup = tk::genEsup( inpoel, 4 );
1674 : 1325 : auto nelem = inpoel.size() / 4;
1675 : :
1676 : : // Combine own and communicated contributions and finish averaging of node
1677 : : // field output in chare boundary nodes
1678 : : const auto& lid = std::get< 2 >( m_outmesh.chunk );
1679 [ + + ]: 51485 : for (const auto& [g,f] : m_uNodefieldsc) {
1680 : : Assert( m_uNodefields.nprop() == f.first.size(), "Size mismatch" );
1681 : 50160 : auto p = tk::cref_find( lid, g );
1682 [ + + ]: 319474 : for (std::size_t i=0; i<f.first.size(); ++i) {
1683 : 269314 : m_uNodefields(p,i) += f.first[i];
1684 : 269314 : m_uNodefields(p,i) /= static_cast< tk::real >(
1685 : 269314 : esup.second[p+1] - esup.second[p] + f.second );
1686 : : }
1687 : : }
1688 : 1325 : tk::destroy( m_uNodefieldsc );
1689 [ + + ]: 51485 : for (const auto& [g,f] : m_pNodefieldsc) {
1690 : : Assert( m_pNodefields.nprop() == f.first.size(), "Size mismatch" );
1691 : 50160 : auto p = tk::cref_find( lid, g );
1692 [ + + ]: 100760 : for (std::size_t i=0; i<f.first.size(); ++i) {
1693 : 50600 : m_pNodefields(p,i) += f.first[i];
1694 : 50600 : m_pNodefields(p,i) /= static_cast< tk::real >(
1695 : 50600 : esup.second[p+1] - esup.second[p] + f.second );
1696 : : }
1697 : : }
1698 : 1325 : tk::destroy( m_pNodefieldsc );
1699 : :
1700 : : // Lambda to decide if a node (global id) is on a chare boundary of the field
1701 : : // output mesh. p - global node id, return true if node is on the chare
1702 : : // boundary.
1703 : : auto chbnd = [ this ]( std::size_t p ) {
1704 : : return
1705 : : std::any_of( m_outmesh.nodeCommMap.cbegin(), m_outmesh.nodeCommMap.cend(),
1706 : : [&](const auto& s) { return s.second.find(p) != s.second.cend(); } );
1707 : : };
1708 : :
1709 : : // Finish computing node field output averages in internal nodes
1710 : 1325 : auto npoin = m_outmesh.coord[0].size();
1711 : : auto& gid = std::get< 1 >( m_outmesh.chunk );
1712 [ + + ]: 264258 : for (std::size_t p=0; p<npoin; ++p) {
1713 [ + + ]: 262933 : if (!chbnd(gid[p])) {
1714 : 212773 : auto n = static_cast< tk::real >( esup.second[p+1] - esup.second[p] );
1715 [ + + ]: 1060697 : for (std::size_t i=0; i<m_uNodefields.nprop(); ++i)
1716 : 847924 : m_uNodefields(p,i) /= n;
1717 [ + + ]: 320807 : for (std::size_t i=0; i<m_pNodefields.nprop(); ++i)
1718 : 108034 : m_pNodefields(p,i) /= n;
1719 : : }
1720 : : }
1721 : :
1722 : : // Collect field output from numerical solution requested by user
1723 : 1325 : auto elemfields = numericFieldOutput( m_uElemfields, tk::Centering::ELEM,
1724 [ + - ][ + - ]: 2650 : g_dgpde[Disc()->MeshId()].OutVarFn(), m_pElemfields );
[ + - ]
1725 : 1325 : auto nodefields = numericFieldOutput( m_uNodefields, tk::Centering::NODE,
1726 [ + - ][ + - ]: 2650 : g_dgpde[Disc()->MeshId()].OutVarFn(), m_pNodefields );
[ + - ]
1727 : :
1728 : : // Collect field output from analytical solutions (if exist)
1729 : 1325 : const auto& coord = m_outmesh.coord;
1730 [ + - ]: 1325 : auto geoElem = tk::genGeoElemTet( inpoel, coord );
1731 [ + - ]: 1325 : auto t = Disc()->T();
1732 [ + - ]: 1325 : analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::ELEM,
1733 [ + - ][ + - ]: 5300 : geoElem.extract_comp(1), geoElem.extract_comp(2), geoElem.extract_comp(3),
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
1734 : : t, elemfields );
1735 [ + - ]: 1325 : analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::NODE, coord[0],
1736 : : coord[1], coord[2], t, nodefields );
1737 : :
1738 : : // Add adaptive indicator array to element-centered field output
1739 [ + + ]: 1325 : if (g_inputdeck.get< tag::pref, tag::pref >()) {
1740 [ + - ]: 402 : std::vector< tk::real > ndof( begin(m_ndof), end(m_ndof) );
1741 [ + - ]: 402 : ndof.resize( nelem );
1742 [ + + ]: 182529 : for(std::size_t k = 0; k < nelem; k++) {
1743 : : // Mark the cell with THINC reconstruction as 0 for output
1744 [ + + ]: 182127 : if(m_interface[k] == 1) ndof[k] = 0;
1745 : : }
1746 [ + + ]: 87834 : for (const auto& [child,parent] : addedTets)
1747 : 87432 : ndof[child] = static_cast< tk::real >( m_ndof[parent] );
1748 [ + - ]: 402 : elemfields.push_back( ndof );
1749 : : }
1750 : :
1751 : : // Add shock detection marker array to element-centered field output
1752 [ + - ][ - - ]: 1325 : std::vector< tk::real > shockmarker( begin(m_shockmarker), end(m_shockmarker) );
1753 : : // Here m_shockmarker has a size of m_u.nunk() which is the number of the
1754 : : // elements within this partition (nelem) plus the ghost partition cells. In
1755 : : // terms of output purpose, we only need the solution data within this
1756 : : // partition. Therefore, resizing it to nelem removes the extra partition
1757 : : // boundary allocations in the shockmarker vector. Since the code assumes that
1758 : : // the boundary elements are on the top, the resize operation keeps the lower
1759 : : // portion.
1760 [ + - ]: 1325 : shockmarker.resize( nelem );
1761 [ + + ]: 158837 : for (const auto& [child,parent] : addedTets)
1762 : 157512 : shockmarker[child] = static_cast< tk::real >(m_shockmarker[parent]);
1763 [ + - ]: 1325 : elemfields.push_back( shockmarker );
1764 : :
1765 : : // Add rho0*det(g)/rho to make sure it is staying close to 1,
1766 : : // averaged for all materials
1767 [ + - ][ - - ]: 1325 : std::vector< tk::real > densityConstr(nelem);
1768 [ + - ]: 1325 : g_dgpde[d->MeshId()].computeDensityConstr(nelem, m_u, densityConstr);
1769 [ + + ]: 158837 : for (const auto& [child,parent] : addedTets)
1770 : 157512 : densityConstr[child] = 0.0;
1771 [ + + ][ + - ]: 1325 : if (densityConstr.size() > 0) elemfields.push_back( densityConstr );
1772 : :
1773 : : // Query fields names requested by user
1774 [ + - ]: 2650 : auto elemfieldnames = numericFieldNames( tk::Centering::ELEM );
1775 [ + - ]: 1325 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
1776 : :
1777 : : // Collect field output names for analytical solutions
1778 [ + - ]: 1325 : analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::ELEM, elemfieldnames );
1779 [ + - ]: 1325 : analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::NODE, nodefieldnames );
1780 : :
1781 [ + + ]: 1325 : if (g_inputdeck.get< tag::pref, tag::pref >()) {
1782 [ + - ]: 804 : elemfieldnames.push_back( "NDOF" );
1783 : : }
1784 : :
1785 [ + - ]: 1325 : elemfieldnames.push_back( "shock_marker" );
1786 : :
1787 [ + + ]: 1325 : if (densityConstr.size() > 0)
1788 [ + - ]: 358 : elemfieldnames.push_back( "density_constraint" );
1789 : :
1790 : : Assert( elemfieldnames.size() == elemfields.size(), "Size mismatch" );
1791 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
1792 : :
1793 : : // Output chare mesh and fields metadata to file
1794 : 1325 : const auto& triinpoel = m_outmesh.triinpoel;
1795 [ + - ][ + - ]: 3975 : d->write( inpoel, m_outmesh.coord, m_outmesh.bface, {},
[ + + ][ - - ]
1796 [ + - ]: 2650 : tk::remap( triinpoel, lid ), elemfieldnames, nodefieldnames,
1797 : : {}, {}, elemfields, nodefields, {}, {}, c );
1798 : 1325 : }
1799 : :
1800 : : void
1801 : 10380 : DG::comnodeout( const std::vector< std::size_t >& gid,
1802 : : const std::vector< std::size_t >& nesup,
1803 : : const std::vector< std::vector< tk::real > >& Lu,
1804 : : const std::vector< std::vector< tk::real > >& Lp )
1805 : : // *****************************************************************************
1806 : : // Receive chare-boundary nodal solution (for field output) contributions from
1807 : : // neighboring chares
1808 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
1809 : : //! \param[in] nesup Number of elements surrounding points
1810 : : //! \param[in] Lu Partial contributions of solution nodal fields to
1811 : : //! chare-boundary nodes
1812 : : //! \param[in] Lp Partial contributions of primitive quantity nodal fields to
1813 : : //! chare-boundary nodes
1814 : : // *****************************************************************************
1815 : : {
1816 : : Assert( gid.size() == nesup.size(), "Size mismatch" );
1817 : : Assert(Lu.size() == m_uNodefields.nprop(), "Fields size mismatch");
1818 : : Assert(Lp.size() == m_pNodefields.nprop(), "Fields size mismatch");
1819 [ + + ]: 69172 : for (std::size_t f=0; f<Lu.size(); ++f)
1820 : : Assert( gid.size() == Lu[f].size(), "Size mismatch" );
1821 [ + + ]: 21008 : for (std::size_t f=0; f<Lp.size(); ++f)
1822 : : Assert( gid.size() == Lp[f].size(), "Size mismatch" );
1823 : :
1824 [ + + ]: 107480 : for (std::size_t i=0; i<gid.size(); ++i) {
1825 : : auto& nfu = m_uNodefieldsc[ gid[i] ];
1826 : 97100 : nfu.first.resize( Lu.size() );
1827 [ + + ]: 641666 : for (std::size_t f=0; f<Lu.size(); ++f) nfu.first[f] += Lu[f][i];
1828 : 97100 : nfu.second += nesup[i];
1829 : 97100 : auto& nfp = m_pNodefieldsc[ gid[i] ];
1830 : 97100 : nfp.first.resize( Lp.size() );
1831 [ + + ]: 206580 : for (std::size_t f=0; f<Lp.size(); ++f) nfp.first[f] += Lp[f][i];
1832 : 97100 : nfp.second += nesup[i];
1833 : : }
1834 : :
1835 : : // When we have heard from all chares we communicate with, this chare is done
1836 [ + + ]: 10380 : if (++m_nnod == Disc()->NodeCommMap().size()) {
1837 : 1164 : m_nnod = 0;
1838 : 1164 : comnodeout_complete();
1839 : : }
1840 : 10380 : }
1841 : :
1842 : : void
1843 : 38265 : DG::stage()
1844 : : // *****************************************************************************
1845 : : // Evaluate whether to continue with next time step stage
1846 : : // *****************************************************************************
1847 : : {
1848 : : // Increment Runge-Kutta stage counter
1849 : 38265 : ++m_stage;
1850 : :
1851 : : // if not all Runge-Kutta stages complete, continue to next time stage,
1852 : : // otherwise prepare for nodal field output
1853 [ + + ]: 38265 : if (m_stage < m_nstage)
1854 : 25510 : next();
1855 : : else
1856 [ + - ][ + - ]: 38265 : startFieldOutput( CkCallback(CkIndex_DG::step(), thisProxy[thisIndex]) );
[ - + ][ - - ]
1857 : 38265 : }
1858 : :
1859 : : void
1860 : 12084 : DG::evalLB( int nrestart )
1861 : : // *****************************************************************************
1862 : : // Evaluate whether to do load balancing
1863 : : //! \param[in] nrestart Number of times restarted
1864 : : // *****************************************************************************
1865 : : {
1866 : 12084 : auto d = Disc();
1867 : :
1868 : : // Detect if just returned from a checkpoint and if so, zero timers
1869 : 12084 : d->restarted( nrestart );
1870 : :
1871 : 12084 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1872 : 12084 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
1873 : :
1874 : : // Load balancing if user frequency is reached or after the second time-step
1875 [ - + ][ - - ]: 12084 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
1876 : :
1877 : 12084 : AtSync();
1878 [ - + ]: 12084 : if (nonblocking) next();
1879 : :
1880 : : } else {
1881 : :
1882 : 0 : next();
1883 : :
1884 : : }
1885 : 12084 : }
1886 : :
1887 : : void
1888 : 12084 : DG::evalRestart()
1889 : : // *****************************************************************************
1890 : : // Evaluate whether to save checkpoint/restart
1891 : : // *****************************************************************************
1892 : : {
1893 : 12084 : auto d = Disc();
1894 : :
1895 : 12084 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
1896 : 12084 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1897 : :
1898 [ + + ][ - + ]: 12084 : if (not benchmark and not (d->It() % rsfreq)) {
1899 : :
1900 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1901 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1902 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
[ - - ]
1903 : :
1904 : : } else {
1905 : :
1906 : 12084 : evalLB( /* nrestart = */ -1 );
1907 : :
1908 : : }
1909 : 12084 : }
1910 : :
1911 : : void
1912 : 12755 : DG::step()
1913 : : // *****************************************************************************
1914 : : // Evaluate wether to continue with next time step
1915 : : // *****************************************************************************
1916 : : {
1917 : 12755 : auto d = Disc();
1918 : :
1919 : : // Output time history
1920 [ + - ][ + - ]: 12755 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
1921 : 0 : std::vector< std::vector< tk::real > > hist;
1922 [ - - ][ - - ]: 0 : auto h = g_dgpde[d->MeshId()].histOutput( d->Hist(), myGhosts()->m_inpoel,
1923 [ - - ][ - - ]: 0 : myGhosts()->m_coord, m_u, m_p );
1924 [ - - ]: 0 : hist.insert( end(hist), begin(h), end(h) );
1925 [ - - ]: 0 : d->history( std::move(hist) );
1926 : : }
1927 : :
1928 : : // Free memory storing output mesh
1929 : 12755 : m_outmesh.destroy();
1930 : :
1931 : : // Output one-liner status report to screen
1932 : 12755 : d->status();
1933 : : // Reset Runge-Kutta stage counter
1934 : 12755 : m_stage = 0;
1935 : :
1936 : 12755 : const auto term = g_inputdeck.get< tag::term >();
1937 : 12755 : const auto nstep = g_inputdeck.get< tag::nstep >();
1938 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
1939 : :
1940 : : // If neither max iterations nor max time reached, continue, otherwise finish
1941 [ + - ][ + + ]: 12755 : if (std::fabs(d->T()-term) > eps && d->It() < nstep) {
1942 : :
1943 : 12084 : evalRestart();
1944 : :
1945 : : } else {
1946 : :
1947 : 671 : auto meshid = d->MeshId();
1948 [ + - ]: 1342 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1949 : 1342 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1950 : :
1951 : : }
1952 : 12755 : }
1953 : :
1954 : : void
1955 : 0 : DG::imex_integrate()
1956 : : // *****************************************************************************
1957 : : // Perform the Implicit-Explicit Runge-Kutta stage update
1958 : : //
1959 : : //! \details
1960 : : //! Performs the Implicit-Explicit Runge-Kutta step. Scheme taken from
1961 : : //! Cavaglieri, D., & Bewley, T. (2015). Low-storage implicit/explicit
1962 : : //! Runge–Kutta schemes for the simulation of stiff high-dimensional ODE
1963 : : //! systems. Journal of Computational Physics, 286, 172-193.
1964 : : //!
1965 : : //! Scheme given by equations (25a,b):
1966 : : //!
1967 : : //! u[0] = u[n] + dt * (expl_rkcoef[1,0]*R_ex(u[n])+impl_rkcoef[1,1]*R_im(u[0]))
1968 : : //!
1969 : : //! u[1] = u[n] + dt * (expl_rkcoef[2,1]*R_ex(u[0])+impl_rkcoef[2,1]*R_im(u[0])
1970 : : //! +impl_rkcoef[2,2]*R_im(u[1]))
1971 : : //!
1972 : : //! u[n+1] = u[n] + dt * (expl_rkcoef[3,1]*R_ex(u[0])+impl_rkcoef[3,1]*R_im(u[0])
1973 : : //! expl_rkcoef[3,2]*R_ex(u[1])+impl_rkcoef[3,2]*R_im(u[1]))
1974 : : //!
1975 : : //! In order to solve the first two equations we need to solve a series of
1976 : : //! systems of non-linear equations:
1977 : : //!
1978 : : //! F1(u[0]) = B1 + R1(u[0]) = 0, and
1979 : : //! F2(u[1]) = B2 + R2(u[1]) = 0,
1980 : : //!
1981 : : //! where
1982 : : //!
1983 : : //! B1 = u[n] + dt * expl_rkcoef[1,0]*R_ex(u[n]),
1984 : : //! R1 = dt * impl_rkcoef[1,1]*R_im(u[0]) - u([0]),
1985 : : //! B2 = u[n] + dt * (expl_rkcoef[2,1]*R_ex(u[0])+impl_rkcoef[2,1]*R_im(u[0])),
1986 : : //! R2 = dt * impl_rkcoef[2,2]*R_im(u[1]) - u([1]).
1987 : : //!
1988 : : //! In order to solve the non-linear system F(U) = 0, we employ Broyden's method.
1989 : : //! Taken from https://en.wikipedia.org/wiki/Broyden%27s_method.
1990 : : //! The method consists in obtaining an approximation for the inverse of the
1991 : : //! Jacobian H = J^(-1) and advancing in a quasi-newton step:
1992 : : //!
1993 : : //! U[k+1] = U[k] - H[k]*F(U[k]),
1994 : : //!
1995 : : //! until F(U) is close enough to zero.
1996 : : //!
1997 : : //! The approximation H[k] is improved at every iteration following
1998 : : //!
1999 : : //! H[k] = H[k-1] + (DU[k]-H[k-1]*DF[k])/(DU[k]^T*H[k-1]*DF[k]) * DU[k]^T*H[k-1],
2000 : : //!
2001 : : //! where DU[k] = U[k] - U[k-1] and DF[k] = F(U[k]) - F(U[k-1)).
2002 : : //!
2003 : : //! This function performs the following main algorithmic steps:
2004 : : //! - If stage == 0 or stage == 1:
2005 : : //! - Take Initial value:
2006 : : //! U[0] = U[n] + dt * expl_rkcoef[1,0]*R_ex(U[n]) (for stage 0)
2007 : : //! U[1] = U[n] + dt * (expl_rkcoef[2,1]*R_ex(U[0])
2008 : : //! +impl_rkcoef[2,1]*R_im(U[0])) (for stage 1)
2009 : : //! - Loop over the Elements (e++)
2010 : : //! - Initialize Jacobian inverse approximation as the identity
2011 : : //! - Compute implicit right-hand-side (F_im) with current U
2012 : : //! - Iterate for the solution (iter++)
2013 : : //! - Compute new solution U[k+1] = U[k] - H[k]*F(U[k])
2014 : : //! - Compute implicit right-hand-side (F_im) with current U
2015 : : //! - Compute DU and DF
2016 : : //! - Update inverse Jacobian approximation by:
2017 : : //! - Compute V1 = H[k-1]*DF[k] and V2 = DU[k]^T*H[k-1]
2018 : : //! - Compute d = DU[k]^T*V1 and V3 = DU[k]-V1
2019 : : //! - Compute V4 = V3/d
2020 : : //! - Update H[k] = H[k-1] + V4*V2
2021 : : //! - Save old U and F
2022 : : //! - Compute absolute and relative errors
2023 : : //! - Break iterations if error < tol or iter == max_iter
2024 : : //! - Update explicit equations using only the explicit terms.
2025 : : //! - Else if stage == 2:
2026 : : //! - Update explicit equations using only the explicit terms.
2027 : : //! - Update implicit equations using:
2028 : : //! u[n+1] = u[n]+dt*(expl_rkcoef[3,1]*R_ex(u[0])+impl_rkcoef[3,1]*R_im(u[0])
2029 : : //! expl_rkcoef[3,2]*R_ex(u[1])+impl_rkcoef[3,2]*R_im(u[1]))
2030 : : // *****************************************************************************
2031 : : {
2032 : 0 : auto d = Disc();
2033 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
2034 : 0 : const auto ndof = g_inputdeck.get< tag::ndof >();
2035 [ - - ]: 0 : if (m_stage < m_nstage-1) {
2036 : : // Save previous stiff_rhs
2037 : : m_stiffrhsprev = m_stiffrhs;
2038 : :
2039 : : // Compute the imex update
2040 : :
2041 : : // Integrate explicitly on the imex equations
2042 : : // (To use as initial values)
2043 [ - - ]: 0 : for (std::size_t e=0; e<myGhosts()->m_nunk; ++e)
2044 [ - - ]: 0 : for (std::size_t c=0; c<m_nstiffeq; ++c)
2045 : : {
2046 [ - - ]: 0 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
2047 : : {
2048 [ - - ]: 0 : auto rmark = m_stiffEqIdx[c]*rdof+k;
2049 [ - - ]: 0 : auto mark = m_stiffEqIdx[c]*ndof+k;
2050 : 0 : m_u(e, rmark) = m_un(e, rmark) + d->Dt() * (
2051 : 0 : expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/m_lhs(e, mark)
2052 : 0 : + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/m_lhs(e, mark)
2053 : 0 : + impl_rkcoef[0][m_stage]
2054 : 0 : * m_stiffrhsprev(e,c*ndof+k)/m_lhs(e, mark) );
2055 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2056 : 0 : m_u(e, rmark) = 0;
2057 : : }
2058 : : }
2059 : :
2060 : : // Solve for implicit-explicit equations
2061 : 0 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
2062 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
2063 : : {
2064 : : // Non-linear system f(u) = 0 to be solved
2065 : : // Broyden's method
2066 : : // Control parameters
2067 : 0 : std::size_t max_iter = g_inputdeck.get< tag::imex_maxiter >();
2068 : 0 : tk::real rel_tol = g_inputdeck.get< tag::imex_reltol >();
2069 : 0 : tk::real abs_tol = g_inputdeck.get< tag::imex_abstol >();
2070 : : tk::real rel_err = rel_tol+1;
2071 : : tk::real abs_err = abs_tol+1;
2072 : 0 : std::size_t nstiff = m_nstiffeq*ndof;
2073 : :
2074 : : // Initialize Jacobian to be the identity
2075 : : std::vector< std::vector< tk::real > >
2076 [ - - ][ - - ]: 0 : approx_jacob(nstiff, std::vector< tk::real >(nstiff, 0.0));
2077 [ - - ]: 0 : for (std::size_t i=0; i<nstiff; ++i)
2078 : 0 : approx_jacob[i][i] = 1.0e+00;
2079 : :
2080 : : // Save explicit terms to be re-used
2081 [ - - ]: 0 : std::vector< tk::real > expl_terms(nstiff, 0.0);
2082 [ - - ]: 0 : for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2083 [ - - ]: 0 : for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2084 : : {
2085 : 0 : auto stiffmark = m_stiffEqIdx[ieq]*ndof+idof;
2086 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2087 : 0 : expl_terms[ieq*ndof+idof] = m_un(e, stiffrmark)
2088 : 0 : + d->Dt() * ( expl_rkcoef[0][m_stage]
2089 : 0 : * m_rhsprev(e,stiffmark)/m_lhs(e,stiffmark)
2090 : 0 : + expl_rkcoef[1][m_stage]
2091 : 0 : * m_rhs(e,stiffmark)/m_lhs(e,stiffmark)
2092 : 0 : + impl_rkcoef[0][m_stage]
2093 : 0 : * m_stiffrhsprev(e,ieq*ndof+idof)/m_lhs(e,stiffmark) );
2094 : : }
2095 : :
2096 : : // Compute stiff_rhs with initial u
2097 [ - - ][ - - ]: 0 : g_dgpde[d->MeshId()].stiff_rhs( e, myGhosts()->m_geoElem,
2098 [ - - ][ - - ]: 0 : myGhosts()->m_inpoel, myGhosts()->m_coord,
2099 [ - - ]: 0 : m_u, m_p, m_ndof, m_stiffrhs );
2100 : :
2101 : : // Make auxiliary u_old and f_old to store previous values
2102 [ - - ][ - - ]: 0 : std::vector< tk::real > u_old(nstiff, 0.0), f_old(nstiff, 0.0);
[ - - ][ - - ]
2103 : : // Make delta_u and delta_f
2104 [ - - ][ - - ]: 0 : std::vector< tk::real > delta_u(nstiff, 0.0), delta_f(nstiff, 0.0);
[ - - ][ - - ]
2105 : : // Store f
2106 [ - - ][ - - ]: 0 : std::vector< tk::real > f(nstiff, 0.0);
2107 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2108 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2109 : : {
2110 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2111 : 0 : auto stiffmark = m_stiffEqIdx[ieq]*ndof+idof;
2112 : 0 : f[ieq*ndof+idof] = expl_terms[ieq*ndof+idof]
2113 : 0 : + d->Dt() * impl_rkcoef[1][m_stage]
2114 : 0 : * m_stiffrhs(e,ieq*ndof+idof)/m_lhs(e,stiffmark)
2115 : 0 : - m_u(e, stiffrmark);
2116 : : }
2117 : :
2118 : : // Initialize u_old and f_old
2119 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2120 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2121 : : {
2122 : 0 : u_old[ieq*ndof+idof] = m_u(e, m_stiffEqIdx[ieq]*rdof+idof);
2123 : 0 : f_old[ieq*ndof+idof] = f[ieq*ndof+idof];
2124 : : }
2125 : :
2126 : : // Store the norm of f initially, for relative error measure
2127 : : tk::real err0 = 0.0;
2128 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2129 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2130 : 0 : err0 += f[ieq*ndof+idof]*f[ieq*ndof+idof];
2131 : 0 : err0 = std::sqrt(err0);
2132 : :
2133 : : // Iterate for the solution if err0 > 0
2134 [ - - ]: 0 : if (err0 > abs_tol)
2135 [ - - ]: 0 : for (size_t iter=0; iter<max_iter; ++iter)
2136 : : {
2137 : :
2138 : : // Compute new solution
2139 : : tk::real delta;
2140 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2141 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2142 : : {
2143 : : delta = 0.0;
2144 [ - - ]: 0 : for (std::size_t jeq=0; jeq<m_nstiffeq; ++jeq)
2145 [ - - ]: 0 : for (std::size_t jdof=0; jdof<m_numEqDof[jeq]; ++jdof)
2146 : 0 : delta +=
2147 : 0 : approx_jacob[ieq*ndof+idof][jeq*ndof+jdof] * f[jeq*ndof+jdof];
2148 : : // Update u
2149 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2150 : 0 : m_u(e, stiffrmark) -= delta;
2151 : : }
2152 : :
2153 : : // Compute new stiff_rhs
2154 [ - - ][ - - ]: 0 : g_dgpde[d->MeshId()].stiff_rhs( e, myGhosts()->m_geoElem,
2155 [ - - ][ - - ]: 0 : myGhosts()->m_inpoel, myGhosts()->m_coord,
[ - - ]
2156 : : m_u, m_p, m_ndof, m_stiffrhs );
2157 : :
2158 : : // Compute new f(u)
2159 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2160 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2161 : : {
2162 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2163 : 0 : auto stiffmark = m_stiffEqIdx[ieq]*ndof+idof;
2164 : 0 : f[ieq*ndof+idof] = expl_terms[ieq*ndof+idof]
2165 : 0 : + d->Dt() * impl_rkcoef[1][m_stage]
2166 : 0 : * m_stiffrhs(e,ieq*ndof+idof)/m_lhs(e,stiffmark)
2167 : 0 : - m_u(e, stiffrmark);
2168 : : }
2169 : :
2170 : : // Compute delta_u and delta_f
2171 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2172 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2173 : : {
2174 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2175 : 0 : delta_u[ieq*ndof+idof] = m_u(e, stiffrmark) - u_old[ieq*ndof+idof];
2176 : 0 : delta_f[ieq*ndof+idof] = f[ieq*ndof+idof] - f_old[ieq*ndof+idof];
2177 : : }
2178 : :
2179 : : // Update inverse Jacobian approximation
2180 : :
2181 : : // 1. Compute approx_jacob*delta_f and delta_u*jacob_approx
2182 : : tk::real sum1, sum2;
2183 [ - - ][ - - ]: 0 : std::vector< tk::real > auxvec1(nstiff, 0.0), auxvec2(nstiff, 0.0);
[ - - ][ - - ]
2184 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2185 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2186 : : {
2187 : : sum1 = 0.0;
2188 : : sum2 = 0.0;
2189 [ - - ]: 0 : for (std::size_t jeq=0; jeq<m_nstiffeq; ++jeq)
2190 [ - - ]: 0 : for (std::size_t jdof=0; jdof<m_numEqDof[jeq]; ++jdof)
2191 : : {
2192 : 0 : sum1 += approx_jacob[ieq*ndof+idof][jeq*ndof+jdof] *
2193 : 0 : delta_f[jeq*ndof+jdof];
2194 : 0 : sum2 += delta_u[jeq*ndof+jdof] *
2195 : 0 : approx_jacob[jeq*ndof+jdof][ieq*ndof+idof];
2196 : : }
2197 : 0 : auxvec1[ieq*ndof+idof] = sum1;
2198 : 0 : auxvec2[ieq*ndof+idof] = sum2;
2199 : : }
2200 : :
2201 : : // 2. Compute delta_u*approx_jacob*delta_f
2202 : : // and delta_u-approx_jacob*delta_f
2203 : : tk::real denom = 0.0;
2204 [ - - ]: 0 : for (std::size_t jeq=0; jeq<m_nstiffeq; ++jeq)
2205 [ - - ]: 0 : for (std::size_t jdof=0; jdof<m_numEqDof[jeq]; ++jdof)
2206 : : {
2207 : 0 : denom += delta_u[jeq*ndof+jdof]*auxvec1[jeq*ndof+jdof];
2208 : 0 : auxvec1[jeq*ndof+jdof] =
2209 : 0 : delta_u[jeq*ndof+jdof]-auxvec1[jeq*ndof+jdof];
2210 : : }
2211 : :
2212 : : // 3. Divide delta_u+approx_jacob*delta_f
2213 : : // by delta_u*(approx_jacob*delta_f)
2214 [ - - ]: 0 : if (std::abs(denom) < 1.0e-18)
2215 : : {
2216 [ - - ]: 0 : if (denom < 0.0)
2217 : : {
2218 [ - - ]: 0 : for (std::size_t jeq=0; jeq<m_nstiffeq; ++jeq)
2219 [ - - ]: 0 : for (std::size_t jdof=0; jdof<m_numEqDof[jeq]; ++jdof)
2220 : 0 : auxvec1[jeq*ndof+jdof] /= -1.0e-18;
2221 : : }
2222 : : else
2223 : : {
2224 [ - - ]: 0 : for (std::size_t jeq=0; jeq<m_nstiffeq; ++jeq)
2225 [ - - ]: 0 : for (std::size_t jdof=0; jdof<m_numEqDof[jeq]; ++jdof)
2226 : 0 : auxvec1[jeq*ndof+jdof] /= 1.0e-18;
2227 : : }
2228 : : }
2229 : : else
2230 : : {
2231 [ - - ]: 0 : for (std::size_t jeq=0; jeq<m_nstiffeq; ++jeq)
2232 [ - - ]: 0 : for (std::size_t jdof=0; jdof<m_numEqDof[jeq]; ++jdof)
2233 : 0 : auxvec1[jeq*ndof+jdof] /= denom;
2234 : : }
2235 : :
2236 : : // 4. Perform outter product between the two arrays and
2237 : : // add that quantity to the new jacobian approximation
2238 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2239 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2240 [ - - ]: 0 : for (std::size_t jeq=0; jeq<m_nstiffeq; ++jeq)
2241 [ - - ]: 0 : for (std::size_t jdof=0; jdof<m_numEqDof[jeq]; ++jdof)
2242 : 0 : approx_jacob[ieq*ndof+idof][jeq*ndof+jdof] +=
2243 : 0 : auxvec1[ieq*ndof+idof] * auxvec2[jeq*ndof+jdof];
2244 : :
2245 : : // Save solution and f
2246 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2247 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2248 : : {
2249 : 0 : u_old[ieq*ndof+idof] = m_u(e, m_stiffEqIdx[ieq]*rdof+idof);
2250 : 0 : f_old[ieq*ndof+idof] = f[ieq*ndof+idof];
2251 : : }
2252 : :
2253 : : // Compute a measure of error, use norm of f
2254 : : tk::real err = 0.0;
2255 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2256 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2257 : 0 : err += f[ieq*ndof+idof]*f[ieq*ndof+idof];
2258 : 0 : abs_err = std::sqrt(err);
2259 : 0 : rel_err = abs_err/err0;
2260 : :
2261 : : // Check if error condition is met and loop back
2262 [ - - ][ - - ]: 0 : if (rel_err < rel_tol || abs_err < abs_tol)
2263 : : break;
2264 : :
2265 : : // If we did not converge, print a message
2266 [ - - ]: 0 : if (iter == max_iter-1)
2267 : : {
2268 [ - - ]: 0 : printf("\nIMEX-RK: Non-linear solver did not converge in %lu iterations\n", max_iter);
2269 [ - - ]: 0 : printf("Element #%lu\n", e);
2270 [ - - ]: 0 : printf("Relative error: %e\n", rel_err);
2271 [ - - ]: 0 : printf("Absolute error: %e\n\n", abs_err);
2272 : : }
2273 : : }
2274 : : }
2275 : :
2276 : : // Then, integrate explicitly on the remaining equations
2277 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
2278 [ - - ]: 0 : for (std::size_t c=0; c<m_nnonstiffeq; ++c)
2279 : : {
2280 [ - - ]: 0 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
2281 : : {
2282 [ - - ]: 0 : auto rmark = m_nonStiffEqIdx[c]*rdof+k;
2283 [ - - ]: 0 : auto mark = m_nonStiffEqIdx[c]*ndof+k;
2284 : 0 : m_u(e, rmark) = m_un(e, rmark) + d->Dt() * (
2285 : 0 : expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/m_lhs(e, mark)
2286 : 0 : + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/m_lhs(e, mark));
2287 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2288 : 0 : m_u(e, rmark) = 0;
2289 : : }
2290 : : }
2291 : : }
2292 : : else {
2293 : : // For last stage just use all previously computed stages
2294 : 0 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
2295 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
2296 : : {
2297 : : // First integrate explicitly on nonstiff equations
2298 [ - - ]: 0 : for (std::size_t c=0; c<m_nnonstiffeq; ++c)
2299 : : {
2300 [ - - ]: 0 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
2301 : : {
2302 [ - - ]: 0 : auto rmark = m_nonStiffEqIdx[c]*rdof+k;
2303 [ - - ]: 0 : auto mark = m_nonStiffEqIdx[c]*ndof+k;
2304 : 0 : m_u(e, rmark) = m_un(e, rmark) + d->Dt() * (
2305 : 0 : expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/m_lhs(e, mark)
2306 : 0 : + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/m_lhs(e, mark));
2307 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2308 : 0 : m_u(e, rmark) = 0;
2309 : : }
2310 : : }
2311 : : // Then, integrate the imex-equations
2312 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2313 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2314 : : {
2315 [ - - ]: 0 : auto rmark = m_stiffEqIdx[ieq]*rdof+idof;
2316 [ - - ]: 0 : auto mark = m_stiffEqIdx[ieq]*ndof+idof;
2317 : 0 : m_u(e, rmark) = m_un(e, rmark)
2318 : 0 : + d->Dt() * (expl_rkcoef[0][m_stage]
2319 : 0 : * m_rhsprev(e,mark)/m_lhs(e,mark)
2320 : 0 : + expl_rkcoef[1][m_stage]
2321 : 0 : * m_rhs(e,mark)/m_lhs(e,mark)
2322 : 0 : + impl_rkcoef[0][m_stage]
2323 : 0 : * m_stiffrhsprev(e,ieq*ndof+idof)/m_lhs(e,mark)
2324 : 0 : + impl_rkcoef[1][m_stage]
2325 : 0 : * m_stiffrhs(e,ieq*ndof+idof)/m_lhs(e,mark) );
2326 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2327 : 0 : m_u(e, rmark) = 0;
2328 : : }
2329 : : }
2330 : : }
2331 : 0 : }
2332 : :
2333 : : #include "NoWarning/dg.def.h"
|
