Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/DGPDE.hpp
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 Partial differential equation base for discontinuous Galerkin PDEs
9 : : \details This file defines a generic partial differential equation (PDE)
10 : : class for PDEs that use discontinuous Galerkin spatial discretization.
11 : : The class uses runtime polymorphism without client-side inheritance:
12 : : inheritance is confined to the internals of the class, invisible to
13 : : client-code. The class exclusively deals with ownership enabling client-side
14 : : value semantics. Credit goes to Sean Parent at Adobe:
15 : : https://github.com/sean-parent/sean-parent.github.com/wiki/
16 : : Papers-and-Presentations.
17 : : */
18 : : // *****************************************************************************
19 : : #ifndef DGPDE_h
20 : : #define DGPDE_h
21 : :
22 : : #include <array>
23 : : #include <string>
24 : : #include <vector>
25 : : #include <memory>
26 : : #include <unordered_set>
27 : : #include <functional>
28 : :
29 : : #include "Types.hpp"
30 : : #include "Fields.hpp"
31 : : #include "FaceData.hpp"
32 : : #include "UnsMesh.hpp"
33 : : #include "Inciter/InputDeck/InputDeck.hpp"
34 : : #include "FunctionPrototypes.hpp"
35 : : #include "History.hpp"
36 : :
37 : : namespace inciter {
38 : :
39 : : extern ctr::InputDeck g_inputdeck;
40 : :
41 : : using ncomp_t = tk::ncomp_t;
42 : : using BCStateFn =
43 : : std::vector< std::tuple< std::vector< std::size_t >, tk::StateFn, tk::StateFn > >;
44 : :
45 : : //! Extract BC configuration ignoring if BC not specified
46 : : //! \note A more preferable way of catching errors such as this function
47 : : //! hides is during parsing, so that we don't even get here if BCs are
48 : : //! not correctly specified. For now we simply ignore if BCs are not
49 : : //! specified by allowing empty BC vectors from the user input.
50 : : struct ConfigBC {
51 : : BCStateFn& state; //!< BC state config: sidesets + statefn
52 : : const std::vector< tk::StateFn >& fn; //!< BC state functions
53 : : const std::vector< tk::StateFn >& gfn; //!< BC gradient functions
54 : : std::size_t c; //!< Counts BC types configured
55 : : //! Constructor
56 : : ConfigBC( BCStateFn& s,
57 : : const std::vector< tk::StateFn >& f,
58 : 472 : const std::vector< tk::StateFn >& gf ) :
59 [ + - ][ + - ]: 472 : state(s), fn(f), gfn(gf), c(0) {}
[ + - ][ + - ]
[ - - ][ + - ]
[ + - ][ - - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ + - ][ + - ]
[ + - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ + - ][ + - ]
[ + - ][ - - ]
60 : : //! Function to call for each BC type
61 : 2832 : template< typename U > void operator()( brigand::type_<U> ) {
62 : : std::vector< std::size_t > cfg, v;
63 : : // collect sidesets across all meshes
64 [ + + ]: 5664 : for (const auto& ibc : g_inputdeck.get< tag::bc >()) {
65 [ + - ]: 2832 : v.insert(v.end(), ibc.get< U >().begin(), ibc.get< U >().end());
66 : : }
67 [ + + ][ + - ]: 2832 : if (v.size() > 0) cfg = v;
68 : : Assert( fn.size() > c, "StateFn missing for BC type" );
69 [ + - ][ - - ]: 2832 : state.push_back( { cfg, fn[c], gfn[c] } );
70 [ + + ]: 2832 : ++c;
71 : 2832 : }
72 : : };
73 : :
74 : : //! \brief Extract information on inlet BCs, which have a different structure
75 : : //! than other BCs
76 : : void ConfigInletBC( BCStateFn&,
77 : : const tk::StateFn&,
78 : : const tk::StateFn& );
79 : :
80 : : //! State function for invalid/un-configured boundary conditions
81 : : [[noreturn]] tk::StateFn::result_type
82 : : invalidBC( ncomp_t, const std::vector< EOS >&,
83 : : const std::vector< tk::real >&, tk::real, tk::real, tk::real,
84 : : tk::real, const std::array< tk::real, 3> & );
85 : :
86 : : //! \brief Partial differential equation base for discontinuous Galerkin PDEs
87 : : //! \details This class uses runtime polymorphism without client-side
88 : : //! inheritance: inheritance is confined to the internals of the this class,
89 : : //! invisible to client-code. The class exclusively deals with ownership
90 : : //! enabling client-side value semantics. Credit goes to Sean Parent at Adobe:
91 : : //! https://github.com/sean-parent/sean-parent.github.com/wiki/
92 : : //! Papers-and-Presentations. For example client code that models a DGPDE,
93 : : //! see inciter::CompFlow.
94 : 396 : class DGPDE {
95 : :
96 : : private:
97 : : using ncomp_t = tk::ncomp_t;
98 : :
99 : : public:
100 : : //! Default constructor taking no arguments for Charm++
101 : : explicit DGPDE() = default;
102 : :
103 : : //! \brief Constructor taking an object modeling Concept.
104 : : //! \details The object of class T comes pre-constructed.
105 : : //! \param[in] x Instantiated object of type T given by the template
106 : : //! argument.
107 : : template< typename T > explicit DGPDE( T x ) :
108 : : self( std::make_unique< Model<T> >( std::move(x) ) ) {}
109 : :
110 : : //! \brief Constructor taking a function pointer to a constructor of an
111 : : //! object modeling Concept.
112 : : //! \details Passing std::function allows late execution of the constructor,
113 : : //! i.e., as late as inside this class' constructor, and thus usage from
114 : : //! a factory. Note that there are at least two different ways of using
115 : : //! this constructor:
116 : : //! - Bind T's constructor arguments and place it in std::function<T()>
117 : : //! and passing no arguments as args.... This case then instantiates the
118 : : //! model via its constructor and stores it in here.
119 : : //! - Bind a single placeholder argument to T's constructor and pass it in
120 : : //! as host's args..., which then forwards it to model's constructor. This
121 : : //! allows late binding, i.e., binding the argument only here.
122 : : //! \see See also the wrapper tk::recordModel() which does the former and
123 : : //! tk::recordModelLate() which does the latter, both defined in
124 : : //! src/Base/Factory.h.
125 : : //! \param[in] x Function pointer to a constructor of an object modeling
126 : : //! Concept.
127 : : //! \param[in] args Zero or more constructor arguments
128 : : template< typename T, typename...Args >
129 [ - + ]: 396 : explicit DGPDE( std::function<T(Args...)> x, Args&&... args ) :
130 : : self( std::make_unique< Model<T> >(
131 [ + - ]: 544 : std::move( x( std::forward<Args>(args)... ) ) ) ) {}
132 : :
133 : : //! Public interface to find number of primitive quantities for the diff eq
134 : : std::size_t nprim() const
135 [ + - ][ - - ]: 889 : { return self->nprim(); }
136 : :
137 : : //! Public interface to find number of materials for the diff eq
138 : : std::size_t nmat() const
139 : : { return self->nmat(); }
140 : :
141 : : //! Public interface to find Dofs for each equation in pde system
142 : : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
143 [ + - ][ - - ]: 889 : { return self->numEquationDofs(numEqDof); }
144 : :
145 : : //! Public interface to find how 'stiff equations', which are the inverse
146 : : //! deformation equations because of plasticity
147 : : std::size_t nstiffeq() const
148 [ + - ][ + - ]: 3556 : { return self->nstiffeq(); }
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
149 : :
150 : : //! Public interface to find how 'nonstiff equations', which are the inverse
151 : : //! deformation equations because of plasticity
152 : : std::size_t nnonstiffeq() const
153 [ + - ][ + - ]: 1778 : { return self->nnonstiffeq(); }
[ - - ][ - - ]
154 : :
155 : : //! Public function to locate the stiff equations
156 : : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
157 : 0 : { return self->setStiffEqIdx( stiffEqIdx ); }
158 : :
159 : : //! Public function to locate the nonstiff equations
160 : : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
161 : 0 : { return self->setNonStiffEqIdx( nonStiffEqIdx ); }
162 : :
163 : : //! Public interface to determine elements that lie inside the IC box
164 : : void IcBoxElems( const tk::Fields& geoElem,
165 : : std::size_t nielem,
166 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
167 : 889 : { self->IcBoxElems( geoElem, nielem, inbox ); }
168 : :
169 : : //! Public interface to setting the initial conditions for the diff eq
170 : : void initialize(
171 : : const tk::Fields& L,
172 : : const std::vector< std::size_t >& inpoel,
173 : : const tk::UnsMesh::Coords& coord,
174 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
175 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
176 : : elemblkid,
177 : : tk::Fields& unk,
178 : : tk::real t,
179 : : const std::size_t nielem ) const
180 [ + - ][ + - ]: 1021 : { self->initialize( L, inpoel, coord, inbox, elemblkid, unk, t, nielem ); }
181 : :
182 : : //! Public interface for computing density constraint
183 : : void computeDensityConstr( std::size_t nelem,
184 : : tk::Fields& unk,
185 : : std::vector< tk::real >& densityConstr) const
186 [ + - ]: 2198 : { self->computeDensityConstr( nelem, unk, densityConstr); }
187 : :
188 : : //! Public interface to computing the left-hand side matrix for the diff eq
189 : : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const
190 [ + - ][ + - ]: 1021 : { self->lhs( geoElem, l ); }
191 : :
192 : : //! Public interface to updating the interface cells for the diff eq
193 : : void updateInterfaceCells( tk::Fields& unk,
194 : : std::size_t nielem,
195 : : std::vector< std::size_t >& ndofel,
196 : : std::vector< std::size_t >& interface ) const
197 : 71805 : { self->updateInterfaceCells( unk, nielem, ndofel, interface ); }
198 : :
199 : : //! Public interface to updating the primitives for the diff eq
200 : : void updatePrimitives( const tk::Fields& unk,
201 : : const tk::Fields& L,
202 : : const tk::Fields& geoElem,
203 : : tk::Fields& prim,
204 : : std::size_t nielem,
205 : : std::vector< std::size_t >& ndofel ) const
206 : 72694 : { self->updatePrimitives( unk, L, geoElem, prim, nielem, ndofel ); }
207 : :
208 : : //! Public interface to cleaning up trace materials for the diff eq
209 : : void cleanTraceMaterial( tk::real t,
210 : : const tk::Fields& geoElem,
211 : : tk::Fields& unk,
212 : : tk::Fields& prim,
213 : : std::size_t nielem ) const
214 : 71805 : { self->cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }
215 : :
216 : : //! Public interface to reconstructing the second-order solution
217 : : void reconstruct( tk::real t,
218 : : const tk::Fields& geoFace,
219 : : const tk::Fields& geoElem,
220 : : const inciter::FaceData& fd,
221 : : const std::map< std::size_t, std::vector< std::size_t > >&
222 : : esup,
223 : : const std::vector< std::size_t >& inpoel,
224 : : const tk::UnsMesh::Coords& coord,
225 : : tk::Fields& U,
226 : : tk::Fields& P,
227 : : const bool pref,
228 : : const std::vector< std::size_t >& ndofel ) const
229 : : {
230 : 43965 : self->reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P,
231 : 43965 : pref, ndofel );
232 : 43965 : }
233 : :
234 : : //! Public interface to limiting the second-order solution
235 : : void limit( tk::real t,
236 : : const bool pref,
237 : : const tk::Fields& geoFace,
238 : : const tk::Fields& geoElem,
239 : : const inciter::FaceData& fd,
240 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
241 : : const std::vector< std::size_t >& inpoel,
242 : : const tk::UnsMesh::Coords& coord,
243 : : const std::vector< std::size_t >& ndofel,
244 : : const std::vector< std::size_t >& gid,
245 : : const std::unordered_map< std::size_t, std::size_t >& bid,
246 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
247 : : const std::vector< std::vector<tk::real> >& pNodalExtrm,
248 : : const std::vector< std::vector<tk::real> >& mtInv,
249 : : tk::Fields& U,
250 : : tk::Fields& P,
251 : : std::vector< std::size_t >& shockmarker ) const
252 : : {
253 : 43965 : self->limit( t, pref, geoFace, geoElem, fd, esup, inpoel, coord, ndofel,
254 : 43965 : gid, bid, uNodalExtrm, pNodalExtrm, mtInv, U, P, shockmarker );
255 : : }
256 : :
257 : : //! Public interface to update the conservative variable solution
258 : : void CPL( const tk::Fields& prim,
259 : : const tk::Fields& geoElem,
260 : : const std::vector< std::size_t >& inpoel,
261 : : const tk::UnsMesh::Coords& coord,
262 : : tk::Fields& unk,
263 : : std::size_t nielem ) const
264 : : {
265 : 39915 : self->CPL( prim, geoElem, inpoel, coord, unk, nielem );
266 : 39915 : }
267 : :
268 : : //! Public interface to reset the high order solution for p-adaptive scheme
269 : : void resetAdapSol( const inciter::FaceData& fd,
270 : : tk::Fields& unk,
271 : : tk::Fields& prim,
272 : : const std::vector< std::size_t >& ndofel ) const
273 : : {
274 : 1770 : self->resetAdapSol( fd, unk, prim, ndofel );
275 : 1770 : }
276 : :
277 : : //! Public interface to getting the cell-averaged deformation gradients
278 : : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
279 : : const tk::Fields& U,
280 : : std::size_t nielem ) const
281 : : {
282 : : return self->cellAvgDeformGrad( U, nielem );
283 : : }
284 : :
285 : : //! Public interface to computing the P1 right-hand side vector
286 : : void rhs( tk::real t,
287 : : const bool pref,
288 : : const tk::Fields& geoFace,
289 : : const tk::Fields& geoElem,
290 : : const inciter::FaceData& fd,
291 : : const std::vector< std::size_t >& inpoel,
292 : : const std::vector< std::unordered_set< std::size_t > >& boxelems,
293 : : const tk::UnsMesh::Coords& coord,
294 : : const tk::Fields& U,
295 : : const tk::Fields& P,
296 : : const std::vector< std::size_t >& ndofel,
297 : : const tk::real dt,
298 : : tk::Fields& R ) const
299 : : {
300 : 71805 : self->rhs( t, pref, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
301 : 71805 : ndofel, dt, R );
302 : : }
303 : :
304 : : //! Evaluate the adaptive indicator and mark the ndof for each element
305 : : void eval_ndof( std::size_t nunk,
306 : : const tk::UnsMesh::Coords& coord,
307 : : const std::vector< std::size_t >& inpoel,
308 : : const inciter::FaceData& fd,
309 : : const tk::Fields& unk,
310 : : const tk::Fields& prim,
311 : : inciter::ctr::PrefIndicatorType indicator,
312 : : std::size_t ndof,
313 : : std::size_t ndofmax,
314 : : tk::real tolref,
315 : : std::vector< std::size_t >& ndofel ) const
316 : : {
317 : 1636 : self->eval_ndof( nunk, coord, inpoel, fd, unk, prim, indicator, ndof,
318 : 1636 : ndofmax, tolref, ndofel );
319 : 1636 : }
320 : :
321 : : //! Public interface for computing the minimum time step size
322 : : tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
323 : : const std::vector< std::size_t >& inpoel,
324 : : const inciter::FaceData& fd,
325 : : const tk::Fields& geoFace,
326 : : const tk::Fields& geoElem,
327 : : const std::vector< std::size_t >& ndofel,
328 : : const tk::Fields& U,
329 : : const tk::Fields& P,
330 : : const std::size_t nielem ) const
331 : 2330 : { return self->dt( coord, inpoel, fd, geoFace, geoElem, ndofel, U,
332 : 2330 : P, nielem ); }
333 : :
334 : : //! Public interface for computing stiff terms for an element
335 : : void stiff_rhs( std::size_t e,
336 : : const tk::Fields& geoElem,
337 : : const std::vector< std::size_t >& inpoel,
338 : : const tk::UnsMesh::Coords& coord,
339 : : const tk::Fields& U,
340 : : const tk::Fields& P,
341 : : const std::vector< std::size_t >& ndofel,
342 : : tk::Fields& R ) const
343 [ - - ][ - - ]: 0 : { return self->stiff_rhs( e, geoElem, inpoel, coord, U, P, ndofel, R); }
344 : :
345 : : //! Public interface to returning maps of output var functions
346 : : std::map< std::string, tk::GetVarFn > OutVarFn() const
347 [ + - ][ + - ]: 4396 : { return self->OutVarFn(); }
348 : :
349 : : //! Public interface to returning analytic field output labels
350 : : std::vector< std::string > analyticFieldNames() const
351 : 1505 : { return self->analyticFieldNames(); }
352 : :
353 : : //! Public interface to returning time history field output labels
354 [ - - ]: 0 : std::vector< std::string > histNames() const { return self->histNames(); }
355 : :
356 : : //! Public interface to returning variable names
357 : : std::vector< std::string > names() const { return self->names(); }
358 : :
359 : : //! Public interface to returning surface field output
360 : : std::vector< std::vector< tk::real > >
361 : : surfOutput( const std::map< int, std::vector< std::size_t > >& bnd,
362 : : tk::Fields& U ) const
363 : : { return self->surfOutput( bnd, U ); }
364 : :
365 : : //! Public interface to return point history output
366 : : std::vector< std::vector< tk::real > >
367 : : histOutput( const std::vector< HistData >& h,
368 : : const std::vector< std::size_t >& inpoel,
369 : : const tk::UnsMesh::Coords& coord,
370 : : const tk::Fields& U,
371 : : const tk::Fields& P ) const
372 [ - - ]: 0 : { return self->histOutput( h, inpoel, coord, U, P ); }
373 : :
374 : : //! Public interface to returning analytic solution
375 : : tk::InitializeFn::result_type
376 : : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
377 : 1188068 : { return self->analyticSolution( xi, yi, zi, t ); }
378 : :
379 : : //! Public interface to returning the analytic solution for conserved vars
380 : : tk::InitializeFn::result_type
381 : : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
382 [ + - ]: 6231923 : { return self->solution( xi, yi, zi, t ); }
383 : :
384 : : //! Public interface to returning the specific total energy
385 : : tk::real
386 : : sp_totalenergy( std::size_t e, const tk::Fields& unk ) const
387 [ + - ]: 3380569 : { return self->sp_totalenergy( e, unk ); }
388 : :
389 : : //! Copy assignment
390 : : DGPDE& operator=( const DGPDE& x )
391 : : { DGPDE tmp(x); *this = std::move(tmp); return *this; }
392 : : //! Copy constructor
393 : : DGPDE( const DGPDE& x ) : self( x.self->copy() ) {}
394 : : //! Move assignment
395 : : DGPDE& operator=( DGPDE&& ) noexcept = default;
396 : : //! Move constructor
397 [ - - ]: 0 : DGPDE( DGPDE&& ) noexcept = default;
398 : :
399 : : private:
400 : : //! \brief Concept is a pure virtual base class specifying the requirements
401 : : //! of polymorphic objects deriving from it
402 : : struct Concept {
403 : : Concept() = default;
404 : 0 : Concept( const Concept& ) = default;
405 : : virtual ~Concept() = default;
406 : : virtual Concept* copy() const = 0;
407 : : virtual std::size_t nprim() const = 0;
408 : : virtual std::size_t nmat() const = 0;
409 : : virtual void numEquationDofs(std::vector< std::size_t >&) const = 0;
410 : : virtual std::size_t nstiffeq() const = 0;
411 : : virtual std::size_t nnonstiffeq() const = 0;
412 : : virtual void setStiffEqIdx( std::vector< std::size_t >& ) const = 0;
413 : : virtual void setNonStiffEqIdx( std::vector< std::size_t >& ) const = 0;
414 : : virtual void IcBoxElems( const tk::Fields&,
415 : : std::size_t,
416 : : std::vector< std::unordered_set< std::size_t > >& ) const = 0;
417 : : virtual void initialize(
418 : : const tk::Fields&,
419 : : const std::vector< std::size_t >&,
420 : : const tk::UnsMesh::Coords&,
421 : : const std::vector< std::unordered_set< std::size_t > >&,
422 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&,
423 : : tk::Fields&,
424 : : tk::real,
425 : : const std::size_t nielem ) const = 0;
426 : : virtual void computeDensityConstr( std::size_t nelem,
427 : : tk::Fields& unk,
428 : : std::vector< tk::real >& densityConstr)
429 : : const = 0;
430 : : virtual void lhs( const tk::Fields&, tk::Fields& ) const = 0;
431 : : virtual void updateInterfaceCells( tk::Fields&,
432 : : std::size_t,
433 : : std::vector< std::size_t >&,
434 : : std::vector< std::size_t >& ) const = 0;
435 : : virtual void updatePrimitives( const tk::Fields&,
436 : : const tk::Fields&,
437 : : const tk::Fields&,
438 : : tk::Fields&,
439 : : std::size_t,
440 : : std::vector< std::size_t >& ) const = 0;
441 : : virtual void cleanTraceMaterial( tk::real,
442 : : const tk::Fields&,
443 : : tk::Fields&,
444 : : tk::Fields&,
445 : : std::size_t ) const = 0;
446 : : virtual void reconstruct( tk::real,
447 : : const tk::Fields&,
448 : : const tk::Fields&,
449 : : const inciter::FaceData&,
450 : : const std::map< std::size_t,
451 : : std::vector< std::size_t > >&,
452 : : const std::vector< std::size_t >&,
453 : : const tk::UnsMesh::Coords&,
454 : : tk::Fields&,
455 : : tk::Fields&,
456 : : const bool,
457 : : const std::vector< std::size_t >& ) const = 0;
458 : : virtual void limit( tk::real,
459 : : const bool,
460 : : const tk::Fields&,
461 : : const tk::Fields&,
462 : : const inciter::FaceData&,
463 : : const std::map< std::size_t,
464 : : std::vector< std::size_t > >&,
465 : : const std::vector< std::size_t >&,
466 : : const tk::UnsMesh::Coords&,
467 : : const std::vector< std::size_t >&,
468 : : const std::vector< std::size_t >&,
469 : : const std::unordered_map< std::size_t, std::size_t >&,
470 : : const std::vector< std::vector<tk::real> >&,
471 : : const std::vector< std::vector<tk::real> >&,
472 : : const std::vector< std::vector<tk::real> >&,
473 : : tk::Fields&,
474 : : tk::Fields&,
475 : : std::vector< std::size_t >& ) const = 0;
476 : : virtual void CPL( const tk::Fields&,
477 : : const tk::Fields&,
478 : : const std::vector< std::size_t >&,
479 : : const tk::UnsMesh::Coords&,
480 : : tk::Fields&,
481 : : std::size_t ) const = 0;
482 : : virtual std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
483 : : const tk::Fields&,
484 : : std::size_t ) const = 0;
485 : : virtual void rhs( tk::real,
486 : : const bool,
487 : : const tk::Fields&,
488 : : const tk::Fields&,
489 : : const inciter::FaceData&,
490 : : const std::vector< std::size_t >&,
491 : : const std::vector< std::unordered_set< std::size_t > >&,
492 : : const tk::UnsMesh::Coords&,
493 : : const tk::Fields&,
494 : : const tk::Fields&,
495 : : const std::vector< std::size_t >&,
496 : : const tk::real,
497 : : tk::Fields& ) const = 0;
498 : : virtual void resetAdapSol( const inciter::FaceData&,
499 : : tk::Fields&,
500 : : tk::Fields&,
501 : : const std::vector< std::size_t >& )
502 : : const = 0;
503 : : virtual void eval_ndof( std::size_t,
504 : : const tk::UnsMesh::Coords&,
505 : : const std::vector< std::size_t >&,
506 : : const inciter::FaceData&,
507 : : const tk::Fields&,
508 : : const tk::Fields&,
509 : : inciter::ctr::PrefIndicatorType,
510 : : std::size_t,
511 : : std::size_t,
512 : : tk::real,
513 : : std::vector< std::size_t >& ) const = 0;
514 : : virtual tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
515 : : const std::vector< std::size_t >&,
516 : : const inciter::FaceData&,
517 : : const tk::Fields&,
518 : : const tk::Fields&,
519 : : const std::vector< std::size_t >&,
520 : : const tk::Fields&,
521 : : const tk::Fields&,
522 : : const std::size_t ) const = 0;
523 : : virtual void stiff_rhs( std::size_t,
524 : : const tk::Fields&,
525 : : const std::vector< std::size_t >&,
526 : : const tk::UnsMesh::Coords&,
527 : : const tk::Fields&,
528 : : const tk::Fields&,
529 : : const std::vector< std::size_t >&,
530 : : tk::Fields& ) const = 0;
531 : : virtual std::map< std::string, tk::GetVarFn > OutVarFn() const = 0;
532 : : virtual std::vector< std::string > analyticFieldNames() const = 0;
533 : : virtual std::vector< std::string > histNames() const = 0;
534 : : virtual std::vector< std::string > names() const = 0;
535 : : virtual std::vector< std::vector< tk::real > > surfOutput(
536 : : const std::map< int, std::vector< std::size_t > >&,
537 : : tk::Fields& ) const = 0;
538 : : virtual std::vector< std::vector< tk::real > > histOutput(
539 : : const std::vector< HistData >&,
540 : : const std::vector< std::size_t >&,
541 : : const tk::UnsMesh::Coords&,
542 : : const tk::Fields&,
543 : : const tk::Fields& ) const = 0;
544 : : virtual tk::InitializeFn::result_type analyticSolution(
545 : : tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
546 : : virtual tk::InitializeFn::result_type solution(
547 : : tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
548 : : virtual tk::real sp_totalenergy(
549 : : std::size_t, const tk::Fields& ) const = 0;
550 : : };
551 : :
552 : : //! \brief Model models the Concept above by deriving from it and overriding
553 : : //! the virtual functions required by Concept
554 : : template< typename T >
555 [ - - ][ - - ]: 0 : struct Model : Concept {
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
556 : 396 : explicit Model( T x ) : data( std::move(x) ) {}
557 : 0 : Concept* copy() const override { return new Model( *this ); }
558 : 889 : std::size_t nprim() const override
559 : 889 : { return data.nprim(); }
560 : 0 : std::size_t nmat() const override
561 : 0 : { return data.nmat(); }
562 : 889 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const override
563 : 889 : { data.numEquationDofs(numEqDof); }
564 : 3556 : std::size_t nstiffeq() const override
565 : 3556 : { return data.nstiffeq(); }
566 : 1778 : std::size_t nnonstiffeq() const override
567 : 1778 : { return data.nnonstiffeq(); }
568 : 0 : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const override
569 : 0 : { data.setStiffEqIdx(stiffEqIdx); }
570 : 0 : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const override
571 : 0 : { data.setNonStiffEqIdx(nonStiffEqIdx); }
572 : 889 : void IcBoxElems( const tk::Fields& geoElem,
573 : : std::size_t nielem,
574 : : std::vector< std::unordered_set< std::size_t > >& inbox )
575 : 889 : const override { data.IcBoxElems( geoElem, nielem, inbox ); }
576 : 1021 : void initialize(
577 : : const tk::Fields& L,
578 : : const std::vector< std::size_t >& inpoel,
579 : : const tk::UnsMesh::Coords& coord,
580 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
581 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
582 : : elemblkid,
583 : : tk::Fields& unk,
584 : : tk::real t,
585 : : const std::size_t nielem )
586 : 1021 : const override { data.initialize( L, inpoel, coord, inbox, elemblkid, unk,
587 : 1021 : t, nielem ); }
588 : 2198 : void computeDensityConstr( std::size_t nelem,
589 : : tk::Fields& unk,
590 : : std::vector< tk::real >& densityConstr)
591 : : const override
592 : 2198 : { data.computeDensityConstr( nelem, unk, densityConstr ); }
593 : 1021 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const override
594 : 1021 : { data.lhs( geoElem, l ); }
595 : 71805 : void updateInterfaceCells( tk::Fields& unk,
596 : : std::size_t nielem,
597 : : std::vector< std::size_t >& ndofel,
598 : : std::vector< std::size_t >& interface )
599 : 71805 : const override { data.updateInterfaceCells( unk, nielem, ndofel, interface ); }
600 : 72694 : void updatePrimitives( const tk::Fields& unk,
601 : : const tk::Fields& L,
602 : : const tk::Fields& geoElem,
603 : : tk::Fields& prim,
604 : : std::size_t nielem,
605 : : std::vector< std::size_t >& ndofel )
606 : : const override {
607 : 5270 : data.updatePrimitives( unk, L, geoElem, prim, nielem, ndofel );
608 : 72694 : }
609 : 71805 : void cleanTraceMaterial( tk::real t,
610 : : const tk::Fields& geoElem,
611 : : tk::Fields& unk,
612 : : tk::Fields& prim,
613 : : std::size_t nielem )
614 : 71805 : const override { data.cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }
615 : 43965 : void reconstruct( tk::real t,
616 : : const tk::Fields& geoFace,
617 : : const tk::Fields& geoElem,
618 : : const inciter::FaceData& fd,
619 : : const std::map< std::size_t,
620 : : std::vector< std::size_t > >& esup,
621 : : const std::vector< std::size_t >& inpoel,
622 : : const tk::UnsMesh::Coords& coord,
623 : : tk::Fields& U,
624 : : tk::Fields& P,
625 : : const bool pref,
626 : : const std::vector< std::size_t >& ndofel )const override
627 : : {
628 : 43965 : data.reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P,
629 : : pref, ndofel );
630 : 43965 : }
631 : 43965 : void limit( tk::real t,
632 : : const bool pref,
633 : : const tk::Fields& geoFace,
634 : : const tk::Fields& geoElem,
635 : : const inciter::FaceData& fd,
636 : : const std::map< std::size_t, std::vector< std::size_t > >&
637 : : esup,
638 : : const std::vector< std::size_t >& inpoel,
639 : : const tk::UnsMesh::Coords& coord,
640 : : const std::vector< std::size_t >& ndofel,
641 : : const std::vector< std::size_t >& gid,
642 : : const std::unordered_map< std::size_t, std::size_t >& bid,
643 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
644 : : const std::vector< std::vector<tk::real> >& pNodalExtrm,
645 : : const std::vector< std::vector<tk::real> >& mtInv,
646 : : tk::Fields& U,
647 : : tk::Fields& P,
648 : : std::vector< std::size_t >& shockmarker ) const override
649 : : {
650 : 43965 : data.limit( t, pref, geoFace, geoElem, fd, esup, inpoel, coord, ndofel, gid,
651 : : bid, uNodalExtrm, pNodalExtrm, mtInv, U, P, shockmarker );
652 : 43965 : }
653 : 39915 : void CPL( const tk::Fields& prim,
654 : : const tk::Fields& geoElem,
655 : : const std::vector< std::size_t >& inpoel,
656 : : const tk::UnsMesh::Coords& coord,
657 : : tk::Fields& unk,
658 : : std::size_t nielem ) const override
659 : : {
660 : : data.CPL( prim, geoElem, inpoel, coord, unk, nielem );
661 : 39915 : }
662 : 1770 : void resetAdapSol( const inciter::FaceData& fd,
663 : : tk::Fields& unk,
664 : : tk::Fields& prim,
665 : : const std::vector< std::size_t >& ndofel ) const override
666 : : {
667 : 1770 : data.resetAdapSol( fd, unk, prim, ndofel );
668 : 1770 : }
669 : 0 : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
670 : : const tk::Fields& U,
671 : : std::size_t nielem ) const override
672 : : {
673 : 0 : return data.cellAvgDeformGrad( U, nielem );
674 : : }
675 : 71805 : void rhs(
676 : : tk::real t,
677 : : const bool pref,
678 : : const tk::Fields& geoFace,
679 : : const tk::Fields& geoElem,
680 : : const inciter::FaceData& fd,
681 : : const std::vector< std::size_t >& inpoel,
682 : : const std::vector< std::unordered_set< std::size_t > >& boxelems,
683 : : const tk::UnsMesh::Coords& coord,
684 : : const tk::Fields& U,
685 : : const tk::Fields& P,
686 : : const std::vector< std::size_t >& ndofel,
687 : : const tk::real dt,
688 : : tk::Fields& R ) const override
689 : : {
690 : 71805 : data.rhs( t, pref, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
691 : : ndofel, dt, R );
692 : 71805 : }
693 : 1636 : void eval_ndof( std::size_t nunk,
694 : : const tk::UnsMesh::Coords& coord,
695 : : const std::vector< std::size_t >& inpoel,
696 : : const inciter::FaceData& fd,
697 : : const tk::Fields& unk,
698 : : const tk::Fields& prim,
699 : : inciter::ctr::PrefIndicatorType indicator,
700 : : std::size_t ndof,
701 : : std::size_t ndofmax,
702 : : tk::real tolref,
703 : : std::vector< std::size_t >& ndofel ) const override
704 : 1636 : { data.eval_ndof( nunk, coord, inpoel, fd, unk, prim, indicator, ndof,
705 : 1636 : ndofmax, tolref, ndofel ); }
706 : 2330 : tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
707 : : const std::vector< std::size_t >& inpoel,
708 : : const inciter::FaceData& fd,
709 : : const tk::Fields& geoFace,
710 : : const tk::Fields& geoElem,
711 : : const std::vector< std::size_t >& ndofel,
712 : : const tk::Fields& U,
713 : : const tk::Fields& P,
714 : : const std::size_t nielem ) const override
715 : : { return data.dt( coord, inpoel, fd, geoFace, geoElem, ndofel,
716 : 2330 : U, P, nielem ); }
717 : 0 : void stiff_rhs( std::size_t e,
718 : : const tk::Fields& geoElem,
719 : : const std::vector< std::size_t >& inpoel,
720 : : const tk::UnsMesh::Coords& coord,
721 : : const tk::Fields& U,
722 : : const tk::Fields& P,
723 : : const std::vector< std::size_t >& ndofel,
724 : : tk::Fields& R ) const override
725 : 0 : { return data.stiff_rhs( e, geoElem, inpoel, coord, U, P, ndofel, R ); }
726 : 4396 : std::map< std::string, tk::GetVarFn > OutVarFn() const override
727 : 4396 : { return data.OutVarFn(); }
728 : 1505 : std::vector< std::string > analyticFieldNames() const override
729 : 1505 : { return data.analyticFieldNames(); }
730 : 0 : std::vector< std::string > histNames() const override
731 : 0 : { return data.histNames(); }
732 : 88 : std::vector< std::string > names() const override
733 : 88 : { return data.names(); }
734 : 0 : std::vector< std::vector< tk::real > > surfOutput(
735 : : const std::map< int, std::vector< std::size_t > >& bnd,
736 : : tk::Fields& U ) const override
737 : 0 : { return data.surfOutput( bnd, U ); }
738 : 0 : std::vector< std::vector< tk::real > > histOutput(
739 : : const std::vector< HistData >& h,
740 : : const std::vector< std::size_t >& inpoel,
741 : : const tk::UnsMesh::Coords& coord,
742 : : const tk::Fields& U,
743 : : const tk::Fields& P ) const override
744 : 0 : { return data.histOutput( h, inpoel, coord, U, P ); }
745 : : tk::InitializeFn::result_type
746 : 1188068 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
747 : 1188068 : const override { return data.analyticSolution( xi, yi, zi, t ); }
748 : : tk::InitializeFn::result_type
749 : 6231923 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
750 : 6231923 : const override { return data.solution( xi, yi, zi, t ); }
751 : 3380569 : tk::real sp_totalenergy( std::size_t e, const tk::Fields& unk )
752 : 3380569 : const override { return data.sp_totalenergy( e, unk ); }
753 : : T data;
754 : : };
755 : :
756 : : std::unique_ptr< Concept > self; //!< Base pointer used polymorphically
757 : : };
758 : :
759 : : } // inciter::
760 : :
761 : : #endif // DGPDE_h
|
