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 : : #include "Integrate/Quadrature.hpp"
37 : :
38 : : #include <fstream>
39 : :
40 : : // ignore old-style-casts required for lapack/blas calls
41 : : #if defined(__clang__)
42 : : #pragma clang diagnostic ignored "-Wold-style-cast"
43 : : #endif
44 : :
45 : : // Lapacke forward declarations
46 : : extern "C" {
47 : :
48 : : using lapack_int = long;
49 : :
50 : : #define LAPACK_ROW_MAJOR 101
51 : :
52 : : extern lapack_int LAPACKE_dgesv( int, lapack_int, lapack_int, double*,
53 : : lapack_int, lapack_int*, double*, lapack_int );
54 : :
55 : : }
56 : :
57 : : namespace inciter {
58 : :
59 : : extern ctr::InputDeck g_inputdeck;
60 : : extern std::vector< DGPDE > g_dgpde;
61 : :
62 : : //! Runge-Kutta coefficients
63 : : static const std::array< std::array< tk::real, 3 >, 2 >
64 : : rkcoef{{ {{ 0.0, 3.0/4.0, 1.0/3.0 }}, {{ 1.0, 1.0/4.0, 2.0/3.0 }} }};
65 : :
66 : : //! Implicit-Explicit Runge-Kutta Coefficients
67 : : static const tk::real rk_gamma = (2.0-std::sqrt(2.0))/2.0;
68 : : static const tk::real rk_delta = -2.0*std::sqrt(2.0)/3.0;
69 : : static const tk::real c2 =
70 : : (27.0 + std::pow(2187.0-1458.0*std::sqrt(2.0),1.0/3.0)
71 : : + 9.0*std::pow(3.0+2.0*std::sqrt(2.0),1.0/3.0))/54.0;
72 : : static const tk::real c3 = c2/(6.0*std::pow(c2,2.0)-3.0*c2+1.0);
73 : : static const tk::real b2 = (3.0*c2-1.0)/(6.0*std::pow(c2,2.0));
74 : : static const tk::real b3 =
75 : : (6.0*std::pow(c2,2.0)-3.0*c2+1.0)/(6.0*std::pow(c2,2.0));
76 : : static const tk::real a22_impl = c2;
77 : : static const tk::real a21_expl = c2;
78 : : static const tk::real a32_expl = c3;
79 : : static const tk::real a33_impl =
80 : : (1.0/6.0-b2*std::pow(c2,2.0)-b3*c2*c3)/(b3*(c3-c2));
81 : : static const tk::real a32_impl = a33_impl-c3;
82 : : static const std::array< std::array< tk::real, 3 >, 2 >
83 : : expl_rkcoef{{ {{ 0.0, 0.0, b2 }},
84 : : {{ a21_expl, a32_expl, b3 }} }};
85 : : static const std::array< std::array< tk::real, 3 >, 2>
86 : : impl_rkcoef{{ {{ 0.0, a32_impl, b2 }},
87 : : {{ a22_impl, a33_impl, b3}} }};
88 : : static const std::array< tk::real, 10 > mass_dubiner( tk::massMatrixDubiner() );
89 : :
90 : : } // inciter::
91 : :
92 : : extern tk::CProxy_ChareStateCollector stateProxy;
93 : :
94 : : using inciter::DG;
95 : :
96 : 665 : DG::DG( const CProxy_Discretization& disc,
97 : : const CProxy_Ghosts& ghostsproxy,
98 : : const std::map< int, std::vector< std::size_t > >& bface,
99 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
100 : 665 : const std::vector< std::size_t >& triinpoel ) :
101 : : m_disc( disc ),
102 : : m_ghosts( ghostsproxy ),
103 : : m_ndof_NodalExtrm( 3 ), // for the first order derivatives in 3 directions
104 : : m_nsol( 0 ),
105 : : m_ninitsol( 0 ),
106 : : m_nlim( 0 ),
107 : : m_nnod( 0 ),
108 : : m_nrefine( 0 ),
109 : : m_nsmooth( 0 ),
110 : : m_nreco( 0 ),
111 : : m_nnodalExtrema( 0 ),
112 [ + - ]: 665 : m_nstiffeq( g_dgpde[Disc()->MeshId()].nstiffeq() ),
113 [ + - ]: 665 : m_nnonstiffeq( g_dgpde[Disc()->MeshId()].nnonstiffeq() ),
114 [ + - ]: 665 : m_u( Disc()->Inpoel().size()/4,
115 : 665 : g_inputdeck.get< tag::rdof >()*
116 : 665 : g_inputdeck.get< tag::ncomp >() ),
117 : : m_un( m_u.nunk(), m_u.nprop() ),
118 : 665 : m_p( m_u.nunk(), g_inputdeck.get< tag::rdof >()*
119 [ + - ][ + - ]: 665 : g_dgpde[Disc()->MeshId()].nprim() ),
120 : : m_rhs( m_u.nunk(),
121 : 665 : g_inputdeck.get< tag::ndof >()*
122 : 665 : g_inputdeck.get< tag::ncomp >() ),
123 : : m_rhsprev( m_u.nunk(), m_rhs.nprop() ),
124 : 665 : m_stiffrhs( m_u.nunk(), g_inputdeck.get< tag::ndof >()*
125 [ + - ][ + - ]: 665 : g_dgpde[Disc()->MeshId()].nstiffeq() ),
126 : 665 : m_stiffrhsprev( m_u.nunk(), g_inputdeck.get< tag::ndof >()*
127 [ + - ][ + - ]: 665 : g_dgpde[Disc()->MeshId()].nstiffeq() ),
128 [ + - ]: 665 : m_stiffEqIdx( g_dgpde[Disc()->MeshId()].nstiffeq() ),
129 [ + - ]: 665 : m_nonStiffEqIdx( g_dgpde[Disc()->MeshId()].nnonstiffeq() ),
130 : : m_mtInv(
131 : : tk::invMassMatTaylorRefEl(g_inputdeck.get< tag::rdof >()) ),
132 : : m_uNodalExtrm(),
133 : : m_pNodalExtrm(),
134 : : m_uNodalExtrmc(),
135 : : m_pNodalExtrmc(),
136 [ + - ]: 665 : m_npoin( Disc()->Coord()[0].size() ),
137 : : m_diag(),
138 : : m_nstage( 3 ),
139 : : m_stage( 0 ),
140 : : m_ndof(),
141 : : m_interface(),
142 : : m_numEqDof(),
143 : : m_uc(),
144 : : m_pc(),
145 : : m_ndofc(),
146 : : m_interfacec(),
147 : : m_initial( 1 ),
148 : : m_uElemfields( m_u.nunk(),
149 : : g_inputdeck.get< tag::ncomp >() ),
150 : : m_pElemfields( m_u.nunk(),
151 : 665 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
152 : : m_uNodefields( m_npoin,
153 : : g_inputdeck.get< tag::ncomp >() ),
154 : : m_pNodefields( m_npoin,
155 : 665 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
156 : : m_uNodefieldsc(),
157 : : m_pNodefieldsc(),
158 : : m_outmesh(),
159 : : m_boxelems(),
160 [ + - ][ + - ]: 7315 : m_shockmarker(m_u.nunk(), 1)
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
161 : : // *****************************************************************************
162 : : // Constructor
163 : : //! \param[in] disc Discretization proxy
164 : : //! \param[in] bface Boundary-faces mapped to side set ids
165 : : //! \param[in] triinpoel Boundary-face connectivity
166 : : // *****************************************************************************
167 : : {
168 [ + + ]: 665 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
169 [ + + ]: 627 : g_inputdeck.get< tag::cmd, tag::quiescence >())
170 [ + - ][ + - ]: 880 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
171 : : "DG" );
172 : :
173 : : // assign number of dofs for each equation in all pde systems
174 [ + - ][ + - ]: 665 : g_dgpde[Disc()->MeshId()].numEquationDofs(m_numEqDof);
175 : :
176 : : // Allocate storage for the vector of nodal extrema
177 [ + - ][ + - ]: 665 : m_uNodalExtrm.resize( Disc()->Bid().size(),
178 : 0 : std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
179 [ + - ]: 665 : g_inputdeck.get< tag::ncomp >() ) );
180 [ + - ][ + - ]: 665 : m_pNodalExtrm.resize( Disc()->Bid().size(),
181 : 0 : std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
182 [ + - ]: 665 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ) );
183 : :
184 : : // Initialization for the buffer vector of nodal extrema
185 [ + - ]: 665 : resizeNodalExtremac();
186 : :
187 : 665 : usesAtSync = true; // enable migration at AtSync
188 : :
189 : 665 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
190 : :
191 : : // Enable SDAG wait for initially building the solution vector and limiting
192 [ + - ]: 665 : if (m_initial) {
193 [ + - ][ + - ]: 665 : thisProxy[ thisIndex ].wait4sol();
194 [ + + ][ + - ]: 799 : if (pref) thisProxy[ thisIndex ].wait4refine();
[ + - ]
195 [ + - ][ + - ]: 665 : thisProxy[ thisIndex ].wait4smooth();
196 [ + - ][ + - ]: 665 : thisProxy[ thisIndex ].wait4lim();
197 [ + - ][ + - ]: 665 : thisProxy[ thisIndex ].wait4nod();
198 [ + - ][ + - ]: 665 : thisProxy[ thisIndex ].wait4reco();
199 [ + - ][ + - ]: 1330 : thisProxy[ thisIndex ].wait4nodalExtrema();
200 : : }
201 : :
202 [ + - ][ + - ]: 1330 : m_ghosts[thisIndex].insert(m_disc, bface, triinpoel, m_u.nunk(),
203 [ + - ][ + - ]: 1995 : CkCallback(CkIndex_DG::resizeSolVectors(), thisProxy[thisIndex]));
[ - + ][ - + ]
[ - - ][ - - ]
204 : :
205 : : // insert array-element into the implicit solver chare array
206 [ - + ]: 665 : if (g_inputdeck.get< tag::implicit_timestepping >()) {
207 : : // Single-stage BDF1 for implicit solver
208 : 0 : m_nstage = 1;
209 : :
210 [ - - ]: 0 : const auto& inpoel = myGhosts()->m_inpoel;
211 : : // TODO: linear solver:
212 : : // modify CSR to handle element-based structures (or create new one)
213 [ - - ][ - - ]: 0 : tk::CSR A(m_rhs.nprop(), tk::genPsup(inpoel,4,tk::genEsup(inpoel,4)));
[ - - ]
214 [ - - ]: 0 : std::vector< tk::real > x(m_u.nunk()*m_rhs.nprop(), 0.0),
215 [ - - ][ - - ]: 0 : b(m_u.nunk()*m_rhs.nprop(), 0.0);
216 : :
217 [ - - ][ - - ]: 0 : Disc()->ImplicitSolver()[ thisIndex ].insert(std::move(A), std::move(x),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
218 : : std::move(b), Disc()->Gid(), Disc()->Lid(), Disc()->NodeCommMap());
219 : : }
220 : :
221 : : // global-sync to call doneInserting on m_ghosts
222 [ + - ]: 665 : auto meshid = Disc()->MeshId();
223 [ + - ]: 665 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
224 [ + - ][ - - ]: 665 : CkCallback(CkReductionTarget(Transporter,doneInsertingGhosts),
225 [ + - ]: 665 : Disc()->Tr()) );
226 : 665 : }
227 : :
228 : : void
229 : 523 : DG::registerReducers()
230 : : // *****************************************************************************
231 : : // Configure Charm++ reduction types
232 : : //! \details Since this is a [initnode] routine, the runtime system executes the
233 : : //! routine exactly once on every logical node early on in the Charm++ init
234 : : //! sequence. Must be static as it is called without an object. See also:
235 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
236 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
237 : : // *****************************************************************************
238 : : {
239 : 523 : ElemDiagnostics::registerReducers();
240 : 523 : }
241 : :
242 : : void
243 : 11805 : DG::ResumeFromSync()
244 : : // *****************************************************************************
245 : : // Return from migration
246 : : //! \details This is called when load balancing (LB) completes. The presence of
247 : : //! this function does not affect whether or not we block on LB.
248 : : // *****************************************************************************
249 : : {
250 [ - + ][ - - ]: 11805 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
251 : :
252 [ + - ]: 11805 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
253 : 11805 : }
254 : :
255 : : void
256 : 665 : DG::resizeSolVectors()
257 : : // *****************************************************************************
258 : : // Resize solution vectors after extension due to Ghosts and continue with setup
259 : : // *****************************************************************************
260 : : {
261 : : // Resize solution vectors, lhs and rhs by the number of ghost tets
262 : 665 : m_u.resize( myGhosts()->m_nunk );
263 : 665 : m_un.resize( myGhosts()->m_nunk );
264 : 665 : m_p.resize( myGhosts()->m_nunk );
265 : 665 : m_rhs.resize( myGhosts()->m_nunk );
266 : 665 : m_rhsprev.resize( myGhosts()->m_nunk );
267 : 665 : m_stiffrhs.resize( myGhosts()->m_nunk );
268 : 665 : m_stiffrhsprev.resize( myGhosts()->m_nunk );
269 : :
270 : : // Size communication buffer for solution and number of degrees of freedom
271 [ + + ]: 2660 : for (auto& n : m_ndofc) n.resize( myGhosts()->m_bid.size() );
272 [ + + ]: 2660 : for (auto& u : m_uc) u.resize( myGhosts()->m_bid.size() );
273 [ + + ]: 2660 : for (auto& p : m_pc) p.resize( myGhosts()->m_bid.size() );
274 [ + + ]: 1330 : for (auto& i : m_interfacec) i.resize( myGhosts()->m_bid.size() );
275 : :
276 : : // Initialize number of degrees of freedom in mesh elements
277 : 665 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
278 [ + + ]: 665 : if( pref )
279 : : {
280 : 134 : const auto ndofmax = g_inputdeck.get< tag::pref, tag::ndofmax >();
281 : 134 : m_ndof.resize( myGhosts()->m_nunk, ndofmax );
282 : : }
283 : : else
284 : : {
285 : 531 : const auto ndof = g_inputdeck.get< tag::ndof >();
286 : 531 : m_ndof.resize( myGhosts()->m_nunk, ndof );
287 : : }
288 : 665 : m_interface.resize( myGhosts()->m_nunk, 0 );
289 : :
290 : : // Ensure that we also have all the geometry and connectivity data
291 : : // (including those of ghosts)
292 : : Assert( myGhosts()->m_geoElem.nunk() == m_u.nunk(),
293 : : "GeoElem unknowns size mismatch" );
294 : :
295 : : // Signal the runtime system that all workers have received their adjacency
296 : 665 : std::vector< std::size_t > meshdata{ myGhosts()->m_initial, Disc()->MeshId() };
297 : 665 : contribute( meshdata, CkReduction::sum_ulong,
298 [ + - ][ + - ]: 1995 : CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
[ + - ][ - - ]
299 : 665 : }
300 : :
301 : : void
302 : 665 : DG::setup()
303 : : // *****************************************************************************
304 : : // Set initial conditions, generate lhs, output mesh
305 : : // *****************************************************************************
306 : : {
307 [ + + ]: 665 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
308 [ + + ]: 627 : g_inputdeck.get< tag::cmd, tag::quiescence >())
309 [ + - ][ + - ]: 880 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
310 : : "setup" );
311 : :
312 : 665 : auto d = Disc();
313 : :
314 : : // Compute left-hand side of discrete PDEs
315 : : lhs();
316 : :
317 : : // Determine elements inside user-defined IC box
318 : 665 : g_dgpde[d->MeshId()].IcBoxElems( myGhosts()->m_geoElem,
319 : 665 : myGhosts()->m_fd.Esuel().size()/4, m_boxelems );
320 : :
321 : : // Compute volume of user-defined box IC
322 [ + - ][ - + ]: 665 : d->boxvol( {}, {}, 0 ); // punt for now
323 : :
324 : : // Query time history field output labels from all PDEs integrated
325 : : const auto& hist_points = g_inputdeck.get< tag::history_output, tag::point >();
326 [ - + ]: 665 : if (!hist_points.empty()) {
327 : 0 : std::vector< std::string > histnames;
328 [ - - ]: 0 : auto n = g_dgpde[d->MeshId()].histNames();
329 [ - - ]: 0 : histnames.insert( end(histnames), begin(n), end(n) );
330 [ - - ]: 0 : d->histheader( std::move(histnames) );
331 : : }
332 : :
333 : : // If working with IMEX-RK, Store stiff equations into m_stiffEqIdx
334 [ - + ]: 665 : if (g_inputdeck.get< tag::imex_runge_kutta >())
335 : : {
336 : 0 : g_dgpde[Disc()->MeshId()].setStiffEqIdx(m_stiffEqIdx);
337 : 0 : g_dgpde[Disc()->MeshId()].setNonStiffEqIdx(m_nonStiffEqIdx);
338 : : }
339 : 665 : }
340 : :
341 : : void
342 : 665 : DG::box( tk::real v, const std::vector< tk::real >& )
343 : : // *****************************************************************************
344 : : // Receive total box IC volume and set conditions in box
345 : : //! \param[in] v Total volume within user-specified box
346 : : // *****************************************************************************
347 : : {
348 : 665 : auto d = Disc();
349 : :
350 : : // Store user-defined box IC volume
351 : 665 : d->Boxvol() = v;
352 : :
353 : : // Set initial conditions for all PDEs
354 : 1330 : g_dgpde[d->MeshId()].initialize( myGhosts()->m_geoElem, myGhosts()->m_inpoel,
355 : 665 : myGhosts()->m_coord, m_boxelems, d->ElemBlockId(), m_u, d->T(),
356 : 665 : myGhosts()->m_fd.Esuel().size()/4 );
357 : 665 : g_dgpde[d->MeshId()].updatePrimitives( m_u, myGhosts()->m_geoElem, m_p,
358 : 665 : myGhosts()->m_fd.Esuel().size()/4, m_ndof );
359 : :
360 : : m_un = m_u;
361 : :
362 : : // Output initial conditions to file (regardless of whether it was requested)
363 [ + - ][ + - ]: 1995 : startFieldOutput( CkCallback(CkIndex_DG::start(), thisProxy[thisIndex]) );
[ - + ][ - - ]
364 : 665 : }
365 : :
366 : : void
367 : 665 : DG::start()
368 : : // *****************************************************************************
369 : : // Start time stepping
370 : : // *****************************************************************************
371 : : {
372 : : // Free memory storing output mesh
373 : 665 : m_outmesh.destroy();
374 : :
375 : : // Start timer measuring time stepping wall clock time
376 : 665 : Disc()->Timer().zero();
377 : : // Zero grind-timer
378 : 665 : Disc()->grindZero();
379 : : // Start time stepping by computing the size of the next time step)
380 : 665 : next();
381 : 665 : }
382 : :
383 : : void
384 : 13135 : DG::startFieldOutput( CkCallback c )
385 : : // *****************************************************************************
386 : : // Start preparing fields for output to file
387 : : //! \param[in] c Function to continue with after the write
388 : : // *****************************************************************************
389 : : {
390 : : // No field output in benchmark mode or if field output frequency not hit
391 [ + + ][ + + ]: 13135 : if (g_inputdeck.get< tag::cmd, tag::benchmark >() || !fieldOutput()) {
392 : :
393 : 11857 : c.send();
394 : :
395 : : } else {
396 : :
397 : : // Optionally refine mesh for field output
398 : 1278 : auto d = Disc();
399 : :
400 [ + + ]: 1278 : if (refinedOutput()) {
401 : :
402 : 33 : const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
403 [ + - ][ + - ]: 66 : d->Ref()->outref( myGhosts()->m_fd.Bface(), {}, tr, c );
[ + - ][ + - ]
[ - - ]
404 : :
405 : : } else {
406 : :
407 : : // cut off ghosts from mesh connectivity and coordinates
408 : 1245 : const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
409 [ + - ][ + - ]: 2490 : extractFieldOutput( {}, d->Chunk(), d->Coord(), {}, {},
[ + + ][ - - ]
410 [ + - ]: 1245 : d->NodeCommMap(), myGhosts()->m_fd.Bface(), {}, tr, c );
411 : :
412 : : }
413 : :
414 : : }
415 : 13135 : }
416 : :
417 : : void
418 : 37410 : DG::next()
419 : : // *****************************************************************************
420 : : // Advance equations to next time step
421 : : // *****************************************************************************
422 : : {
423 : 37410 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
424 : :
425 : 37410 : auto d = Disc();
426 : :
427 [ + + ][ + + ]: 37410 : if (pref && m_stage == 0 && d->T() > 0)
[ + + ]
428 : 3272 : g_dgpde[d->MeshId()].eval_ndof( myGhosts()->m_nunk, myGhosts()->m_coord,
429 : 1636 : myGhosts()->m_inpoel,
430 : 1636 : myGhosts()->m_fd, m_u, m_p,
431 : : g_inputdeck.get< tag::pref, tag::indicator >(),
432 : : g_inputdeck.get< tag::ndof >(),
433 : : g_inputdeck.get< tag::pref, tag::ndofmax >(),
434 : : g_inputdeck.get< tag::pref, tag::tolref >(),
435 : 1636 : m_ndof );
436 : :
437 : : // communicate solution ghost data (if any)
438 [ + + ]: 37410 : if (myGhosts()->m_sendGhost.empty())
439 : 3195 : comsol_complete();
440 : : else
441 [ + + ]: 451500 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
442 [ + - ]: 383070 : std::vector< std::size_t > tetid( ghostdata.size() );
443 [ + - ][ + - ]: 766140 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
444 [ + - ][ + - ]: 766140 : prim( ghostdata.size() );
445 [ + - ][ + - ]: 383070 : std::vector< std::size_t > interface( ghostdata.size() );
446 [ + - ][ - - ]: 383070 : std::vector< std::size_t > ndof( ghostdata.size() );
447 : : std::size_t j = 0;
448 [ + + ]: 6733830 : for(const auto& i : ghostdata) {
449 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
450 : : "Sending solution ghost data" );
451 [ + - ]: 6350760 : tetid[j] = i;
452 [ + - ][ - + ]: 6350760 : u[j] = m_u[i];
453 [ + - ][ - + ]: 6350760 : prim[j] = m_p[i];
454 [ + + ][ + + ]: 6350760 : if (pref && m_stage == 0) {
455 : 395110 : ndof[j] = m_ndof[i];
456 : 395110 : interface[j] = m_interface[i];
457 : : }
458 : 6350760 : ++j;
459 : : }
460 [ + - ][ + - ]: 766140 : thisProxy[ cid ].comsol( thisIndex, m_stage, tetid, u, prim, interface, ndof );
[ + - ][ - - ]
461 : : }
462 : :
463 : 37410 : ownsol_complete();
464 : 37410 : }
465 : :
466 : : void
467 : 383070 : DG::comsol( int fromch,
468 : : std::size_t fromstage,
469 : : const std::vector< std::size_t >& tetid,
470 : : const std::vector< std::vector< tk::real > >& u,
471 : : const std::vector< std::vector< tk::real > >& prim,
472 : : const std::vector< std::size_t >& interface,
473 : : const std::vector< std::size_t >& ndof )
474 : : // *****************************************************************************
475 : : // Receive chare-boundary solution ghost data from neighboring chares
476 : : //! \param[in] fromch Sender chare id
477 : : //! \param[in] fromstage Sender chare time step stage
478 : : //! \param[in] tetid Ghost tet ids we receive solution data for
479 : : //! \param[in] u Solution ghost data
480 : : //! \param[in] prim Primitive variables in ghost cells
481 : : //! \param[in] interface Interface marker in ghost cells
482 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
483 : : //! \details This function receives contributions to the unlimited solution
484 : : //! from fellow chares.
485 : : // *****************************************************************************
486 : : {
487 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comsol()" );
488 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comsol()" );
489 : :
490 : 383070 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
491 : :
492 : 383070 : if (pref && fromstage == 0) {
493 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsol()" );
494 : : Assert( interface.size() == tetid.size(), "Size mismatch in DG::comsol()" );
495 : : }
496 : :
497 : : // Find local-to-ghost tet id map for sender chare
498 : 383070 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
499 : :
500 [ + + ]: 6733830 : for (std::size_t i=0; i<tetid.size(); ++i) {
501 : 6350760 : auto j = tk::cref_find( n, tetid[i] );
502 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
503 : : "Receiving solution non-ghost data" );
504 : 6350760 : auto b = tk::cref_find( myGhosts()->m_bid, j );
505 : : Assert( b < m_uc[0].size(), "Indexing out of bounds" );
506 : : Assert( b < m_pc[0].size(), "Indexing out of bounds" );
507 : 6350760 : m_uc[0][b] = u[i];
508 : 6350760 : m_pc[0][b] = prim[i];
509 [ + + ]: 6350760 : if (pref && fromstage == 0) {
510 : : Assert( b < m_ndofc[0].size(), "Indexing out of bounds" );
511 : 395110 : m_ndofc[0][b] = ndof[i];
512 : : Assert( b < m_interfacec[0].size(), "Indexing out of bounds" );
513 : 395110 : m_interfacec[0][b] = interface[i];
514 : : }
515 : : }
516 : :
517 : : // if we have received all solution ghost contributions from neighboring
518 : : // chares (chares we communicate along chare-boundary faces with), and
519 : : // contributed our solution to these neighbors, proceed to reconstructions
520 [ + + ]: 383070 : if (++m_nsol == myGhosts()->m_sendGhost.size()) {
521 : 34215 : m_nsol = 0;
522 : 34215 : comsol_complete();
523 : : }
524 : 383070 : }
525 : :
526 : : void
527 : 1278 : DG::extractFieldOutput(
528 : : const std::vector< std::size_t >& /*ginpoel*/,
529 : : const tk::UnsMesh::Chunk& chunk,
530 : : const tk::UnsMesh::Coords& coord,
531 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
532 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
533 : : const tk::NodeCommMap& nodeCommMap,
534 : : const std::map< int, std::vector< std::size_t > >& bface,
535 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
536 : : const std::vector< std::size_t >& triinpoel,
537 : : CkCallback c )
538 : : // *****************************************************************************
539 : : // Extract field output going to file
540 : : //! \param[in] chunk Field-output mesh chunk (connectivity and global<->local
541 : : //! id maps)
542 : : //! \param[in] coord Field-output mesh node coordinates
543 : : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
544 : : //! \param[in] nodeCommMap Field-output mesh node communication map
545 : : //! \param[in] bface Field-output meshndary-faces mapped to side set ids
546 : : //! \param[in] triinpoel Field-output mesh boundary-face connectivity
547 : : //! \param[in] c Function to continue with after the write
548 : : // *****************************************************************************
549 : : {
550 : : m_outmesh.chunk = chunk;
551 : : m_outmesh.coord = coord;
552 : 1278 : m_outmesh.triinpoel = triinpoel;
553 : : m_outmesh.bface = bface;
554 : : m_outmesh.nodeCommMap = nodeCommMap;
555 : :
556 : : const auto& inpoel = std::get< 0 >( chunk );
557 : :
558 : : // Evaluate element solution on incoming mesh
559 : 1278 : evalSolution( *Disc(), inpoel, coord, addedTets, m_ndof, m_u, m_p,
560 : 1278 : m_uElemfields, m_pElemfields, m_uNodefields, m_pNodefields );
561 : :
562 : : // Send node fields contributions to neighbor chares
563 [ + + ]: 1278 : if (nodeCommMap.empty())
564 : 150 : comnodeout_complete();
565 : : else {
566 : : const auto& lid = std::get< 2 >( chunk );
567 : 2256 : auto esup = tk::genEsup( inpoel, 4 );
568 [ + + ]: 11400 : for(const auto& [ch,nodes] : nodeCommMap) {
569 : : // Pack node field data in chare boundary nodes
570 : : std::vector< std::vector< tk::real > >
571 [ + - ][ + - ]: 30816 : lu( m_uNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
[ + - ]
572 : : std::vector< std::vector< tk::real > >
573 [ + - ][ + - ]: 30816 : lp( m_pNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
574 [ + + ]: 68524 : for (std::size_t f=0; f<m_uNodefields.nprop(); ++f) {
575 : : std::size_t j = 0;
576 [ + + ]: 594448 : for (auto g : nodes)
577 : 536196 : lu[f][j++] = m_uNodefields(tk::cref_find(lid,g),f);
578 : : }
579 [ + + ]: 20900 : for (std::size_t f=0; f<m_pNodefields.nprop(); ++f) {
580 : : std::size_t j = 0;
581 [ + + ]: 120108 : for (auto g : nodes)
582 : 109480 : lp[f][j++] = m_pNodefields(tk::cref_find(lid,g),f);
583 : : }
584 : : // Pack (partial) number of elements surrounding chare boundary nodes
585 [ + - ]: 10272 : std::vector< std::size_t > nesup( nodes.size() );
586 : : std::size_t j = 0;
587 [ + + ]: 105698 : for (auto g : nodes) {
588 : 95426 : auto i = tk::cref_find( lid, g );
589 : 95426 : nesup[j++] = esup.second[i+1] - esup.second[i];
590 : : }
591 [ + - ][ + - ]: 30816 : thisProxy[ch].comnodeout(
[ + - ][ - - ]
592 [ + - ][ - + ]: 20544 : std::vector<std::size_t>(begin(nodes),end(nodes)), nesup, lu, lp );
593 : : }
594 : : }
595 : :
596 [ + - ]: 1278 : ownnod_complete( c, addedTets );
597 : 1278 : }
598 : :
599 : 37410 : void DG::p_refine()
600 : : // *****************************************************************************
601 : : // Add the protective layer for ndof refinement
602 : : // *****************************************************************************
603 : : {
604 : 37410 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
605 : :
606 : : // Combine own and communicated contributions of unreconstructed solution and
607 : : // degrees of freedom in cells (if p-adaptive)
608 [ + + ]: 6425580 : for (const auto& b : myGhosts()->m_bid) {
609 : : Assert( m_uc[0][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
610 : : Assert( m_pc[0][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
611 [ + + ]: 167849685 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
612 : 161498925 : m_u(b.first,c) = m_uc[0][b.second][c];
613 : : }
614 [ + + ]: 33384885 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
615 : 27034125 : m_p(b.first,c) = m_pc[0][b.second][c];
616 : : }
617 [ + + ][ + + ]: 6350760 : if (pref && m_stage == 0) {
618 : 395110 : m_ndof[ b.first ] = m_ndofc[0][ b.second ];
619 : 395110 : m_interface[ b.first ] = m_interfacec[0][ b.second ];
620 : : }
621 : : }
622 : :
623 [ + + ][ + + ]: 37410 : if (pref && m_stage==0) refine_ndof();
624 : :
625 [ + + ]: 37410 : if (!pref) {
626 : : // if p-refinement is not configured, proceed directly to reconstructions
627 : 32100 : reco();
628 : : }
629 : : else {
630 : : // if p-refinement is configured, do refine-smoothing before reconstruction
631 [ + + ]: 5310 : if (myGhosts()->m_sendGhost.empty())
632 : 150 : comrefine_complete();
633 : : else
634 [ + + ]: 50160 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
635 [ + - ]: 39840 : std::vector< std::size_t > tetid( ghostdata.size() );
636 [ + - ][ + - ]: 79680 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
637 [ + - ][ + - ]: 79680 : prim( ghostdata.size() );
638 [ + - ]: 39840 : std::vector< std::size_t > ndof( ghostdata.size() );
639 : : std::size_t j = 0;
640 [ + + ]: 1225170 : for(const auto& i : ghostdata) {
641 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending refined ndof "
642 : : "data" );
643 [ + + ]: 1185330 : tetid[j] = i;
644 [ + + ]: 1185330 : if (pref && m_stage == 0) ndof[j] = m_ndof[i];
645 : 1185330 : ++j;
646 : : }
647 [ + - ][ + - ]: 79680 : thisProxy[ cid ].comrefine( thisIndex, tetid, ndof );
[ + - ][ - - ]
648 : : }
649 : :
650 : 5310 : ownrefine_complete();
651 : : }
652 : 37410 : }
653 : :
654 : : void
655 : 39840 : DG::comrefine( int fromch,
656 : : const std::vector< std::size_t >& tetid,
657 : : const std::vector< std::size_t >& ndof )
658 : : // *****************************************************************************
659 : : // Receive chare-boundary ghost data from neighboring chares
660 : : //! \param[in] fromch Sender chare id
661 : : //! \param[in] tetid Ghost tet ids we receive solution data for
662 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
663 : : //! \details This function receives contributions to the refined ndof data
664 : : //! from fellow chares.
665 : : // *****************************************************************************
666 : : {
667 : 39840 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
668 : :
669 : : if (pref && m_stage == 0)
670 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comrefine()" );
671 : :
672 : : // Find local-to-ghost tet id map for sender chare
673 : 39840 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
674 : :
675 [ + + ]: 1225170 : for (std::size_t i=0; i<tetid.size(); ++i) {
676 : 1185330 : auto j = tk::cref_find( n, tetid[i] );
677 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
678 : : "Receiving solution non-ghost data" );
679 : 1185330 : auto b = tk::cref_find( myGhosts()->m_bid, j );
680 [ + - ][ + + ]: 1185330 : if (pref && m_stage == 0) {
681 : : Assert( b < m_ndofc[1].size(), "Indexing out of bounds" );
682 : 395110 : m_ndofc[1][b] = ndof[i];
683 : : }
684 : : }
685 : :
686 : : // if we have received all solution ghost contributions from neighboring
687 : : // chares (chares we communicate along chare-boundary faces with), and
688 : : // contributed our solution to these neighbors, proceed to limiting
689 [ + + ]: 39840 : if (++m_nrefine == myGhosts()->m_sendGhost.size()) {
690 : 5160 : m_nrefine = 0;
691 : 5160 : comrefine_complete();
692 : : }
693 : 39840 : }
694 : :
695 : : void
696 : 5310 : DG::smooth()
697 : : // *****************************************************************************
698 : : // Smooth the refined ndof distribution
699 : : // *****************************************************************************
700 : : {
701 : 5310 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
702 : :
703 [ + + ]: 1195950 : for (const auto& b : myGhosts()->m_bid) {
704 [ + - ][ + + ]: 1185330 : if (pref && m_stage == 0)
705 : 395110 : m_ndof[ b.first ] = m_ndofc[1][ b.second ];
706 : : }
707 : :
708 [ + - ][ + + ]: 5310 : if (pref && m_stage==0) smooth_ndof();
709 : :
710 [ + + ]: 5310 : if (myGhosts()->m_sendGhost.empty())
711 : 150 : comsmooth_complete();
712 : : else
713 [ + + ]: 50160 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
714 : 39840 : std::vector< std::size_t > tetid( ghostdata.size() );
715 : : std::vector< std::size_t > ndof;
716 : : std::size_t j = 0;
717 [ + + ]: 1225170 : for(const auto& i : ghostdata) {
718 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending ndof data" );
719 [ + - ]: 1185330 : tetid[j] = i;
720 [ + - ][ + + ]: 1185330 : if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
[ + - ]
721 : 1185330 : ++j;
722 : : }
723 [ + - ][ + - ]: 79680 : thisProxy[ cid ].comsmooth( thisIndex, tetid, ndof );
[ + + ][ - - ]
724 : : }
725 : :
726 : 5310 : ownsmooth_complete();
727 : 5310 : }
728 : :
729 : : void
730 : 39840 : DG::comsmooth( int fromch,
731 : : const std::vector< std::size_t >& tetid,
732 : : const std::vector< std::size_t >& ndof )
733 : : // *****************************************************************************
734 : : // Receive chare-boundary ghost data from neighboring chares
735 : : //! \param[in] fromch Sender chare id
736 : : //! \param[in] tetid Ghost tet ids we receive solution data for
737 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
738 : : //! \details This function receives contributions to the smoothed ndof data
739 : : //! from fellow chares.
740 : : // *****************************************************************************
741 : : {
742 : 39840 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
743 : :
744 : : if (pref && m_stage == 0)
745 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsmooth()" );
746 : :
747 : 39840 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
748 : :
749 [ + + ]: 1225170 : for (std::size_t i=0; i<tetid.size(); ++i) {
750 : 1185330 : auto j = tk::cref_find( n, tetid[i] );
751 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4, "Receiving ndof data" );
752 : 1185330 : auto b = tk::cref_find( myGhosts()->m_bid, j );
753 [ + - ][ + + ]: 1185330 : if (pref && m_stage == 0) {
754 : : Assert( b < m_ndofc[2].size(), "Indexing out of bounds" );
755 : 395110 : m_ndofc[2][b] = ndof[i];
756 : : }
757 : : }
758 : :
759 [ + + ]: 39840 : if (++m_nsmooth == myGhosts()->m_sendGhost.size()) {
760 : 5160 : m_nsmooth = 0;
761 : 5160 : comsmooth_complete();
762 : : }
763 : 39840 : }
764 : :
765 : : void
766 : 37410 : DG::reco()
767 : : // *****************************************************************************
768 : : // Compute reconstructions
769 : : // *****************************************************************************
770 : : {
771 : 37410 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
772 : 37410 : const auto rdof = g_inputdeck.get< tag::rdof >();
773 : :
774 : : // Combine own and communicated contributions of unreconstructed solution and
775 : : // degrees of freedom in cells (if p-adaptive)
776 [ + + ]: 37410 : if (pref) {
777 [ + + ]: 1195950 : for (const auto& b : myGhosts()->m_bid) {
778 [ + + ]: 1185330 : if (m_stage == 0) {
779 : 395110 : m_ndof[ b.first ] = m_ndofc[2][ b.second ];
780 : : }
781 : : }
782 : : }
783 : :
784 : 37410 : auto d = Disc();
785 [ + + ][ + + ]: 37410 : if (pref && m_stage==0) {
786 : 1770 : g_dgpde[d->MeshId()].resetAdapSol( myGhosts()->m_fd, m_u, m_p, m_ndof );
787 : : }
788 : :
789 [ + + ]: 37410 : if (rdof > 1)
790 : : // Reconstruct second-order solution and primitive quantities
791 : 47970 : g_dgpde[d->MeshId()].reconstruct( d->T(), myGhosts()->m_geoFace,
792 : 23985 : myGhosts()->m_geoElem,
793 : 23985 : myGhosts()->m_fd, myGhosts()->m_esup, myGhosts()->m_inpoel,
794 : 23985 : myGhosts()->m_coord, m_u, m_p, pref, m_ndof );
795 : :
796 : : // Send reconstructed solution to neighboring chares
797 [ + + ]: 37410 : if (myGhosts()->m_sendGhost.empty())
798 : 3195 : comreco_complete();
799 : : else
800 [ + + ]: 451500 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
801 [ + - ]: 383070 : std::vector< std::size_t > tetid( ghostdata.size() );
802 [ + - ][ + - ]: 766140 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
803 [ + - ]: 383070 : prim( ghostdata.size() );
804 : : std::size_t j = 0;
805 [ + + ]: 6733830 : for(const auto& i : ghostdata) {
806 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending reconstructed ghost "
807 : : "data" );
808 [ + - ]: 6350760 : tetid[j] = i;
809 [ + - ][ - + ]: 6350760 : u[j] = m_u[i];
810 [ + - ][ - + ]: 6350760 : prim[j] = m_p[i];
811 : 6350760 : ++j;
812 : : }
813 [ + - ][ + - ]: 766140 : thisProxy[ cid ].comreco( thisIndex, tetid, u, prim );
814 : : }
815 : :
816 : 37410 : ownreco_complete();
817 : 37410 : }
818 : :
819 : : void
820 : 383070 : DG::comreco( int fromch,
821 : : const std::vector< std::size_t >& tetid,
822 : : const std::vector< std::vector< tk::real > >& u,
823 : : const std::vector< std::vector< tk::real > >& prim )
824 : : // *****************************************************************************
825 : : // Receive chare-boundary reconstructed ghost data from neighboring chares
826 : : //! \param[in] fromch Sender chare id
827 : : //! \param[in] tetid Ghost tet ids we receive solution data for
828 : : //! \param[in] u Reconstructed high-order solution
829 : : //! \param[in] prim Limited high-order primitive quantities
830 : : //! \details This function receives contributions to the reconstructed solution
831 : : //! from fellow chares.
832 : : // *****************************************************************************
833 : : {
834 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comreco()" );
835 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comreco()" );
836 : :
837 : : // Find local-to-ghost tet id map for sender chare
838 : 383070 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
839 : :
840 [ + + ]: 6733830 : for (std::size_t i=0; i<tetid.size(); ++i) {
841 : 6350760 : auto j = tk::cref_find( n, tetid[i] );
842 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
843 : : "Receiving solution non-ghost data" );
844 : 6350760 : auto b = tk::cref_find( myGhosts()->m_bid, j );
845 : : Assert( b < m_uc[1].size(), "Indexing out of bounds" );
846 : : Assert( b < m_pc[1].size(), "Indexing out of bounds" );
847 : 6350760 : m_uc[1][b] = u[i];
848 : 6350760 : m_pc[1][b] = prim[i];
849 : : }
850 : :
851 : : // if we have received all solution ghost contributions from neighboring
852 : : // chares (chares we communicate along chare-boundary faces with), and
853 : : // contributed our solution to these neighbors, proceed to limiting
854 [ + + ]: 383070 : if (++m_nreco == myGhosts()->m_sendGhost.size()) {
855 : 34215 : m_nreco = 0;
856 : 34215 : comreco_complete();
857 : : }
858 : 383070 : }
859 : :
860 : : void
861 : 37410 : DG::nodalExtrema()
862 : : // *****************************************************************************
863 : : // Compute nodal extrema at chare-boundary nodes. Extrema at internal nodes
864 : : // are calculated in limiter function.
865 : : // *****************************************************************************
866 : : {
867 : 37410 : auto d = Disc();
868 : 37410 : auto gid = d->Gid();
869 : : auto bid = d->Bid();
870 : 37410 : const auto rdof = g_inputdeck.get< tag::rdof >();
871 : 37410 : const auto ncomp = m_u.nprop() / rdof;
872 : 37410 : const auto nprim = m_p.nprop() / rdof;
873 : :
874 : : // Combine own and communicated contributions of unlimited solution, and
875 : : // if a p-adaptive algorithm is used, degrees of freedom in cells
876 [ + - ][ + + ]: 6425580 : for (const auto& [boundary, localtet] : myGhosts()->m_bid) {
877 : : Assert( m_uc[1][localtet].size() == m_u.nprop(), "ncomp size mismatch" );
878 : : Assert( m_pc[1][localtet].size() == m_p.nprop(), "ncomp size mismatch" );
879 [ + + ]: 167849685 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
880 : 161498925 : m_u(boundary,c) = m_uc[1][localtet][c];
881 : : }
882 [ + + ]: 33384885 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
883 : 27034125 : m_p(boundary,c) = m_pc[1][localtet][c];
884 : : }
885 : : }
886 : :
887 : : // Initialize nodal extrema vector
888 : : auto large = std::numeric_limits< tk::real >::max();
889 [ + + ]: 1288830 : for(std::size_t i = 0; i<bid.size(); i++)
890 : : {
891 [ + + ]: 9094320 : for (std::size_t c=0; c<ncomp; ++c)
892 : : {
893 [ + + ]: 31371600 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
894 : : {
895 : 23528700 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
896 : 23528700 : auto min_mark = max_mark + 1;
897 : 23528700 : m_uNodalExtrm[i][max_mark] = -large;
898 : 23528700 : m_uNodalExtrm[i][min_mark] = large;
899 : : }
900 : : }
901 [ + + ]: 3380070 : for (std::size_t c=0; c<nprim; ++c)
902 : : {
903 [ + + ]: 8514600 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
904 : : {
905 : 6385950 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
906 : 6385950 : auto min_mark = max_mark + 1;
907 : 6385950 : m_pNodalExtrm[i][max_mark] = -large;
908 : 6385950 : m_pNodalExtrm[i][min_mark] = large;
909 : : }
910 : : }
911 : : }
912 : :
913 : : // Evaluate the max/min value for the chare-boundary nodes
914 [ + + ]: 37410 : if(rdof > 4) {
915 [ + - ]: 13020 : evalNodalExtrmRefEl(ncomp, nprim, m_ndof_NodalExtrm, d->bndel(),
916 [ + - ][ + - ]: 6510 : myGhosts()->m_inpoel, gid, bid, m_u, m_p, m_uNodalExtrm, m_pNodalExtrm);
917 : : }
918 : :
919 : : // Communicate extrema at nodes to other chares on chare-boundary
920 [ + + ]: 37410 : if (d->NodeCommMap().empty()) // in serial we are done
921 [ + - ]: 3195 : comnodalExtrema_complete();
922 : : else // send nodal extrema to chare-boundary nodes to fellow chares
923 : : {
924 [ + + ]: 417285 : for (const auto& [c,n] : d->NodeCommMap()) {
925 [ + - ][ + - ]: 766140 : std::vector< std::vector< tk::real > > g1( n.size() ), g2( n.size() );
926 : : std::size_t j = 0;
927 [ + + ]: 3233220 : for (auto i : n)
928 : : {
929 : 2850150 : auto p = tk::cref_find(d->Bid(),i);
930 [ + - ]: 2850150 : g1[ j ] = m_uNodalExtrm[ p ];
931 [ + - ]: 2850150 : g2[ j++ ] = m_pNodalExtrm[ p ];
932 : : }
933 [ + - ][ + - ]: 766140 : thisProxy[c].comnodalExtrema( std::vector<std::size_t>(begin(n),end(n)),
[ + - ][ - + ]
934 : : g1, g2 );
935 : : }
936 : : }
937 [ + - ]: 37410 : ownnodalExtrema_complete();
938 : 37410 : }
939 : :
940 : : void
941 : 383070 : DG::comnodalExtrema( const std::vector< std::size_t >& gid,
942 : : const std::vector< std::vector< tk::real > >& G1,
943 : : const std::vector< std::vector< tk::real > >& G2 )
944 : : // *****************************************************************************
945 : : // Receive contributions to nodal extrema on chare-boundaries
946 : : //! \param[in] gid Global mesh node IDs at which we receive grad contributions
947 : : //! \param[in] G1 Partial contributions of extrema for conservative variables to
948 : : //! chare-boundary nodes
949 : : //! \param[in] G2 Partial contributions of extrema for primitive variables to
950 : : //! chare-boundary nodes
951 : : //! \details This function receives contributions to m_uNodalExtrm/m_pNodalExtrm
952 : : //! , which stores nodal extrems at mesh chare-boundary nodes. While
953 : : //! m_uNodalExtrm/m_pNodalExtrm stores own contributions, m_uNodalExtrmc
954 : : //! /m_pNodalExtrmc collects the neighbor chare contributions during
955 : : //! communication.
956 : : // *****************************************************************************
957 : : {
958 : : Assert( G1.size() == gid.size(), "Size mismatch" );
959 : : Assert( G2.size() == gid.size(), "Size mismatch" );
960 : :
961 : 383070 : const auto rdof = g_inputdeck.get< tag::rdof >();
962 : 383070 : const auto ncomp = m_u.nprop() / rdof;
963 : 383070 : const auto nprim = m_p.nprop() / rdof;
964 : :
965 [ + + ]: 3233220 : for (std::size_t i=0; i<gid.size(); ++i)
966 : : {
967 : : auto& u = m_uNodalExtrmc[gid[i]];
968 : 2850150 : auto& p = m_pNodalExtrmc[gid[i]];
969 [ + + ]: 20275440 : for (std::size_t c=0; c<ncomp; ++c)
970 : : {
971 [ + + ]: 69701160 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
972 : : {
973 : 52275870 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
974 : 52275870 : auto min_mark = max_mark + 1;
975 [ + + ][ + + ]: 52786939 : u[max_mark] = std::max( G1[i][max_mark], u[max_mark] );
976 : 52275870 : u[min_mark] = std::min( G1[i][min_mark], u[min_mark] );
977 : : }
978 : : }
979 [ + + ]: 7200000 : for (std::size_t c=0; c<nprim; ++c)
980 : : {
981 [ + + ]: 17399400 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
982 : : {
983 : 13049550 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
984 : 13049550 : auto min_mark = max_mark + 1;
985 [ - + ][ - + ]: 13049550 : p[max_mark] = std::max( G2[i][max_mark], p[max_mark] );
986 : 13049550 : p[min_mark] = std::min( G2[i][min_mark], p[min_mark] );
987 : : }
988 : : }
989 : : }
990 : :
991 [ + + ]: 383070 : if (++m_nnodalExtrema == Disc()->NodeCommMap().size())
992 : : {
993 : 34215 : m_nnodalExtrema = 0;
994 : 34215 : comnodalExtrema_complete();
995 : : }
996 : 383070 : }
997 : :
998 : 38075 : void DG::resizeNodalExtremac()
999 : : // *****************************************************************************
1000 : : // Resize the buffer vector of nodal extrema
1001 : : // *****************************************************************************
1002 : : {
1003 : 38075 : const auto rdof = g_inputdeck.get< tag::rdof >();
1004 : 38075 : const auto ncomp = m_u.nprop() / rdof;
1005 : 38075 : const auto nprim = m_p.nprop() / rdof;
1006 : :
1007 : 38075 : auto large = std::numeric_limits< tk::real >::max();
1008 [ + - ][ + + ]: 466232 : for (const auto& [c,n] : Disc()->NodeCommMap())
1009 : : {
1010 [ + + ][ + - ]: 3286530 : for (auto i : n) {
1011 : : auto& u = m_uNodalExtrmc[i];
1012 : : auto& p = m_pNodalExtrmc[i];
1013 [ + - ]: 2896448 : u.resize( 2*m_ndof_NodalExtrm*ncomp, large );
1014 [ + - ]: 2896448 : p.resize( 2*m_ndof_NodalExtrm*nprim, large );
1015 : :
1016 : : // Initialize the minimum nodal extrema
1017 [ + + ]: 11585792 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1018 : : {
1019 [ + + ]: 61756014 : for(std::size_t k = 0; k < ncomp; k++)
1020 : 53066670 : u[2*k*m_ndof_NodalExtrm+2*idof] = -large;
1021 [ + + ]: 21885156 : for(std::size_t k = 0; k < nprim; k++)
1022 : 13195812 : p[2*k*m_ndof_NodalExtrm+2*idof] = -large;
1023 : : }
1024 : : }
1025 : : }
1026 : 38075 : }
1027 : :
1028 : 6510 : void DG::evalNodalExtrmRefEl(
1029 : : const std::size_t ncomp,
1030 : : const std::size_t nprim,
1031 : : const std::size_t ndof_NodalExtrm,
1032 : : const std::vector< std::size_t >& bndel,
1033 : : const std::vector< std::size_t >& inpoel,
1034 : : const std::vector< std::size_t >& gid,
1035 : : const std::unordered_map< std::size_t, std::size_t >& bid,
1036 : : const tk::Fields& U,
1037 : : const tk::Fields& P,
1038 : : std::vector< std::vector<tk::real> >& uNodalExtrm,
1039 : : std::vector< std::vector<tk::real> >& pNodalExtrm )
1040 : : // *****************************************************************************
1041 : : // Compute the nodal extrema of ref el derivatives for chare-boundary nodes
1042 : : //! \param[in] ncomp Number of conservative variables
1043 : : //! \param[in] nprim Number of primitive variables
1044 : : //! \param[in] ndof_NodalExtrm Degree of freedom for nodal extrema
1045 : : //! \param[in] bndel List of elements contributing to chare-boundary nodes
1046 : : //! \param[in] inpoel Element-node connectivity for element e
1047 : : //! \param[in] gid Local->global node id map
1048 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
1049 : : //! global node ids (key)
1050 : : //! \param[in] U Vector of conservative variables
1051 : : //! \param[in] P Vector of primitive variables
1052 : : //! \param[in,out] uNodalExtrm Chare-boundary nodal extrema for conservative
1053 : : //! variables
1054 : : //! \param[in,out] pNodalExtrm Chare-boundary nodal extrema for primitive
1055 : : //! variables
1056 : : // *****************************************************************************
1057 : : {
1058 : 6510 : const auto rdof = g_inputdeck.get< tag::rdof >();
1059 : :
1060 [ + + ]: 599550 : for (auto e : bndel)
1061 : : {
1062 : : // access node IDs
1063 : : const std::vector<std::size_t> N
1064 [ + - ]: 593040 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
1065 : :
1066 : : // Loop over nodes of element e
1067 [ + + ]: 2965200 : for(std::size_t ip=0; ip<4; ++ip)
1068 : : {
1069 : 2372160 : auto i = bid.find( gid[N[ip]] );
1070 [ + + ]: 2372160 : if (i != end(bid)) // If ip is the chare boundary point
1071 : : {
1072 : : // If DG(P2) is applied, find the nodal extrema of the gradients of
1073 : : // conservative/primitive variables in the reference element
1074 : :
1075 : : // Vector used to store the first order derivatives for both
1076 : : // conservative and primitive variables
1077 [ + - ][ - - ]: 1623930 : std::vector< std::array< tk::real, 3 > > gradc(ncomp, {0.0, 0.0, 0.0});
1078 [ + - ][ - - ]: 1623930 : std::vector< std::array< tk::real, 3 > > gradp(ncomp, {0.0, 0.0, 0.0});
1079 : :
1080 : : // Derivatives of the Dubiner basis
1081 : 1623930 : std::array< tk::real, 3 > center {{0.25, 0.25, 0.25}};
1082 [ + - ]: 1623930 : auto dBdxi = tk::eval_dBdxi(rdof, center);
1083 : :
1084 : : // Evaluate the first order derivative
1085 [ + + ]: 9743580 : for(std::size_t icomp = 0; icomp < ncomp; icomp++)
1086 : : {
1087 : 8119650 : auto mark = icomp * rdof;
1088 [ + + ]: 32478600 : for(std::size_t idir = 0; idir < 3; idir++)
1089 : : {
1090 : 24358950 : gradc[icomp][idir] = 0;
1091 [ + + ]: 243589500 : for(std::size_t idof = 1; idof < rdof; idof++)
1092 : 219230550 : gradc[icomp][idir] += U(e, mark+idof) * dBdxi[idir][idof];
1093 : : }
1094 : : }
1095 [ - + ]: 1623930 : for(std::size_t icomp = 0; icomp < nprim; icomp++)
1096 : : {
1097 : 0 : auto mark = icomp * rdof;
1098 [ - - ]: 0 : for(std::size_t idir = 0; idir < 3; idir++)
1099 : : {
1100 : 0 : gradp[icomp][idir] = 0;
1101 [ - - ]: 0 : for(std::size_t idof = 1; idof < rdof; idof++)
1102 : 0 : gradp[icomp][idir] += P(e, mark+idof) * dBdxi[idir][idof];
1103 : : }
1104 : : }
1105 : :
1106 : : // Store the extrema for the gradients
1107 [ + + ]: 9743580 : for (std::size_t c=0; c<ncomp; ++c)
1108 : : {
1109 [ + + ]: 32478600 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
1110 : : {
1111 : 24358950 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1112 : 24358950 : auto min_mark = max_mark + 1;
1113 [ + + ]: 24358950 : auto& ex = uNodalExtrm[i->second];
1114 [ + + ][ + + ]: 30437923 : ex[max_mark] = std::max(ex[max_mark], gradc[c][idof]);
1115 : 24358950 : ex[min_mark] = std::min(ex[min_mark], gradc[c][idof]);
1116 : : }
1117 : : }
1118 [ - + ]: 1623930 : for (std::size_t c=0; c<nprim; ++c)
1119 : : {
1120 [ - - ]: 0 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
1121 : : {
1122 : 0 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1123 : 0 : auto min_mark = max_mark + 1;
1124 [ - - ]: 0 : auto& ex = pNodalExtrm[i->second];
1125 [ - - ][ - - ]: 0 : ex[max_mark] = std::max(ex[max_mark], gradp[c][idof]);
1126 : 0 : ex[min_mark] = std::min(ex[min_mark], gradp[c][idof]);
1127 : : }
1128 : : }
1129 : : }
1130 : : }
1131 : : }
1132 : 6510 : }
1133 : :
1134 : : void
1135 : 37410 : DG::lim()
1136 : : // *****************************************************************************
1137 : : // Compute limiter function
1138 : : // *****************************************************************************
1139 : : {
1140 : 37410 : auto d = Disc();
1141 : 37410 : const auto rdof = g_inputdeck.get< tag::rdof >();
1142 : 37410 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1143 : 37410 : const auto ncomp = m_u.nprop() / rdof;
1144 : 37410 : const auto nprim = m_p.nprop() / rdof;
1145 : :
1146 : : // Combine own and communicated contributions to nodal extrema
1147 [ + + ]: 1288830 : for (const auto& [gid,g] : m_uNodalExtrmc) {
1148 : 1251420 : auto bid = tk::cref_find( d->Bid(), gid );
1149 [ + + ]: 9094320 : for (ncomp_t c=0; c<ncomp; ++c)
1150 : : {
1151 [ + + ]: 31371600 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1152 : : {
1153 : 23528700 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1154 : 23528700 : auto min_mark = max_mark + 1;
1155 : 23528700 : m_uNodalExtrm[bid][max_mark] =
1156 [ + + ][ + + ]: 24134054 : std::max(g[max_mark], m_uNodalExtrm[bid][max_mark]);
1157 : 23528700 : m_uNodalExtrm[bid][min_mark] =
1158 : 23528700 : std::min(g[min_mark], m_uNodalExtrm[bid][min_mark]);
1159 : : }
1160 : : }
1161 : : }
1162 [ + + ]: 1288830 : for (const auto& [gid,g] : m_pNodalExtrmc) {
1163 : 1251420 : auto bid = tk::cref_find( d->Bid(), gid );
1164 [ + + ]: 3380070 : for (ncomp_t c=0; c<nprim; ++c)
1165 : : {
1166 [ + + ]: 8514600 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1167 : : {
1168 : 6385950 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1169 : 6385950 : auto min_mark = max_mark + 1;
1170 : 6385950 : m_pNodalExtrm[bid][max_mark] =
1171 [ - + ][ - + ]: 6385950 : std::max(g[max_mark], m_pNodalExtrm[bid][max_mark]);
1172 : 6385950 : m_pNodalExtrm[bid][min_mark] =
1173 : 6385950 : std::min(g[min_mark], m_pNodalExtrm[bid][min_mark]);
1174 : : }
1175 : : }
1176 : : }
1177 : :
1178 : : // clear gradients receive buffer
1179 : 37410 : tk::destroy(m_uNodalExtrmc);
1180 : 37410 : tk::destroy(m_pNodalExtrmc);
1181 : :
1182 [ + + ]: 37410 : if (rdof > 1) {
1183 : 47970 : g_dgpde[d->MeshId()].limit( d->T(), pref, myGhosts()->m_geoFace,
1184 : 23985 : myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_esup,
1185 : 23985 : myGhosts()->m_inpoel, myGhosts()->m_coord, m_ndof, d->Gid(),
1186 : 23985 : d->Bid(), m_uNodalExtrm, m_pNodalExtrm, m_mtInv, m_u, m_p,
1187 : 23985 : m_shockmarker );
1188 : :
1189 [ + + ]: 23985 : if (g_inputdeck.get< tag::limsol_projection >())
1190 : 40620 : g_dgpde[d->MeshId()].CPL(m_p, myGhosts()->m_geoElem,
1191 : 20310 : myGhosts()->m_inpoel, myGhosts()->m_coord, m_u,
1192 : 20310 : myGhosts()->m_fd.Esuel().size()/4);
1193 : : }
1194 : :
1195 : : // Send limited solution to neighboring chares
1196 [ + + ]: 37410 : if (myGhosts()->m_sendGhost.empty())
1197 : 3195 : comlim_complete();
1198 : : else
1199 [ + + ]: 451500 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
1200 [ + - ]: 383070 : std::vector< std::size_t > tetid( ghostdata.size() );
1201 [ + - ][ + - ]: 766140 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
1202 [ + - ]: 383070 : prim( ghostdata.size() );
1203 : : std::vector< std::size_t > ndof;
1204 : : std::size_t j = 0;
1205 [ + + ]: 6733830 : for(const auto& i : ghostdata) {
1206 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
1207 : : "Sending limiter ghost data" );
1208 [ + - ]: 6350760 : tetid[j] = i;
1209 [ + - ][ - + ]: 6350760 : u[j] = m_u[i];
1210 [ + - ][ - + ]: 6350760 : prim[j] = m_p[i];
1211 : 6350760 : ++j;
1212 : : }
1213 [ + - ][ + - ]: 766140 : thisProxy[ cid ].comlim( thisIndex, tetid, u, prim );
1214 : : }
1215 : :
1216 : 37410 : ownlim_complete();
1217 : 37410 : }
1218 : :
1219 : : void
1220 : 1770 : DG::refine_ndof()
1221 : : // *****************************************************************************
1222 : : // p-refine all elements that are adjacent to p-refined elements
1223 : : //! \details This function p-refines all the neighbors of an element that has
1224 : : //! been p-refined as a result of an error indicator.
1225 : : // *****************************************************************************
1226 : : {
1227 : 1770 : auto d = Disc();
1228 : : const auto& coord = d->Coord();
1229 : 1770 : const auto& inpoel = d->Inpoel();
1230 : 1770 : const auto npoin = coord[0].size();
1231 : 1770 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
1232 : 1770 : std::vector<std::size_t> node_ndof(npoin, 1);
1233 : :
1234 : : // Mark the max ndof for each node and store in node_ndof
1235 [ + + ]: 188540 : for(std::size_t ip=0; ip<npoin; ip++)
1236 : : {
1237 [ + - ]: 186770 : const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
1238 [ + + ]: 2641810 : for(auto er : pesup)
1239 [ + + ]: 2536278 : node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
1240 : : }
1241 : :
1242 [ + + ]: 414490 : for(std::size_t e = 0; e < nelem; e++)
1243 : : {
1244 : : // Find if any node of this element has p1/p2 ndofs
1245 : : std::size_t counter_p2(0);
1246 : : std::size_t counter_p1(0);
1247 [ + + ]: 2063600 : for(std::size_t inode = 0; inode < 4; inode++)
1248 : : {
1249 [ + + ]: 1650880 : auto node = inpoel[4*e+inode];
1250 [ + + ]: 1650880 : if(node_ndof[node] == 10)
1251 : 320808 : counter_p2++;
1252 [ + + ]: 1330072 : else if (node_ndof[node] == 4)
1253 : 180140 : counter_p1++;
1254 : : }
1255 : :
1256 : : // If there is at least one node with p1/p2 ndofs, all of the elements
1257 : : // around this node are refined to p1/p2.
1258 [ + + ][ + + ]: 412720 : if(counter_p2 > 0 && m_ndof[e] < 10)
1259 : : {
1260 [ + + ]: 16717 : if(m_ndof[e] == 4)
1261 : 15693 : m_ndof[e] = 10;
1262 [ + + ]: 16717 : if(m_ndof[e] == 1)
1263 : 1024 : m_ndof[e] = 4;
1264 : : }
1265 [ + + ][ + + ]: 396003 : else if(counter_p1 > 0 && m_ndof[e] < 4)
1266 : 13452 : m_ndof[e] = 4;
1267 : : }
1268 : 1770 : }
1269 : :
1270 : 1770 : void DG::smooth_ndof()
1271 : : // *****************************************************************************
1272 : : // Smooth the refined ndof distribution to avoid zigzag refinement
1273 : : // *****************************************************************************
1274 : : {
1275 : 1770 : auto d = Disc();
1276 : 1770 : const auto& inpoel = d->Inpoel();
1277 : : const auto& coord = d->Coord();
1278 : 1770 : const auto npoin = coord[0].size();
1279 : 1770 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
1280 : 1770 : std::vector<std::size_t> node_ndof(npoin, 1);
1281 : :
1282 : : // Mark the max ndof for each node and store in node_ndof
1283 [ + + ]: 188540 : for(std::size_t ip=0; ip<npoin; ip++)
1284 : : {
1285 [ + - ]: 186770 : const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
1286 [ + + ]: 2641810 : for(auto er : pesup)
1287 [ + + ]: 2542408 : node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
1288 : : }
1289 : :
1290 [ + + ]: 414490 : for(std::size_t e = 0; e < nelem; e++)
1291 : : {
1292 : : // Find if any node of this element has p1/p2 ndofs
1293 : : std::size_t counter_p2(0);
1294 : : std::size_t counter_p1(0);
1295 [ + + ]: 2063600 : for(std::size_t inode = 0; inode < 4; inode++)
1296 : : {
1297 [ + + ]: 1650880 : auto node = inpoel[4*e+inode];
1298 [ + + ]: 1650880 : if(node_ndof[node] == 10)
1299 : 382503 : counter_p2++;
1300 [ + + ]: 1268377 : else if (node_ndof[node] == 4)
1301 : 182564 : counter_p1++;
1302 : : }
1303 : :
1304 : : // If all the nodes in the element are p1/p2, this element is refined to
1305 : : // p1/p2.
1306 [ + + ][ + + ]: 412720 : if(counter_p2 == 4 && m_ndof[e] == 4)
1307 : 1469 : m_ndof[e] = 10;
1308 [ + + ][ + + ]: 411251 : else if(counter_p1 == 4 && m_ndof[e] == 1)
1309 : 1553 : m_ndof[e] = 4;
1310 : : }
1311 : 1770 : }
1312 : :
1313 : : void
1314 : 383070 : DG::comlim( int fromch,
1315 : : const std::vector< std::size_t >& tetid,
1316 : : const std::vector< std::vector< tk::real > >& u,
1317 : : const std::vector< std::vector< tk::real > >& prim )
1318 : : // *****************************************************************************
1319 : : // Receive chare-boundary limiter ghost data from neighboring chares
1320 : : //! \param[in] fromch Sender chare id
1321 : : //! \param[in] tetid Ghost tet ids we receive solution data for
1322 : : //! \param[in] u Limited high-order solution
1323 : : //! \param[in] prim Limited high-order primitive quantities
1324 : : //! \details This function receives contributions to the limited solution from
1325 : : //! fellow chares.
1326 : : // *****************************************************************************
1327 : : {
1328 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comlim()" );
1329 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comlim()" );
1330 : :
1331 : : // Find local-to-ghost tet id map for sender chare
1332 : 383070 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
1333 : :
1334 [ + + ]: 6733830 : for (std::size_t i=0; i<tetid.size(); ++i) {
1335 : 6350760 : auto j = tk::cref_find( n, tetid[i] );
1336 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
1337 : : "Receiving solution non-ghost data" );
1338 : 6350760 : auto b = tk::cref_find( myGhosts()->m_bid, j );
1339 : : Assert( b < m_uc[2].size(), "Indexing out of bounds" );
1340 : : Assert( b < m_pc[2].size(), "Indexing out of bounds" );
1341 : 6350760 : m_uc[2][b] = u[i];
1342 : 6350760 : m_pc[2][b] = prim[i];
1343 : : }
1344 : :
1345 : : // if we have received all solution ghost contributions from neighboring
1346 : : // chares (chares we communicate along chare-boundary faces with), and
1347 : : // contributed our solution to these neighbors, proceed to limiting
1348 [ + + ]: 383070 : if (++m_nlim == myGhosts()->m_sendGhost.size()) {
1349 : 34215 : m_nlim = 0;
1350 : 34215 : comlim_complete();
1351 : : }
1352 : 383070 : }
1353 : :
1354 : : void
1355 : 37410 : DG::dt()
1356 : : // *****************************************************************************
1357 : : // Compute time step size
1358 : : // *****************************************************************************
1359 : : {
1360 : 37410 : auto d = Disc();
1361 : :
1362 : : // Combine own and communicated contributions of limited solution and degrees
1363 : : // of freedom in cells (if p-adaptive)
1364 [ + + ]: 6425580 : for (const auto& b : myGhosts()->m_bid) {
1365 : : Assert( m_uc[2][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
1366 : : Assert( m_pc[2][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
1367 [ + + ]: 167849685 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1368 : 161498925 : m_u(b.first,c) = m_uc[2][b.second][c];
1369 : : }
1370 [ + + ]: 33384885 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
1371 : 27034125 : m_p(b.first,c) = m_pc[2][b.second][c];
1372 : : }
1373 : : }
1374 : :
1375 : 37410 : auto mindt = std::numeric_limits< tk::real >::max();
1376 : :
1377 [ + + ]: 37410 : if (m_stage == 0)
1378 : : {
1379 [ + + ]: 12470 : auto const_dt = g_inputdeck.get< tag::dt >();
1380 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1381 : :
1382 : : // use constant dt if configured
1383 [ + + ]: 12470 : if (std::abs(const_dt) > eps) {
1384 : :
1385 : 9725 : mindt = const_dt;
1386 : :
1387 : : } else { // compute dt based on CFL
1388 : :
1389 : : // find the minimum dt across all PDEs integrated
1390 : : auto eqdt =
1391 : 5490 : g_dgpde[d->MeshId()].dt( myGhosts()->m_coord, myGhosts()->m_inpoel,
1392 : 2745 : myGhosts()->m_fd,
1393 : 2745 : myGhosts()->m_geoFace, myGhosts()->m_geoElem, m_ndof, m_u, m_p,
1394 : 2745 : myGhosts()->m_fd.Esuel().size()/4 );
1395 [ + - ]: 2745 : if (eqdt < mindt) mindt = eqdt;
1396 : :
1397 : : // time-step suppression for unsteady problems
1398 : : tk::real coeff(1.0);
1399 : 2745 : auto ramp_steps = g_inputdeck.get< tag::cfl_ramping_steps >();
1400 [ - + ][ - - ]: 2745 : if (g_inputdeck.get< tag::cfl_ramping >() && d->It() < ramp_steps)
1401 : 0 : coeff = 1.0/static_cast< tk::real >(ramp_steps)
1402 : 0 : * static_cast< tk::real >(d->It()+1);
1403 : :
1404 : 2745 : mindt *= coeff * g_inputdeck.get< tag::cfl >();
1405 : : }
1406 : : }
1407 : : else
1408 : : {
1409 : 24940 : mindt = d->Dt();
1410 : : }
1411 : :
1412 : : // Resize the buffer vector of nodal extrema
1413 : 37410 : resizeNodalExtremac();
1414 : :
1415 : : // Set up the reduction target for finding minimum dt across chares.
1416 : : // 1. If implicit solver is used, first invoke the solver object via
1417 : : // appropriate entry methods, and then proceed to solve.
1418 : : // 2. If explicit, directly proceed to solve.
1419 : : CkCallback minDtDone;
1420 [ + - ]: 37410 : if (!g_inputdeck.get< tag::implicit_timestepping >())
1421 : 37410 : minDtDone = CkCallback(CkReductionTarget(DG,solve), thisProxy);
1422 : : else
1423 : 0 : minDtDone = CkCallback(CkReductionTarget(DG,initializeLinearSystem),
1424 : : thisProxy);
1425 : :
1426 : : // Contribute to minimum dt across all chares then advance to next step
1427 [ + - ]: 37410 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double, minDtDone );
1428 : 37410 : }
1429 : :
1430 : : void
1431 : 0 : DG::initializeLinearSystem( tk::real newdt )
1432 : : // *****************************************************************************
1433 : : // Initialize the linear solver via the interface BiCG::init()
1434 : : //! \param[in] newdt Size of this new time step
1435 : : // *****************************************************************************
1436 : : {
1437 : 0 : auto d = Disc();
1438 : :
1439 : : // Set new time step size
1440 [ - - ]: 0 : if (m_stage == 0) d->setdt( newdt );
1441 : :
1442 : : // Initialize linear solver, and route to solveLinearSystem()
1443 : : // TODO: linear solver:
1444 : : // 1. jacobian computation (call to e.g. g_dgpde[d->MeshId()].computeJacobian)
1445 : : // 2. the following call is just a stand-in/example- verify correctness
1446 [ - - ][ - - ]: 0 : d->ImplicitSolver()[ thisIndex ].init( m_u.flat(), {}, {}, 1,
[ - - ][ - - ]
1447 [ - - ][ - - ]: 0 : CkCallback(CkIndex_DG::solveLinearSystem(), thisProxy[thisIndex]) );
[ - - ][ - - ]
[ - - ][ - - ]
1448 : 0 : }
1449 : :
1450 : : void
1451 : 0 : DG::solveLinearSystem()
1452 : : // *****************************************************************************
1453 : : // Solve the linear system via the interface BiCG::solve()
1454 : : // *****************************************************************************
1455 : : {
1456 : 0 : auto d = Disc();
1457 : :
1458 : : // Get new time step size to pass along to solve()
1459 : : auto dt = d->Dt();
1460 : :
1461 : : // Solve linear system, and route to solve()
1462 : : // TODO: linear solver:
1463 : : // the following call is just a stand-in/example- verify correctness
1464 [ - - ]: 0 : d->ImplicitSolver()[ thisIndex ].solve(
1465 : : g_inputdeck.get< tag::ale, tag::maxit >(),
1466 : : g_inputdeck.get< tag::residual >(),
1467 [ - - ][ - - ]: 0 : CkCallback(CkIndex_DG::solve(dt), thisProxy[thisIndex]) );
[ - - ][ - - ]
[ - - ][ - - ]
1468 : 0 : }
1469 : :
1470 : : void
1471 : 37410 : DG::solve( tk::real newdt )
1472 : : // *****************************************************************************
1473 : : // Compute right-hand side of discrete transport equations
1474 : : //! \param[in] newdt Size of this new time step
1475 : : // *****************************************************************************
1476 : : {
1477 : 37410 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1478 : :
1479 : : // Enable SDAG wait for building the solution vector during the next stage
1480 [ + - ]: 37410 : thisProxy[ thisIndex ].wait4sol();
1481 [ + + ][ + - ]: 42720 : if (pref) thisProxy[ thisIndex ].wait4refine();
1482 [ + - ]: 37410 : thisProxy[ thisIndex ].wait4smooth();
1483 [ + - ]: 37410 : thisProxy[ thisIndex ].wait4reco();
1484 [ + - ]: 37410 : thisProxy[ thisIndex ].wait4nodalExtrema();
1485 [ + - ]: 37410 : thisProxy[ thisIndex ].wait4lim();
1486 [ + - ]: 37410 : thisProxy[ thisIndex ].wait4nod();
1487 : :
1488 : 37410 : auto d = Disc();
1489 : 37410 : const auto rdof = g_inputdeck.get< tag::rdof >();
1490 : 37410 : const auto ndof = g_inputdeck.get< tag::ndof >();
1491 : 37410 : const auto neq = m_u.nprop()/rdof;
1492 : :
1493 : : // Set new time step size. If implicit solver, time step has already been
1494 : : // set in initializeImplicitSystem()
1495 [ + + ][ + - ]: 37410 : if (m_stage == 0 && !g_inputdeck.get< tag::implicit_timestepping >())
1496 : 12470 : d->setdt( newdt );
1497 : :
1498 : : // Update Un
1499 [ + + ]: 37410 : if (m_stage == 0) m_un = m_u;
1500 : :
1501 : : // Explicit or IMEX
1502 : 37410 : const auto imex_runge_kutta = g_inputdeck.get< tag::imex_runge_kutta >();
1503 : 37410 : const auto implicit_ts = g_inputdeck.get< tag::implicit_timestepping >();
1504 : :
1505 : : // physical time at time-stage for computing exact source terms
1506 : 37410 : tk::real physT(d->T());
1507 [ + + ]: 37410 : if (m_stage == 1) {
1508 : 12470 : physT += d->Dt();
1509 : : }
1510 [ + + ]: 24940 : else if (m_stage == 2) {
1511 : 12470 : physT += 0.5*d->Dt();
1512 : : }
1513 : :
1514 [ - + ]: 37410 : if (imex_runge_kutta) {
1515 [ - - ]: 0 : if (m_stage == 0)
1516 : : {
1517 : : // Save previous rhs
1518 : : m_rhsprev = m_rhs;
1519 : : // Initialize m_stiffrhs to zero
1520 : : m_stiffrhs.fill(0.0);
1521 : : m_stiffrhsprev.fill(0.0);
1522 : : }
1523 : : }
1524 : :
1525 : 74820 : g_dgpde[d->MeshId()].rhs( physT, pref, myGhosts()->m_geoFace,
1526 : 37410 : myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_inpoel, m_boxelems,
1527 : 37410 : myGhosts()->m_coord, m_u, m_p, m_ndof, d->Dt(), m_rhs );
1528 : :
1529 [ - + ]: 37410 : if (imex_runge_kutta) {
1530 : : // Implicit-Explicit time-stepping using RK3 to discretize time-derivative
1531 : 0 : DG::imex_integrate();
1532 : : }
1533 [ + - ]: 37410 : else if (implicit_ts) {
1534 : : // Implicit time-stepping using BDF1 to discretize time-derivative
1535 : 0 : DG::BDF1_integrate();
1536 : : }
1537 : : else {
1538 : : // Explicit time-stepping using RK3 to discretize time-derivative
1539 [ + + ]: 15874950 : for(std::size_t e=0; e<myGhosts()->m_nunk; ++e) {
1540 : 15837540 : auto vole = myGhosts()->m_geoElem(e,0);
1541 [ + + ]: 107664435 : for(std::size_t c=0; c<neq; ++c)
1542 : : {
1543 [ + + ]: 351184560 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
1544 : : {
1545 [ + + ]: 259357665 : if(k < m_ndof[e]) {
1546 : 174605205 : auto rmark = c*rdof+k;
1547 : 174605205 : auto mark = c*ndof+k;
1548 [ + + ]: 174605205 : m_u(e, rmark) = rkcoef[0][m_stage] * m_un(e, rmark)
1549 : 174605205 : + rkcoef[1][m_stage] * ( m_u(e, rmark)
1550 : 174605205 : + d->Dt() * m_rhs(e, mark)/ (vole*mass_dubiner[k]));
1551 [ + + ]: 174605205 : if(fabs(m_u(e, rmark)) < 1e-16)
1552 : 28121194 : m_u(e, rmark) = 0;
1553 : : }
1554 : : }
1555 : : }
1556 : : }
1557 : : }
1558 : :
1559 [ + + ]: 15874950 : for(std::size_t e=0; e<myGhosts()->m_nunk; ++e)
1560 [ + + ]: 107664435 : for(std::size_t c=0; c<neq; ++c)
1561 : : {
1562 : : // zero out unused/reconstructed dofs of equations using reduced dofs
1563 : : // (see DGMultiMat::numEquationDofs())
1564 [ + + ]: 91826895 : if (m_numEqDof[c] < rdof) {
1565 [ + + ]: 95142960 : for (std::size_t k=m_numEqDof[c]; k<rdof; ++k)
1566 : : {
1567 : 71357220 : auto rmark = c*rdof+k;
1568 : 71357220 : m_u(e, rmark) = 0.0;
1569 : : }
1570 : : }
1571 : : }
1572 : :
1573 : : // Update primitives based on the evolved solution
1574 : 37410 : g_dgpde[d->MeshId()].updateInterfaceCells( m_u,
1575 : 37410 : myGhosts()->m_fd.Esuel().size()/4, m_ndof, m_interface );
1576 : 37410 : g_dgpde[d->MeshId()].updatePrimitives( m_u, myGhosts()->m_geoElem, m_p,
1577 : 37410 : myGhosts()->m_fd.Esuel().size()/4, m_ndof );
1578 [ + - ]: 37410 : if (!g_inputdeck.get< tag::accuracy_test >()) {
1579 : 37410 : g_dgpde[d->MeshId()].cleanTraceMaterial( physT, myGhosts()->m_geoElem, m_u,
1580 : 37410 : m_p, myGhosts()->m_fd.Esuel().size()/4 );
1581 : : }
1582 : :
1583 [ + + ]: 37410 : if (m_stage < m_nstage-1) {
1584 : :
1585 : : // continue with next time step stage
1586 : 24940 : stage();
1587 : :
1588 : : } else {
1589 : :
1590 : : // Increase number of iterations and physical time
1591 : 12470 : d->next();
1592 : :
1593 : : // Compute diagnostics, e.g., residuals
1594 : 12470 : auto diag_computed = m_diag.compute( *d,
1595 : 12470 : m_u.nunk()-myGhosts()->m_fd.Esuel().size()/4, myGhosts()->m_geoElem,
1596 : 12470 : m_ndof, m_u, m_un );
1597 : :
1598 : : // Continue to mesh refinement (if configured)
1599 [ + + ][ + - ]: 14677 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 0.0 ) );
[ + - ]
1600 : :
1601 : : }
1602 : 37410 : }
1603 : :
1604 : : void
1605 : 12470 : DG::refine( [[maybe_unused]] const std::vector< tk::real >& l2res )
1606 : : // *****************************************************************************
1607 : : // Optionally refine/derefine mesh
1608 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1609 : : //! computed across the whole problem
1610 : : // *****************************************************************************
1611 : : {
1612 : 12470 : auto d = Disc();
1613 : :
1614 : 12470 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
1615 : 12470 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
1616 : :
1617 : : // if t>0 refinement enabled and we hit the dtref frequency
1618 [ - + ][ - - ]: 12470 : if (dtref && !(d->It() % dtfreq)) { // refine
1619 : :
1620 : 0 : d->startvol();
1621 [ - - ][ - - ]: 0 : d->Ref()->dtref( myGhosts()->m_fd.Bface(), {},
[ - - ][ - - ]
1622 : 0 : tk::remap(myGhosts()->m_fd.Triinpoel(),d->Gid()) );
1623 : 0 : d->refined() = 1;
1624 : :
1625 : : } else { // do not refine
1626 : :
1627 : 12470 : d->refined() = 0;
1628 : 12470 : stage();
1629 : :
1630 : : }
1631 : 12470 : }
1632 : :
1633 : : void
1634 : 0 : DG::resizePostAMR(
1635 : : const std::vector< std::size_t >& /*ginpoel*/,
1636 : : const tk::UnsMesh::Chunk& chunk,
1637 : : const tk::UnsMesh::Coords& coord,
1638 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
1639 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
1640 : : const std::set< std::size_t >& removedNodes,
1641 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
1642 : : const tk::NodeCommMap& nodeCommMap,
1643 : : const std::map< int, std::vector< std::size_t > >& bface,
1644 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
1645 : : const std::vector< std::size_t >& triinpoel,
1646 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid )
1647 : : // *****************************************************************************
1648 : : // Receive new mesh from Refiner
1649 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1650 : : //! \param[in] coord New mesh node coordinates
1651 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1652 : : //! \param[in] removedNodes Newly removed mesh node local ids
1653 : : //! \param[in] amrNodeMap Node id map after amr (local ids)
1654 : : //! \param[in] nodeCommMap New node communication map
1655 : : //! \param[in] bface Boundary-faces mapped to side set ids
1656 : : //! \param[in] triinpoel Boundary-face connectivity
1657 : : //! \param[in] elemblockid Local tet ids associated with mesh block ids
1658 : : // *****************************************************************************
1659 : : {
1660 : 0 : auto d = Disc();
1661 : :
1662 : : // Set flag that indicates that we are during time stepping
1663 : 0 : m_initial = 0;
1664 : 0 : myGhosts()->m_initial = 0;
1665 : :
1666 : : // Zero field output iteration count between two mesh refinement steps
1667 : 0 : d->Itf() = 0;
1668 : :
1669 : : // Increase number of iterations with mesh refinement
1670 : 0 : ++d->Itr();
1671 : :
1672 : : // Save old number of elements
1673 : 0 : [[maybe_unused]] auto old_nelem = myGhosts()->m_inpoel.size()/4;
1674 : :
1675 : : // Resize mesh data structures
1676 : 0 : d->resizePostAMR( chunk, coord, amrNodeMap, nodeCommMap, removedNodes,
1677 : : elemblockid );
1678 : :
1679 : : // Update state
1680 : 0 : myGhosts()->m_inpoel = d->Inpoel();
1681 : 0 : myGhosts()->m_coord = d->Coord();
1682 : 0 : auto nelem = myGhosts()->m_inpoel.size()/4;
1683 : : m_p.resize( nelem );
1684 : : m_u.resize( nelem );
1685 : : m_un.resize( nelem );
1686 : : m_rhs.resize( nelem );
1687 : : m_rhsprev.resize( nelem );
1688 : : m_stiffrhs.resize( nelem );
1689 : : m_stiffrhsprev.resize( nelem );
1690 [ - - ][ - - ]: 0 : m_uNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
1691 [ - - ]: 0 : m_ndof_NodalExtrm*g_inputdeck.get< tag::ncomp >() ) );
1692 [ - - ][ - - ]: 0 : m_pNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
1693 [ - - ]: 0 : m_ndof_NodalExtrm*m_p.nprop()/g_inputdeck.get< tag::rdof >()));
1694 : :
1695 : : // Resize the buffer vector of nodal extrema
1696 : 0 : resizeNodalExtremac();
1697 : :
1698 [ - - ][ - - ]: 0 : myGhosts()->m_fd = FaceData( myGhosts()->m_inpoel, bface,
[ - - ][ - - ]
[ - - ]
1699 : 0 : tk::remap(triinpoel,d->Lid()) );
1700 : :
1701 [ - - ]: 0 : myGhosts()->m_geoFace =
1702 : 0 : tk::Fields( tk::genGeoFaceTri( myGhosts()->m_fd.Nipfac(),
1703 : 0 : myGhosts()->m_fd.Inpofa(), coord ) );
1704 [ - - ]: 0 : myGhosts()->m_geoElem = tk::Fields( tk::genGeoElemTet( myGhosts()->m_inpoel,
1705 : 0 : coord ) );
1706 : :
1707 : 0 : myGhosts()->m_nfac = myGhosts()->m_fd.Inpofa().size()/3;
1708 : 0 : myGhosts()->m_nunk = nelem;
1709 : 0 : m_npoin = coord[0].size();
1710 : 0 : myGhosts()->m_bndFace.clear();
1711 : 0 : myGhosts()->m_exptGhost.clear();
1712 : 0 : myGhosts()->m_sendGhost.clear();
1713 : 0 : myGhosts()->m_ghost.clear();
1714 : 0 : myGhosts()->m_esup.clear();
1715 : :
1716 : : // Update solution on new mesh, P0 (cell center value) only for now
1717 : : m_un = m_u;
1718 : : auto pn = m_p;
1719 : 0 : auto unprop = m_u.nprop();
1720 : : auto pnprop = m_p.nprop();
1721 [ - - ]: 0 : for (const auto& [child,parent] : addedTets) {
1722 : : Assert( child < nelem, "Indexing out of new solution vector" );
1723 : : Assert( parent < old_nelem, "Indexing out of old solution vector" );
1724 [ - - ]: 0 : for (std::size_t i=0; i<unprop; ++i) m_u(child,i) = m_un(parent,i);
1725 [ - - ]: 0 : for (std::size_t i=0; i<pnprop; ++i) m_p(child,i) = pn(parent,i);
1726 : : }
1727 : : m_un = m_u;
1728 : :
1729 : : // Resize communication buffers
1730 [ - - ][ - - ]: 0 : m_ghosts[thisIndex].resizeComm();
[ - - ][ - - ]
1731 : 0 : }
1732 : :
1733 : : bool
1734 : 9923 : DG::fieldOutput() const
1735 : : // *****************************************************************************
1736 : : // Decide wether to output field data
1737 : : //! \return True if field data is output in this step
1738 : : // *****************************************************************************
1739 : : {
1740 : 9923 : auto d = Disc();
1741 : :
1742 : : // Output field data
1743 [ + + ][ + - ]: 9923 : return d->fielditer() or d->fieldtime() or d->fieldrange() or d->finished();
[ + - ][ + + ]
1744 : : }
1745 : :
1746 : : bool
1747 : 1278 : DG::refinedOutput() const
1748 : : // *****************************************************************************
1749 : : // Decide if we write field output using a refined mesh
1750 : : //! \return True if field output will use a refined mesh
1751 : : // *****************************************************************************
1752 : : {
1753 [ + + ]: 1278 : return g_inputdeck.get< tag::field_output, tag::refined >() &&
1754 [ - + ]: 33 : g_inputdeck.get< tag::scheme >() != ctr::SchemeType::DGP0;
1755 : : }
1756 : :
1757 : : void
1758 : 1278 : DG::writeFields(
1759 : : CkCallback c,
1760 : : const std::unordered_map< std::size_t, std::size_t >& addedTets )
1761 : : // *****************************************************************************
1762 : : // Output mesh field data
1763 : : //! \param[in] c Function to continue with after the write
1764 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1765 : : // *****************************************************************************
1766 : : {
1767 : 1278 : auto d = Disc();
1768 : :
1769 : : const auto& inpoel = std::get< 0 >( m_outmesh.chunk );
1770 : 2556 : auto esup = tk::genEsup( inpoel, 4 );
1771 : 1278 : auto nelem = inpoel.size() / 4;
1772 : :
1773 : : // Combine own and communicated contributions and finish averaging of node
1774 : : // field output in chare boundary nodes
1775 : : const auto& lid = std::get< 2 >( m_outmesh.chunk );
1776 [ + + ]: 50115 : for (const auto& [g,f] : m_uNodefieldsc) {
1777 : : Assert( m_uNodefields.nprop() == f.first.size(), "Size mismatch" );
1778 : 48837 : auto p = tk::cref_find( lid, g );
1779 [ + + ]: 311536 : for (std::size_t i=0; i<f.first.size(); ++i) {
1780 : 262699 : m_uNodefields(p,i) += f.first[i];
1781 : 262699 : m_uNodefields(p,i) /= static_cast< tk::real >(
1782 : 262699 : esup.second[p+1] - esup.second[p] + f.second );
1783 : : }
1784 : : }
1785 : 1278 : tk::destroy( m_uNodefieldsc );
1786 [ + + ]: 50115 : for (const auto& [g,f] : m_pNodefieldsc) {
1787 : : Assert( m_pNodefields.nprop() == f.first.size(), "Size mismatch" );
1788 : 48837 : auto p = tk::cref_find( lid, g );
1789 [ + + ]: 99437 : for (std::size_t i=0; i<f.first.size(); ++i) {
1790 : 50600 : m_pNodefields(p,i) += f.first[i];
1791 : 50600 : m_pNodefields(p,i) /= static_cast< tk::real >(
1792 : 50600 : esup.second[p+1] - esup.second[p] + f.second );
1793 : : }
1794 : : }
1795 : 1278 : tk::destroy( m_pNodefieldsc );
1796 : :
1797 : : // Lambda to decide if a node (global id) is on a chare boundary of the field
1798 : : // output mesh. p - global node id, return true if node is on the chare
1799 : : // boundary.
1800 : : auto chbnd = [ this ]( std::size_t p ) {
1801 : : return
1802 : : std::any_of( m_outmesh.nodeCommMap.cbegin(), m_outmesh.nodeCommMap.cend(),
1803 : : [&](const auto& s) { return s.second.find(p) != s.second.cend(); } );
1804 : : };
1805 : :
1806 : : // Finish computing node field output averages in internal nodes
1807 : 1278 : auto npoin = m_outmesh.coord[0].size();
1808 : : auto& gid = std::get< 1 >( m_outmesh.chunk );
1809 [ + + ]: 258216 : for (std::size_t p=0; p<npoin; ++p) {
1810 [ + + ]: 256938 : if (!chbnd(gid[p])) {
1811 : 208101 : auto n = static_cast< tk::real >( esup.second[p+1] - esup.second[p] );
1812 [ + + ]: 1028737 : for (std::size_t i=0; i<m_uNodefields.nprop(); ++i)
1813 : 820636 : m_uNodefields(p,i) /= n;
1814 [ + + ]: 311225 : for (std::size_t i=0; i<m_pNodefields.nprop(); ++i)
1815 : 103124 : m_pNodefields(p,i) /= n;
1816 : : }
1817 : : }
1818 : :
1819 : : // Collect field output from numerical solution requested by user
1820 : 1278 : auto elemfields = numericFieldOutput( m_uElemfields, tk::Centering::ELEM,
1821 [ + - ][ + - ]: 2556 : g_dgpde[Disc()->MeshId()].OutVarFn(), m_pElemfields );
[ + - ]
1822 : 1278 : auto nodefields = numericFieldOutput( m_uNodefields, tk::Centering::NODE,
1823 [ + - ][ + - ]: 2556 : g_dgpde[Disc()->MeshId()].OutVarFn(), m_pNodefields );
[ + - ]
1824 : :
1825 : : // Collect field output from analytical solutions (if exist)
1826 : 1278 : const auto& coord = m_outmesh.coord;
1827 [ + - ]: 1278 : auto geoElem = tk::genGeoElemTet( inpoel, coord );
1828 [ + - ]: 1278 : auto t = Disc()->T();
1829 [ + - ]: 1278 : analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::ELEM,
1830 [ + - ][ + - ]: 5112 : geoElem.extract_comp(1), geoElem.extract_comp(2), geoElem.extract_comp(3),
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
1831 : : t, elemfields );
1832 [ + - ]: 1278 : analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::NODE, coord[0],
1833 : : coord[1], coord[2], t, nodefields );
1834 : :
1835 : : // Add adaptive indicator array to element-centered field output
1836 [ + + ]: 1278 : if (g_inputdeck.get< tag::pref, tag::pref >()) {
1837 [ + - ]: 402 : std::vector< tk::real > ndof( begin(m_ndof), end(m_ndof) );
1838 [ + - ]: 402 : ndof.resize( nelem );
1839 [ + + ]: 182529 : for(std::size_t k = 0; k < nelem; k++) {
1840 : : // Mark the cell with THINC reconstruction as 0 for output
1841 [ + + ]: 182127 : if(m_interface[k] == 1) ndof[k] = 0;
1842 : : }
1843 [ + + ]: 87834 : for (const auto& [child,parent] : addedTets)
1844 : 87432 : ndof[child] = static_cast< tk::real >( m_ndof[parent] );
1845 [ + - ]: 402 : elemfields.push_back( ndof );
1846 : : }
1847 : :
1848 : : // Add shock detection marker array to element-centered field output
1849 [ + - ][ - - ]: 1278 : std::vector< tk::real > shockmarker( begin(m_shockmarker), end(m_shockmarker) );
1850 : : // Here m_shockmarker has a size of m_u.nunk() which is the number of the
1851 : : // elements within this partition (nelem) plus the ghost partition cells. In
1852 : : // terms of output purpose, we only need the solution data within this
1853 : : // partition. Therefore, resizing it to nelem removes the extra partition
1854 : : // boundary allocations in the shockmarker vector. Since the code assumes that
1855 : : // the boundary elements are on the top, the resize operation keeps the lower
1856 : : // portion.
1857 [ + - ]: 1278 : shockmarker.resize( nelem );
1858 [ + + ]: 158790 : for (const auto& [child,parent] : addedTets)
1859 : 157512 : shockmarker[child] = static_cast< tk::real >(m_shockmarker[parent]);
1860 [ + - ]: 1278 : elemfields.push_back( shockmarker );
1861 : :
1862 : : // Add rho0*det(g)/rho to make sure it is staying close to 1,
1863 : : // averaged for all materials
1864 [ + - ][ - - ]: 1278 : std::vector< tk::real > densityConstr(nelem);
1865 [ + - ]: 1278 : g_dgpde[d->MeshId()].computeDensityConstr(nelem, m_u, densityConstr);
1866 [ + + ]: 158790 : for (const auto& [child,parent] : addedTets)
1867 : 157512 : densityConstr[child] = 0.0;
1868 [ + + ][ + - ]: 1278 : if (densityConstr.size() > 0) elemfields.push_back( densityConstr );
1869 : :
1870 : : // Query fields names requested by user
1871 [ + - ]: 2556 : auto elemfieldnames = numericFieldNames( tk::Centering::ELEM );
1872 [ + - ]: 2556 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
1873 : :
1874 : : // Collect field output names for analytical solutions
1875 [ + - ]: 1278 : analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::ELEM, elemfieldnames );
1876 [ + - ]: 1278 : analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::NODE, nodefieldnames );
1877 : :
1878 [ + + ]: 1278 : if (g_inputdeck.get< tag::pref, tag::pref >()) {
1879 [ + - ]: 804 : elemfieldnames.push_back( "NDOF" );
1880 : : }
1881 : :
1882 [ + - ]: 1278 : elemfieldnames.push_back( "shock_marker" );
1883 : :
1884 [ + + ]: 1278 : if (densityConstr.size() > 0)
1885 [ + - ]: 354 : elemfieldnames.push_back( "density_constraint" );
1886 : :
1887 : : Assert( elemfieldnames.size() == elemfields.size(), "Size mismatch" );
1888 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
1889 : :
1890 : : // Collect surface output names
1891 [ + - ]: 2556 : auto surfnames = g_dgpde[d->MeshId()].surfNames();
1892 : :
1893 : : // Collect surface field solution
1894 [ + - ]: 1278 : const auto& fd = myGhosts()->m_fd;
1895 [ + - ]: 1278 : auto elemsurfs = g_dgpde[d->MeshId()].surfOutput(fd, m_u, m_p);
1896 : :
1897 : : // Output chare mesh and fields metadata to file
1898 : 1278 : const auto& triinpoel = m_outmesh.triinpoel;
1899 [ + - ][ + - ]: 3834 : d->write( inpoel, m_outmesh.coord, m_outmesh.bface, {},
[ + + ][ - - ]
1900 [ + - ]: 2556 : tk::remap( triinpoel, lid ), elemfieldnames, nodefieldnames,
1901 : : surfnames, {}, elemfields, nodefields, elemsurfs, {}, c );
1902 : 1278 : }
1903 : :
1904 : : void
1905 : 10272 : DG::comnodeout( const std::vector< std::size_t >& gid,
1906 : : const std::vector< std::size_t >& nesup,
1907 : : const std::vector< std::vector< tk::real > >& Lu,
1908 : : const std::vector< std::vector< tk::real > >& Lp )
1909 : : // *****************************************************************************
1910 : : // Receive chare-boundary nodal solution (for field output) contributions from
1911 : : // neighboring chares
1912 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
1913 : : //! \param[in] nesup Number of elements surrounding points
1914 : : //! \param[in] Lu Partial contributions of solution nodal fields to
1915 : : //! chare-boundary nodes
1916 : : //! \param[in] Lp Partial contributions of primitive quantity nodal fields to
1917 : : //! chare-boundary nodes
1918 : : // *****************************************************************************
1919 : : {
1920 : : Assert( gid.size() == nesup.size(), "Size mismatch" );
1921 : : Assert(Lu.size() == m_uNodefields.nprop(), "Fields size mismatch");
1922 : : Assert(Lp.size() == m_pNodefields.nprop(), "Fields size mismatch");
1923 [ + + ]: 68524 : for (std::size_t f=0; f<Lu.size(); ++f)
1924 : : Assert( gid.size() == Lu[f].size(), "Size mismatch" );
1925 [ + + ]: 20900 : for (std::size_t f=0; f<Lp.size(); ++f)
1926 : : Assert( gid.size() == Lp[f].size(), "Size mismatch" );
1927 : :
1928 [ + + ]: 105698 : for (std::size_t i=0; i<gid.size(); ++i) {
1929 : : auto& nfu = m_uNodefieldsc[ gid[i] ];
1930 : 95426 : nfu.first.resize( Lu.size() );
1931 [ + + ]: 631622 : for (std::size_t f=0; f<Lu.size(); ++f) nfu.first[f] += Lu[f][i];
1932 : 95426 : nfu.second += nesup[i];
1933 : 95426 : auto& nfp = m_pNodefieldsc[ gid[i] ];
1934 : 95426 : nfp.first.resize( Lp.size() );
1935 [ + + ]: 204906 : for (std::size_t f=0; f<Lp.size(); ++f) nfp.first[f] += Lp[f][i];
1936 : 95426 : nfp.second += nesup[i];
1937 : : }
1938 : :
1939 : : // When we have heard from all chares we communicate with, this chare is done
1940 [ + + ]: 10272 : if (++m_nnod == Disc()->NodeCommMap().size()) {
1941 : 1128 : m_nnod = 0;
1942 : 1128 : comnodeout_complete();
1943 : : }
1944 : 10272 : }
1945 : :
1946 : : void
1947 : 37410 : DG::stage()
1948 : : // *****************************************************************************
1949 : : // Evaluate whether to continue with next time step stage
1950 : : // *****************************************************************************
1951 : : {
1952 : : // Increment Runge-Kutta stage counter
1953 : 37410 : ++m_stage;
1954 : :
1955 : : // if not all Runge-Kutta stages complete, continue to next time stage,
1956 : : // otherwise prepare for nodal field output
1957 [ + + ]: 37410 : if (m_stage < m_nstage)
1958 : 24940 : next();
1959 : : else
1960 [ + - ][ + - ]: 37410 : startFieldOutput( CkCallback(CkIndex_DG::step(), thisProxy[thisIndex]) );
[ - + ][ - - ]
1961 : 37410 : }
1962 : :
1963 : : void
1964 : 11805 : DG::evalLB( int nrestart )
1965 : : // *****************************************************************************
1966 : : // Evaluate whether to do load balancing
1967 : : //! \param[in] nrestart Number of times restarted
1968 : : // *****************************************************************************
1969 : : {
1970 : 11805 : auto d = Disc();
1971 : :
1972 : : // Detect if just returned from a checkpoint and if so, zero timers
1973 : 11805 : d->restarted( nrestart );
1974 : :
1975 : 11805 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1976 : 11805 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
1977 : :
1978 : : // Load balancing if user frequency is reached or after the second time-step
1979 [ - + ][ - - ]: 11805 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
1980 : :
1981 : 11805 : AtSync();
1982 [ - + ]: 11805 : if (nonblocking) next();
1983 : :
1984 : : } else {
1985 : :
1986 : 0 : next();
1987 : :
1988 : : }
1989 : 11805 : }
1990 : :
1991 : : void
1992 : 11805 : DG::evalRestart()
1993 : : // *****************************************************************************
1994 : : // Evaluate whether to save checkpoint/restart
1995 : : // *****************************************************************************
1996 : : {
1997 : 11805 : auto d = Disc();
1998 : :
1999 : 11805 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
2000 : 11805 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
2001 : :
2002 [ + + ][ - + ]: 11805 : if (not benchmark and not (d->It() % rsfreq)) {
2003 : :
2004 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
2005 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
2006 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
[ - - ]
2007 : :
2008 : : } else {
2009 : :
2010 : 11805 : evalLB( /* nrestart = */ -1 );
2011 : :
2012 : : }
2013 : 11805 : }
2014 : :
2015 : : void
2016 : 12470 : DG::step()
2017 : : // *****************************************************************************
2018 : : // Evaluate wether to continue with next time step
2019 : : // *****************************************************************************
2020 : : {
2021 : 12470 : auto d = Disc();
2022 : :
2023 : : // Output time history
2024 [ + - ][ + - ]: 12470 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
2025 : 0 : std::vector< std::vector< tk::real > > hist;
2026 [ - - ][ - - ]: 0 : auto h = g_dgpde[d->MeshId()].histOutput( d->Hist(), myGhosts()->m_inpoel,
2027 [ - - ][ - - ]: 0 : myGhosts()->m_coord, m_u, m_p );
2028 [ - - ]: 0 : hist.insert( end(hist), begin(h), end(h) );
2029 [ - - ]: 0 : d->history( std::move(hist) );
2030 : : }
2031 : :
2032 : : // Free memory storing output mesh
2033 : 12470 : m_outmesh.destroy();
2034 : :
2035 : : // Output one-liner status report to screen
2036 : 12470 : d->status();
2037 : : // Reset Runge-Kutta stage counter
2038 : 12470 : m_stage = 0;
2039 : :
2040 : 12470 : const auto term = g_inputdeck.get< tag::term >();
2041 : 12470 : const auto nstep = g_inputdeck.get< tag::nstep >();
2042 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
2043 : :
2044 : : // If neither max iterations nor max time reached, continue, otherwise finish
2045 [ + - ][ + + ]: 12470 : if (std::fabs(d->T()-term) > eps && d->It() < nstep) {
2046 : :
2047 : 11805 : evalRestart();
2048 : :
2049 : : } else {
2050 : :
2051 : 665 : auto meshid = d->MeshId();
2052 [ + - ]: 1330 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2053 : 1330 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
2054 : :
2055 : : }
2056 : 12470 : }
2057 : :
2058 : : void
2059 : 0 : DG::imex_integrate()
2060 : : // *****************************************************************************
2061 : : // Perform the Implicit-Explicit Runge-Kutta stage update
2062 : : //
2063 : : //! \details
2064 : : //! Performs the Implicit-Explicit Runge-Kutta step. Scheme taken from
2065 : : //! Cavaglieri, D., & Bewley, T. (2015). Low-storage implicit/explicit
2066 : : //! Runge–Kutta schemes for the simulation of stiff high-dimensional ODE
2067 : : //! systems. Journal of Computational Physics, 286, 172-193.
2068 : : //!
2069 : : //! Scheme given by equations (25a,b):
2070 : : //!
2071 : : //! u[0] = u[n] + dt * (expl_rkcoef[1,0]*R_ex(u[n])+impl_rkcoef[1,1]*R_im(u[0]))
2072 : : //!
2073 : : //! u[1] = u[n] + dt * (expl_rkcoef[2,1]*R_ex(u[0])+impl_rkcoef[2,1]*R_im(u[0])
2074 : : //! +impl_rkcoef[2,2]*R_im(u[1]))
2075 : : //!
2076 : : //! u[n+1] = u[n] + dt * (expl_rkcoef[3,1]*R_ex(u[0])+impl_rkcoef[3,1]*R_im(u[0])
2077 : : //! expl_rkcoef[3,2]*R_ex(u[1])+impl_rkcoef[3,2]*R_im(u[1]))
2078 : : //!
2079 : : //! In order to solve the first two equations we need to solve a series of
2080 : : //! systems of non-linear equations:
2081 : : //!
2082 : : //! F1(u[0]) = B1 + R1(u[0]) = 0, and
2083 : : //! F2(u[1]) = B2 + R2(u[1]) = 0,
2084 : : //!
2085 : : //! where
2086 : : //!
2087 : : //! B1 = u[n] + dt * expl_rkcoef[1,0]*R_ex(u[n]),
2088 : : //! R1 = dt * impl_rkcoef[1,1]*R_im(u[0]) - u([0]),
2089 : : //! B2 = u[n] + dt * (expl_rkcoef[2,1]*R_ex(u[0])+impl_rkcoef[2,1]*R_im(u[0])),
2090 : : //! R2 = dt * impl_rkcoef[2,2]*R_im(u[1]) - u([1]).
2091 : : //!
2092 : : //! In order to solve the non-linear system F(U) = 0, we employ:
2093 : : //! First, Broyden's method.
2094 : : //! If that fails, Newton's method (with FD approximation for jacobian).
2095 : : //!
2096 : : //! Broyden's method:
2097 : : //! ----------------
2098 : : //!
2099 : : //! Taken from https://en.wikipedia.org/wiki/Broyden%27s_method.
2100 : : //! The method consists in obtaining an approximation for the inverse of the
2101 : : //! Jacobian H = J^(-1) and advancing in a quasi-newton step:
2102 : : //!
2103 : : //! U[k+1] = U[k] - H[k]*F(U[k]),
2104 : : //!
2105 : : //! until F(U) is close enough to zero.
2106 : : //!
2107 : : //! The approximation H[k] is improved at every iteration following
2108 : : //!
2109 : : //! 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],
2110 : : //!
2111 : : //! where DU[k] = U[k] - U[k-1] and DF[k] = F(U[k]) - F(U[k-1)).
2112 : : //!
2113 : : //! Newton's method:
2114 : : //! ----------------
2115 : : //!
2116 : : //! Taken from https://en.wikipedia.org/wiki/Newton%27s_method
2117 : : //! The method consists in inverting the jacobian
2118 : : //! Jacobian J and advancing in a Newton step:
2119 : : //!
2120 : : //! U[k+1] = U[k] - J^(-1)[k]*F(U[k]),
2121 : : //!
2122 : : //! until F(U) is close enough to zero.
2123 : : //!
2124 : : //!
2125 : : //! This function performs the following main algorithmic steps:
2126 : : //! - If stage == 0 or stage == 1:
2127 : : //! - Take Initial value:
2128 : : //! U[0] = U[n] + dt * expl_rkcoef[1,0]*R_ex(U[n]) (for stage 0)
2129 : : //! U[1] = U[n] + dt * (expl_rkcoef[2,1]*R_ex(U[0])
2130 : : //! +impl_rkcoef[2,1]*R_im(U[0])) (for stage 1)
2131 : : //! - Loop over the Elements (e++)
2132 : : //! - Broyden steps:
2133 : : //! - Initialize Jacobian inverse approximation using FD
2134 : : //! - Compute implicit right-hand-side (F_im) with current U
2135 : : //! - Iterate for the solution (iter++)
2136 : : //! - Perform line search prior to solution update
2137 : : //! - Compute new solution U[k+1] = U[k] - H[k]*F(U[k])
2138 : : //! - Compute implicit right-hand-side (F_im) with current U
2139 : : //! - Compute DU and DF
2140 : : //! - Update inverse Jacobian approximation by:
2141 : : //! - Compute V1 = H[k-1]*DF[k] and V2 = DU[k]^T*H[k-1]
2142 : : //! - Compute d = DU[k]^T*V1 and V3 = DU[k]-V1
2143 : : //! - Compute V4 = V3/d
2144 : : //! - Update H[k] = H[k-1] + V4*V2
2145 : : //! - Save old U and F
2146 : : //! - Compute absolute and relative errors
2147 : : //! - Break iterations if error < tol or iter == max_iter
2148 : : //! - Newton steps:
2149 : : //! - Initialize Jacobian using FD approximation.
2150 : : //! - Compute implicit right-hand-side (F_im) with current U
2151 : : //! - Iterate for the solution (iter++)
2152 : : //! - Perform line search prior to solution update
2153 : : //! - Compute new solution U[k+1] = U[k] - J^(-1)[k]*F(U[k])
2154 : : //! - Compute implicit right-hand-side (F_im) with current U
2155 : : //! - Compute DU and DF
2156 : : //! - Save old U and F
2157 : : //! - Compute absolute and relative errors
2158 : : //! - Break iterations if error < tol or iter == max_iter
2159 : : //! - Update explicit equations using only the explicit terms.
2160 : : //! - Else if stage == 2:
2161 : : //! - Update explicit equations using only the explicit terms.
2162 : : //! - Update implicit equations using:
2163 : : //! u[n+1] = u[n]+dt*(expl_rkcoef[3,1]*R_ex(u[0])+impl_rkcoef[3,1]*R_im(u[0])
2164 : : //! expl_rkcoef[3,2]*R_ex(u[1])+impl_rkcoef[3,2]*R_im(u[1]))
2165 : : // *****************************************************************************
2166 : : {
2167 : 0 : auto d = Disc();
2168 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
2169 : 0 : const auto ndof = g_inputdeck.get< tag::ndof >();
2170 [ - - ]: 0 : if (m_stage < m_nstage-1) {
2171 : : // Save previous stiff_rhs
2172 : : m_stiffrhsprev = m_stiffrhs;
2173 : :
2174 : : // Compute the imex update
2175 : :
2176 : : // Integrate explicitly on the imex equations
2177 : : // (To use as initial values)
2178 [ - - ]: 0 : for (std::size_t e=0; e<myGhosts()->m_nunk; ++e) {
2179 : 0 : auto vole = myGhosts()->m_geoElem(e,0);
2180 [ - - ]: 0 : for (std::size_t c=0; c<m_nstiffeq; ++c)
2181 : : {
2182 [ - - ]: 0 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
2183 : : {
2184 [ - - ]: 0 : auto rmark = m_stiffEqIdx[c]*rdof+k;
2185 [ - - ]: 0 : auto mark = m_stiffEqIdx[c]*ndof+k;
2186 : 0 : m_u(e, rmark) = m_un(e, rmark) + d->Dt() * (
2187 : 0 : expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/(vole*mass_dubiner[k])
2188 : 0 : + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/(vole*mass_dubiner[k])
2189 : 0 : + impl_rkcoef[0][m_stage]
2190 : 0 : * m_stiffrhsprev(e,c*ndof+k)/(vole*mass_dubiner[k]) );
2191 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2192 : 0 : m_u(e, rmark) = 0;
2193 : : }
2194 : : }
2195 : : }
2196 : :
2197 : : // Solve for implicit-explicit equations
2198 : 0 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
2199 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
2200 : : {
2201 : :
2202 : : // Non-linear solver solves for x.
2203 : : // Copy the relevant variables from the state array into x.
2204 : 0 : std::vector< tk::real > x(m_nstiffeq*ndof, 0.0);
2205 [ - - ]: 0 : for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2206 [ - - ]: 0 : for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2207 : : {
2208 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2209 : 0 : x[ieq*ndof+idof] = m_u(e, stiffrmark);
2210 : : }
2211 : :
2212 : : // Save all the values of m_u at stiffEqIdx as x_star,
2213 : : // They will serve to balance the energy exchange
2214 : : // from the implicit step
2215 [ - - ]: 0 : auto x_star = x;
2216 : :
2217 : : // Solve nonlinear system, first try broyden
2218 : : bool solver_failed = false;
2219 [ - - ][ - - ]: 0 : x = DG::nonlinear_broyden(e, x, solver_failed);
[ - - ]
2220 : :
2221 : : // If solver_failed, do newton
2222 : : if (solver_failed) {
2223 : : solver_failed = false;
2224 : : x = DG::nonlinear_newton(e, x, solver_failed);
2225 : : }
2226 : :
2227 : : // If newton failed, crash
2228 : : if (solver_failed)
2229 : : Throw("At element " + std::to_string(e) +
2230 : : " nonlinear solvers was not able to converge");
2231 : :
2232 : : // Balance energy
2233 [ - - ][ - - ]: 0 : g_dgpde[d->MeshId()].balance_plastic_energy(e, x_star, x, m_un);
[ - - ][ - - ]
[ - - ][ - - ]
2234 : :
2235 : : // Update the state u with the converged vector x.
2236 [ - - ]: 0 : for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2237 [ - - ]: 0 : for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2238 : : {
2239 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2240 : 0 : m_u(e, stiffrmark) = x[ieq*ndof+idof];
2241 : : }
2242 : :
2243 : : }
2244 : :
2245 : : // Then, integrate explicitly on the remaining equations
2246 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e) {
2247 : 0 : auto vole = myGhosts()->m_geoElem(e,0);
2248 [ - - ]: 0 : for (std::size_t c=0; c<m_nnonstiffeq; ++c)
2249 : : {
2250 [ - - ]: 0 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
2251 : : {
2252 [ - - ]: 0 : auto rmark = m_nonStiffEqIdx[c]*rdof+k;
2253 [ - - ]: 0 : auto mark = m_nonStiffEqIdx[c]*ndof+k;
2254 : 0 : m_u(e, rmark) = m_un(e, rmark) + d->Dt() * (
2255 : 0 : expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/(vole*mass_dubiner[k])
2256 : 0 : + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/(vole*mass_dubiner[k]));
2257 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2258 : 0 : m_u(e, rmark) = 0;
2259 : : }
2260 : : }
2261 : : }
2262 : : }
2263 : : else {
2264 : : // For last stage just use all previously computed stages
2265 : 0 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
2266 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
2267 : : {
2268 : 0 : auto vole = myGhosts()->m_geoElem(e,0);
2269 : : // First integrate explicitly on nonstiff equations
2270 [ - - ]: 0 : for (std::size_t c=0; c<m_nnonstiffeq; ++c)
2271 : : {
2272 [ - - ]: 0 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
2273 : : {
2274 [ - - ]: 0 : auto rmark = m_nonStiffEqIdx[c]*rdof+k;
2275 [ - - ]: 0 : auto mark = m_nonStiffEqIdx[c]*ndof+k;
2276 : 0 : m_u(e, rmark) = m_un(e, rmark) + d->Dt() * (
2277 : 0 : expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/(vole*mass_dubiner[k])
2278 : 0 : + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/(vole*mass_dubiner[k]));
2279 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2280 : 0 : m_u(e, rmark) = 0;
2281 : : }
2282 : : }
2283 : : // Then, integrate the imex-equations
2284 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2285 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2286 : : {
2287 [ - - ]: 0 : auto rmark = m_stiffEqIdx[ieq]*rdof+idof;
2288 [ - - ]: 0 : auto mark = m_stiffEqIdx[ieq]*ndof+idof;
2289 : 0 : m_u(e, rmark) = m_un(e, rmark)
2290 : 0 : + d->Dt() * (expl_rkcoef[0][m_stage]
2291 : 0 : * m_rhsprev(e,mark)/(vole*mass_dubiner[idof])
2292 : 0 : + expl_rkcoef[1][m_stage]
2293 : 0 : * m_rhs(e,mark)/(vole*mass_dubiner[idof])
2294 : 0 : + impl_rkcoef[0][m_stage]
2295 : 0 : * m_stiffrhsprev(e,ieq*ndof+idof)/(vole*mass_dubiner[idof])
2296 : 0 : + impl_rkcoef[1][m_stage]
2297 : 0 : * m_stiffrhs(e,ieq*ndof+idof)/(vole*mass_dubiner[idof]) );
2298 [ - - ]: 0 : if(fabs(m_u(e, rmark)) < 1e-16)
2299 : 0 : m_u(e, rmark) = 0;
2300 : : }
2301 : : }
2302 : : }
2303 : 0 : }
2304 : :
2305 : : void
2306 : 0 : DG::BDF1_integrate()
2307 : : // *****************************************************************************
2308 : : // Perform the BDF1 update
2309 : : //! \details This function updates the solution using the BDF1 (backward Euler)
2310 : : //! time discretization.
2311 : : // *****************************************************************************
2312 : : {
2313 : : //TODO: implicit solver:
2314 : : // update solution m_u
2315 : 0 : }
2316 : :
2317 : 0 : std::vector< tk::real > DG::nonlinear_func(std::size_t e,
2318 : : std::vector< tk::real > x)
2319 : : // *****************************************************************************
2320 : : // Evaluate the stiff RHS and stiff equations f = b - A(x)
2321 : : //! \param[in] e Element number
2322 : : //! \param[in,out] x Array of unknowns to solve for
2323 : : //! \details
2324 : : //! Defines the F(x) function that the non-linear solvers
2325 : : //! look to minimize. Deals with properly calling the stiff
2326 : : //! RHS functions.
2327 : : // *****************************************************************************
2328 : : {
2329 : 0 : auto d = Disc();
2330 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
2331 : 0 : const auto ndof = g_inputdeck.get< tag::ndof >();
2332 : 0 : std::size_t n = x.size();
2333 : :
2334 : : // m_u <- x
2335 [ - - ]: 0 : for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2336 [ - - ]: 0 : for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2337 : : {
2338 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2339 : 0 : m_u(e, stiffrmark) = x[ieq*ndof+idof];
2340 : : }
2341 : :
2342 : 0 : auto vole = myGhosts()->m_geoElem(e,0);
2343 : :
2344 : : // Compute explicit terms (Should be computed once)
2345 : 0 : std::vector< tk::real > expl_terms(n, 0.0);
2346 [ - - ]: 0 : for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2347 [ - - ]: 0 : for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2348 : : {
2349 : 0 : auto stiffmark = m_stiffEqIdx[ieq]*ndof+idof;
2350 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2351 : 0 : expl_terms[ieq*ndof+idof] = m_un(e, stiffrmark)
2352 : 0 : + d->Dt() * ( expl_rkcoef[0][m_stage]
2353 : 0 : * m_rhsprev(e,stiffmark)/(vole*mass_dubiner[idof])
2354 : 0 : + expl_rkcoef[1][m_stage]
2355 : 0 : * m_rhs(e,stiffmark)/(vole*mass_dubiner[idof])
2356 : 0 : + impl_rkcoef[0][m_stage]
2357 : 0 : * m_stiffrhsprev(e,ieq*ndof+idof)/(vole*mass_dubiner[idof]) );
2358 : : }
2359 : :
2360 : : // Compute stiff_rhs
2361 [ - - ][ - - ]: 0 : g_dgpde[d->MeshId()].stiff_rhs( e, myGhosts()->m_geoElem,
2362 [ - - ][ - - ]: 0 : myGhosts()->m_inpoel, myGhosts()->m_coord,
2363 [ - - ]: 0 : m_u, m_p, m_ndof, m_stiffrhs );
2364 : :
2365 : : // Store f
2366 [ - - ][ - - ]: 0 : std::vector< tk::real > f(n, 0.0);
2367 [ - - ]: 0 : for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
2368 [ - - ]: 0 : for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
2369 : : {
2370 : 0 : auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
2371 : 0 : f[ieq*ndof+idof] = expl_terms[ieq*ndof+idof]
2372 : 0 : + d->Dt() * impl_rkcoef[1][m_stage]
2373 : 0 : * m_stiffrhs(e,ieq*ndof+idof)/(vole*mass_dubiner[idof])
2374 : 0 : - m_u(e, stiffrmark);
2375 : : }
2376 : :
2377 : 0 : return f;
2378 : : }
2379 : :
2380 : 0 : std::vector< tk::real > DG::nonlinear_broyden(std::size_t e,
2381 : : std::vector< tk::real > x,
2382 : : bool solver_failed )
2383 : : // *****************************************************************************
2384 : : // Performs Broyden's method to solve a non-linear system on
2385 : : // element e.
2386 : : //! \param[in] e Element number
2387 : : //! \param[in,out] x Array of unknowns to solve for
2388 : : //! \param[out] solver_failed Returns true if solver did not converge
2389 : : //! \details
2390 : : //! Taken from https://en.wikipedia.org/wiki/Broyden%27s_method.
2391 : : //! The method consists in obtaining an approximation for the inverse of the
2392 : : //! Jacobian H = J^(-1) and advancing in a quasi-newton step:
2393 : : //!
2394 : : //! U[k+1] = U[k] - H[k]*F(U[k]),
2395 : : //!
2396 : : //! until F(U) is close enough to zero.
2397 : : // *****************************************************************************
2398 : : {
2399 : : // Broyden's method
2400 : : // Control parameters
2401 : 0 : std::size_t max_iter = g_inputdeck.get< tag::imex_maxiter >();
2402 : 0 : tk::real rel_tol = g_inputdeck.get< tag::imex_reltol >();
2403 : 0 : tk::real abs_tol = g_inputdeck.get< tag::imex_abstol >();
2404 : : tk::real rel_err = rel_tol+1;
2405 : : tk::real abs_err = abs_tol+1;
2406 : 0 : std::size_t n = x.size();
2407 : :
2408 : : // Compute f with initial guess
2409 [ - - ]: 0 : std::vector< tk::real > f = DG::nonlinear_func(e, x);
2410 : :
2411 : : // Initialize x_old and f_old
2412 [ - - ][ - - ]: 0 : std::vector< tk::real > x_old(n, 0.0), f_old(n, 0.0);
[ - - ][ - - ]
2413 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2414 : : {
2415 : 0 : x_old[i] = x[i];
2416 : 0 : f_old[i] = f[i];
2417 : : }
2418 : :
2419 : : // Initialize delta_x and delta_f
2420 [ - - ][ - - ]: 0 : std::vector< tk::real > delta_x(n, 0.0), delta_f(n, 0.0);
[ - - ][ - - ]
2421 : :
2422 : : // Store the norm of f initially, for relative error measure
2423 : : tk::real err0 = 0.0;
2424 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2425 : 0 : err0 += f[i]*f[i];
2426 : 0 : err0 = std::sqrt(err0);
2427 : :
2428 : : // Iterate for the solution if err0 > 0
2429 : : solver_failed = false;
2430 : : tk::real alpha_jacob = 1.0;
2431 [ - - ]: 0 : if (err0 > abs_tol) {
2432 : :
2433 : : // Evaluate finite difference based jacobian
2434 [ - - ][ - - ]: 0 : std::vector< double > jacob(n*n);
2435 : : tk::real dx = 0.0;
2436 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2437 [ - - ]: 0 : for (std::size_t j=0; j<n; ++j)
2438 : : {
2439 : : // Set dx in the order 1% of the unknown
2440 [ - - ]: 0 : dx = std::max(std::abs(0.1*x[j]),1.0e-06);
2441 : : // Derivative of f[i] with respect to x[j]
2442 [ - - ]: 0 : auto x_perturb = x;
2443 [ - - ]: 0 : x_perturb[j] += dx;
2444 [ - - ][ - - ]: 0 : auto f_perturb = DG::nonlinear_func(e, x_perturb);
[ - - ]
2445 [ - - ]: 0 : jacob[i*n+j] = (f_perturb[i]-f[i])/dx;
2446 : : }
2447 : :
2448 : : // Initialize Jacobian to be the inverse of this jacobian
2449 : : lapack_int ln = static_cast< lapack_int >(n);
2450 [ - - ][ - - ]: 0 : std::vector< lapack_int > ipiv(n);
2451 : :
2452 : : #ifndef NDEBUG
2453 : : lapack_int ierr =
2454 : : #endif
2455 [ - - ]: 0 : LAPACKE_dgetrf(LAPACK_ROW_MAJOR, ln, ln, jacob.data(), ln, ipiv.data());
2456 : : Assert(ierr==0, "Lapack error in LU factorization of FD Jacobian");
2457 : :
2458 : : #ifndef NDEBUG
2459 : : lapack_int jerr =
2460 : : #endif
2461 [ - - ]: 0 : LAPACKE_dgetri(LAPACK_ROW_MAJOR, ln, jacob.data(), ln, ipiv.data());
2462 : : Assert(jerr==0, "Lapack error in inverting FD Jacobian");
2463 : :
2464 : : std::vector< std::vector< tk::real > >
2465 [ - - ][ - - ]: 0 : approx_jacob(n, std::vector< tk::real >(n, 0.0));
[ - - ]
2466 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2467 [ - - ]: 0 : for (std::size_t j=0; j<n; ++j)
2468 : 0 : approx_jacob[i][j] = jacob[i*n+j];
2469 : :
2470 [ - - ]: 0 : for (size_t iter=0; iter<max_iter; ++iter)
2471 : : {
2472 : :
2473 : : // Scale inverse of jacobian if things are not going well
2474 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2475 [ - - ]: 0 : for (std::size_t j=0; j<n; ++j)
2476 : 0 : approx_jacob[i][j] *= alpha_jacob;
2477 : :
2478 : : // Compute new solution
2479 [ - - ]: 0 : std::vector < tk::real > delta(n, 0.0);
2480 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2481 : : {
2482 : 0 : delta[i] = 0.0;
2483 [ - - ]: 0 : for (std::size_t j=0; j<n; ++j)
2484 : 0 : delta[i] -= approx_jacob[i][j] * f[j];
2485 : : }
2486 : :
2487 : : // Update x using line search
2488 : : bool ls_failed = false;
2489 : : tk::real alpha_ls = 1.0E+00;
2490 : : std::size_t nline = 25;
2491 [ - - ]: 0 : auto xtest = x;
2492 [ - - ]: 0 : for (std::size_t iline = 0; iline<nline; ++iline)
2493 : : {
2494 : : // Evaluate xtest
2495 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2496 : 0 : xtest[i] = x[i] + alpha_ls*delta[i];
2497 : :
2498 : : // Compute new f(x)
2499 [ - - ][ - - ]: 0 : f = DG::nonlinear_func(e, xtest);
[ - - ]
2500 : :
2501 : : tk::real err = 0.0;
2502 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2503 : 0 : err += f[i]*f[i];
2504 : 0 : abs_err = std::sqrt(err);
2505 : :
2506 : : // If 1. The error went up
2507 : : // or 2. The function f flipped in sign
2508 : : // Reduce the factor alpha_ls
2509 : : bool flipped_sign = false;
2510 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2511 [ - - ]: 0 : if (f_old[i]*f[i] < 0.0) {
2512 : : flipped_sign = true;
2513 : : break;
2514 : : }
2515 : :
2516 [ - - ]: 0 : if (!flipped_sign)
2517 : : {
2518 : : break;
2519 : : }
2520 : : else
2521 : : {
2522 : 0 : alpha_ls *= 0.5;
2523 : : }
2524 [ - - ]: 0 : if (iline == nline-1) {
2525 : : // Try again by reducing the jacobian,
2526 : : // but only a few times, otherwise give up
2527 : 0 : alpha_jacob *= 0.5;
2528 [ - - ]: 0 : if (alpha_jacob < 1.0E-03)
2529 : : solver_failed = true;
2530 : : else
2531 : : ls_failed = true;
2532 : : }
2533 : : }
2534 : :
2535 [ - - ]: 0 : if (solver_failed) {
2536 : : break;
2537 : : }
2538 : :
2539 [ - - ]: 0 : if (!ls_failed) {
2540 : : // Save x
2541 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2542 : 0 : x[i] = xtest[i];
2543 : :
2544 : : // Compute delta_x and delta_f
2545 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2546 : : {
2547 : 0 : delta_x[i] = x[i] - x_old[i];
2548 : 0 : delta_f[i] = f[i] - f_old[i];
2549 : : }
2550 : :
2551 : : // Update inverse Jacobian approximation
2552 : :
2553 : : // 1. Compute approx_jacob*delta_f and delta_x*jacob_approx
2554 : : tk::real sum1, sum2;
2555 [ - - ][ - - ]: 0 : std::vector< tk::real > auxvec1(n, 0.0), auxvec2(n, 0.0);
[ - - ][ - - ]
2556 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2557 : : {
2558 : : sum1 = 0.0;
2559 : : sum2 = 0.0;
2560 [ - - ]: 0 : for (std::size_t j=0; j<n; ++j)
2561 : : {
2562 : 0 : sum1 += approx_jacob[i][j]*delta_f[j];
2563 : 0 : sum2 += delta_x[j]*approx_jacob[j][i];
2564 : : }
2565 : 0 : auxvec1[i] = sum1;
2566 : 0 : auxvec2[i] = sum2;
2567 : : }
2568 : :
2569 : : // 2. Compute delta_x*approx_jacob*delta_f
2570 : : // and delta_x-approx_jacob*delta_f
2571 : : tk::real denom = 0.0;
2572 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2573 : : {
2574 : 0 : denom += delta_x[i]*auxvec1[i];
2575 : 0 : auxvec1[i] = delta_x[i]-auxvec1[i];
2576 : : }
2577 : :
2578 : : // 3. Divide delta_u+approx_jacob*delta_f
2579 : : // by delta_x*(approx_jacob*delta_f)
2580 [ - - ]: 0 : if (std::abs(denom) < 1.0e-18)
2581 : : {
2582 [ - - ]: 0 : if (denom < 0.0)
2583 : : {
2584 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2585 : 0 : auxvec1[i] /= -1.0e-18;
2586 : : }
2587 : : else
2588 : : {
2589 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2590 : 0 : auxvec1[i] /= 1.0e-18;
2591 : : }
2592 : : }
2593 : : else
2594 : : {
2595 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2596 : 0 : auxvec1[i] /= denom;
2597 : : }
2598 : :
2599 : : // 4. Perform outter product between the two arrays and
2600 : : // add that quantity to the new jacobian approximation
2601 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2602 [ - - ]: 0 : for (std::size_t j=0; j<n; ++j)
2603 : 0 : approx_jacob[i][j] += auxvec1[i] * auxvec2[j];
2604 : :
2605 : : // Compute a measure of error, use norm of f
2606 : : tk::real err = 0.0;
2607 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2608 : 0 : err += f[i]*f[i];
2609 : 0 : abs_err = std::sqrt(err);
2610 : 0 : rel_err = abs_err/err0;
2611 : :
2612 : : // Save solution and f
2613 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2614 : : {
2615 : 0 : x_old[i] = x[i];
2616 : 0 : f_old[i] = f[i];
2617 : : }
2618 : :
2619 : : // check if error condition is met and loop back
2620 [ - - ][ - - ]: 0 : if (rel_err < rel_tol || abs_err < abs_tol)
2621 : : break;
2622 : :
2623 : : // If we did not converge, print a message and keep going
2624 [ - - ]: 0 : if (iter == max_iter-1)
2625 : : {
2626 : : solver_failed = true;
2627 : : }
2628 : : }
2629 : : }
2630 : : }
2631 : :
2632 : 0 : return x;
2633 : : }
2634 : :
2635 : 0 : std::vector< tk::real > DG::nonlinear_newton(std::size_t e,
2636 : : std::vector< tk::real > x,
2637 : : bool solver_failed )
2638 : : // *****************************************************************************
2639 : : // Performs Newton's method to solve a non-linear system on
2640 : : // element e.
2641 : : //! \param[in] e Element number
2642 : : //! \param[in,out] x Array of unknowns to solve for
2643 : : //! \param[out] solver_failed Returns true if solver did not converge
2644 : : //! \details
2645 : : //! Taken from https://en.wikipedia.org/wiki/Newton%27s_method
2646 : : //! The method consists in inverting the jacobian
2647 : : //! Jacobian J and advancing in a Newton step:
2648 : : //!
2649 : : //! U[k+1] = U[k] - J^(-1)[k]*F(U[k]),
2650 : : //!
2651 : : //! until F(U) is close enough to zero.
2652 : : // *****************************************************************************
2653 : : {
2654 : : // Newton's method
2655 : : // Control parameters
2656 : 0 : std::size_t max_iter = g_inputdeck.get< tag::imex_maxiter >();
2657 : 0 : tk::real rel_tol = g_inputdeck.get< tag::imex_reltol >();
2658 : 0 : tk::real abs_tol = g_inputdeck.get< tag::imex_abstol >();
2659 : 0 : tk::real rel_err = rel_tol+1;
2660 : 0 : tk::real abs_err = abs_tol+1;
2661 : 0 : std::size_t n = x.size();
2662 : :
2663 : : // Define jacobian
2664 : 0 : std::vector< double > jacob(n*n);
2665 : :
2666 : : // Compute f with initial guess
2667 [ - - ][ - - ]: 0 : std::vector< tk::real > f = DG::nonlinear_func(e, x);
[ - - ]
2668 : :
2669 : : // Store the norm of f initially, for relative error measure
2670 : : tk::real err0 = 0.0;
2671 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2672 : 0 : err0 += f[i]*f[i];
2673 : 0 : err0 = std::sqrt(err0);
2674 : : auto abs_err_old = err0;
2675 : :
2676 : : // Iterate for the solution if err0 > 0
2677 : : solver_failed = false;
2678 : : tk::real alpha_jacob = 1.0;
2679 [ - - ]: 0 : if (err0 > abs_tol)
2680 [ - - ]: 0 : for (std::size_t iter=0; iter<max_iter; ++iter)
2681 : : {
2682 : :
2683 : : // Evaluate jacobian
2684 : : tk::real dx = 0.0;
2685 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2686 [ - - ]: 0 : for (std::size_t j=0; j<n; ++j)
2687 : : {
2688 : : // Set dx in the order 1% of the unknown
2689 [ - - ]: 0 : dx = alpha_jacob*std::max(std::abs(0.1*x[j]),1.0e-06);
2690 : : // Derivative of f[i] with respect to x[j]
2691 [ - - ]: 0 : auto x_perturb = x;
2692 [ - - ]: 0 : x_perturb[j] += dx;
2693 [ - - ][ - - ]: 0 : auto f_perturb = DG::nonlinear_func(e, x_perturb);
[ - - ]
2694 [ - - ]: 0 : jacob[i*n+j] = (f_perturb[i]-f[i])/dx;
2695 : : }
2696 : :
2697 : : // Compute new solution by solving linear system J*dx = -f
2698 : : lapack_int ln = static_cast< lapack_int >(n);
2699 [ - - ][ - - ]: 0 : std::vector< double > delta(n);
2700 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2701 : 0 : delta[i] = -f[i];
2702 : : lapack_int info;
2703 [ - - ][ - - ]: 0 : std::vector< lapack_int > ipiv(n);
2704 [ - - ]: 0 : info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, ln, 1, jacob.data(), ln, ipiv.data(), delta.data(), 1);
2705 : :
2706 [ - - ]: 0 : if (info != 0) {
2707 [ - - ]: 0 : printf("Failed with info: %ld\n", info);
2708 : : }
2709 : :
2710 : : // Save f as fold
2711 [ - - ][ - - ]: 0 : std::vector< tk::real > fold(n);
2712 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2713 : 0 : fold[i] = f[i];
2714 : :
2715 : : // Update x using line search
2716 : : bool ls_failed = false;
2717 : : tk::real alpha_ls = 1.0E+00;
2718 : : std::size_t nline = 25;
2719 [ - - ]: 0 : auto xtest = x;
2720 [ - - ]: 0 : for (std::size_t iline = 0; iline<nline; ++iline)
2721 : : {
2722 : : // Evaluate xtest
2723 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2724 : 0 : xtest[i] = x[i] + alpha_ls*delta[i];
2725 : :
2726 : : // Compute new f(x)
2727 [ - - ][ - - ]: 0 : f = DG::nonlinear_func(e, xtest);
[ - - ]
2728 : :
2729 : : tk::real err = 0.0;
2730 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2731 : 0 : err += f[i]*f[i];
2732 : 0 : abs_err = std::sqrt(err);
2733 : :
2734 : : // If 1. The error went up
2735 : : // or 2. The function f flipped in sign
2736 : : // Reduce the factor alpha_ls
2737 : : bool flipped_sign = false;
2738 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2739 [ - - ]: 0 : if (fold[i]*f[i] < 0.0) {
2740 : : flipped_sign = true;
2741 : : break;
2742 : : }
2743 : :
2744 [ - - ][ - - ]: 0 : if (abs_err < abs_err_old && !flipped_sign)
2745 : : {
2746 : : break;
2747 : : }
2748 : : else
2749 : : {
2750 : 0 : alpha_ls *= 0.5;
2751 : : }
2752 [ - - ]: 0 : if (iline == nline-1) {
2753 : : //printf("Line search failed to decrease f\n");
2754 : : // Try again by reducing the jacobian,
2755 : : // but only a few times, otherwise give up
2756 : 0 : alpha_jacob *= 0.5;
2757 [ - - ]: 0 : if (alpha_jacob < 1.0E-03)
2758 : : solver_failed = true;
2759 : : else
2760 : : ls_failed = true;
2761 : : }
2762 : : }
2763 : :
2764 [ - - ]: 0 : if (solver_failed) {
2765 [ - - ][ - - ]: 0 : f = DG::nonlinear_func(e, x);
[ - - ][ - - ]
2766 [ - - ]: 0 : printf("\nIMEX-RK: Non-linear solver did not converge in %lu iterations\n", iter+1);
2767 [ - - ]: 0 : printf("Element #%lu\n", e);
2768 [ - - ]: 0 : printf("Relative error: %e\n", rel_err);
2769 [ - - ]: 0 : printf("Absolute error: %e\n\n", abs_err);
2770 : : break;
2771 : : }
2772 : :
2773 [ - - ]: 0 : if (!ls_failed) {
2774 : : // Save x
2775 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2776 : 0 : x[i] = xtest[i];
2777 : :
2778 : : // Compute a measure of error, use norm of f
2779 : : tk::real err = 0.0;
2780 [ - - ]: 0 : for (std::size_t i=0; i<n; ++i)
2781 : 0 : err += f[i]*f[i];
2782 : 0 : abs_err = std::sqrt(err);
2783 : 0 : rel_err = abs_err/err0;
2784 : :
2785 : : // check if error condition is met and loop back
2786 [ - - ][ - - ]: 0 : if (rel_err < rel_tol || abs_err < abs_tol) {
2787 : : break;
2788 : : }
2789 : :
2790 : : // If we did not converge, print a message and keep going
2791 [ - - ]: 0 : if (iter == max_iter-1)
2792 : : {
2793 [ - - ]: 0 : printf("\nIMEX-RK: Non-linear solver did not converge in %lu iterations\n", max_iter);
2794 [ - - ]: 0 : printf("Element #%lu\n", e);
2795 [ - - ]: 0 : printf("Relative error: %e\n", rel_err);
2796 [ - - ]: 0 : printf("Absolute error: %e\n\n", abs_err);
2797 : : }
2798 : : }
2799 : : }
2800 : :
2801 : 0 : return x;
2802 : :
2803 : : }
2804 : :
2805 : : #include "NoWarning/dg.def.h"
|
