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