Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/FVPDE.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 finite volume PDEs
9 : : \details This file defines a generic partial differential equation (PDE)
10 : : class for PDEs that use finite volume 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 FVPDE_h
20 : : #define FVPDE_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 : :
43 : : //! \brief Partial differential equation base for discontinuous Galerkin PDEs
44 : : //! \details This class uses runtime polymorphism without client-side
45 : : //! inheritance: inheritance is confined to the internals of the this class,
46 : : //! invisible to client-code. The class exclusively deals with ownership
47 : : //! enabling client-side value semantics. Credit goes to Sean Parent at Adobe:
48 : : //! https://github.com/sean-parent/sean-parent.github.com/wiki/
49 : : //! Papers-and-Presentations. For example client code that models a FVPDE,
50 : : //! see inciter::CompFlow.
51 : : class FVPDE {
52 : :
53 : : private:
54 : : using ncomp_t = tk::ncomp_t;
55 : :
56 : : public:
57 : : //! Default constructor taking no arguments for Charm++
58 : 0 : explicit FVPDE() = default;
59 : :
60 : : //! \brief Constructor taking an object modeling Concept.
61 : : //! \details The object of class T comes pre-constructed.
62 : : //! \param[in] x Instantiated object of type T given by the template
63 : : //! argument.
64 : : template< typename T > explicit FVPDE( T x ) :
65 : : self( std::make_unique< Model<T> >( std::move(x) ) ) {}
66 : :
67 : : //! \brief Constructor taking a function pointer to a constructor of an
68 : : //! object modeling Concept.
69 : : //! \details Passing std::function allows late execution of the constructor,
70 : : //! i.e., as late as inside this class' constructor, and thus usage from
71 : : //! a factory. Note that there are at least two different ways of using
72 : : //! this constructor:
73 : : //! - Bind T's constructor arguments and place it in std::function<T()>
74 : : //! and passing no arguments as args.... This case then instantiates the
75 : : //! model via its constructor and stores it in here.
76 : : //! - Bind a single placeholder argument to T's constructor and pass it in
77 : : //! as host's args..., which then forwards it to model's constructor. This
78 : : //! allows late binding, i.e., binding the argument only here.
79 : : //! \see See also the wrapper tk::recordModel() which does the former and
80 : : //! tk::recordModelLate() which does the latter, both defined in
81 : : //! src/Base/Factory.h.
82 : : //! \param[in] x Function pointer to a constructor of an object modeling
83 : : //! Concept.
84 : : //! \param[in] args Zero or more constructor arguments
85 : : template< typename T, typename...Args >
86 : 81 : explicit FVPDE( std::function<T(Args...)> x, Args&&... args ) :
87 : : self( std::make_unique< Model<T> >(
88 [ + - ]: 81 : std::move( x( std::forward<Args>(args)... ) ) ) ) {}
89 : :
90 : : //! Public interface to find number of primitive quantities for the diff eq
91 : 298 : std::size_t nprim() const
92 : 298 : { return self->nprim(); }
93 : :
94 : : //! Public interface to find number of materials for the diff eq
95 : : std::size_t nmat() const
96 : : { return self->nmat(); }
97 : :
98 : : //! Public interface to determine elements that lie inside the IC box
99 : 298 : void IcBoxElems( const tk::Fields& geoElem,
100 : : std::size_t nielem,
101 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
102 : 298 : { self->IcBoxElems( geoElem, nielem, inbox ); }
103 : :
104 : : //! Public interface to setting the initial conditions for the diff eq
105 : 298 : void initialize(
106 : : const tk::Fields& L,
107 : : const std::vector< std::size_t >& inpoel,
108 : : const tk::UnsMesh::Coords& coord,
109 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
110 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
111 : : elemblkid,
112 : : tk::Fields& unk,
113 : : tk::real t,
114 : : const std::size_t nielem ) const
115 : 298 : { self->initialize( L, inpoel, coord, inbox, elemblkid, unk, t, nielem ); }
116 : :
117 : : //! Public interface to computing the left-hand side matrix for the diff eq
118 : 298 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const
119 : 298 : { self->lhs( geoElem, l ); }
120 : :
121 : : //! Public interface to updating the primitives for the diff eq
122 : 1966 : void updatePrimitives( const tk::Fields& unk,
123 : : tk::Fields& prim,
124 : : std::size_t nielem ) const
125 : 1966 : { self->updatePrimitives( unk, prim, nielem ); }
126 : :
127 : : //! Public interface to cleaning up trace materials for the diff eq
128 : 1668 : void cleanTraceMaterial( tk::real t,
129 : : const tk::Fields& geoElem,
130 : : tk::Fields& unk,
131 : : tk::Fields& prim,
132 : : std::size_t nielem ) const
133 : 1668 : { self->cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }
134 : :
135 : : //! Public interface to reconstructing the second-order solution
136 : 1668 : void reconstruct( const tk::Fields& geoElem,
137 : : const inciter::FaceData& fd,
138 : : const std::map< std::size_t, std::vector< std::size_t > >&
139 : : esup,
140 : : const std::vector< std::size_t >& inpoel,
141 : : const tk::UnsMesh::Coords& coord,
142 : : tk::Fields& U,
143 : : tk::Fields& P ) const
144 : : {
145 : 1668 : self->reconstruct( geoElem, fd, esup, inpoel, coord, U, P );
146 : 1668 : }
147 : :
148 : : //! Public interface to limiting the second-order solution
149 : 1668 : void limit( const tk::Fields& geoFace,
150 : : const inciter::FaceData& fd,
151 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
152 : : const std::vector< std::size_t >& inpoel,
153 : : const tk::UnsMesh::Coords& coord,
154 : : const std::vector< int >& srcFlag,
155 : : tk::Fields& U,
156 : : tk::Fields& P ) const
157 : : {
158 : 1668 : self->limit( geoFace, fd, esup, inpoel, coord, srcFlag, U, P );
159 : 1668 : }
160 : :
161 : : //! Public interface to computing the P1 right-hand side vector
162 : 1668 : void rhs( tk::real t,
163 : : const tk::Fields& geoFace,
164 : : const tk::Fields& geoElem,
165 : : const inciter::FaceData& fd,
166 : : const std::vector< std::size_t >& inpoel,
167 : : const tk::UnsMesh::Coords& coord,
168 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
169 : : elemblkid,
170 : : const tk::Fields& U,
171 : : const tk::Fields& P,
172 : : tk::Fields& R,
173 : : std::vector< int >& srcFlag ) const
174 : : {
175 : 1668 : self->rhs( t, geoFace, geoElem, fd, inpoel, coord, elemblkid, U, P, R,
176 : 1668 : srcFlag );
177 : 1668 : }
178 : :
179 : : //! Public interface for computing the minimum time step size
180 : 250 : tk::real dt( const inciter::FaceData& fd,
181 : : const tk::Fields& geoFace,
182 : : const tk::Fields& geoElem,
183 : : const tk::Fields& U,
184 : : const tk::Fields& P,
185 : : const std::size_t nielem,
186 : : const std::vector< int >& srcFlag,
187 : : std::vector< tk::real >& local_dte ) const
188 : 250 : { return self->dt( fd, geoFace, geoElem, U, P, nielem, srcFlag, local_dte );
189 : : }
190 : :
191 : : //! Public interface to returning maps of output var functions
192 : 56 : std::map< std::string, tk::GetVarFn > OutVarFn() const
193 : 56 : { return self->OutVarFn(); }
194 : :
195 : : //! Public interface to returning analytic field output labels
196 : 0 : std::vector< std::string > analyticFieldNames() const
197 : 0 : { return self->analyticFieldNames(); }
198 : :
199 : : //! Public interface to returning surface output labels
200 : 28 : std::vector< std::string > surfNames() const { return self->surfNames(); }
201 : :
202 : : //! Public interface to returning time history field output labels
203 : 0 : std::vector< std::string > histNames() const { return self->histNames(); }
204 : :
205 : : //! Public interface to returning variable names
206 : 21 : std::vector< std::string > names() const { return self->names(); }
207 : :
208 : : //! Public interface to returning surface field output
209 : : std::vector< std::vector< tk::real > >
210 : 28 : surfOutput( const inciter::FaceData& fd,
211 : : const tk::Fields& U,
212 : : const tk::Fields& P ) const
213 : 28 : { return self->surfOutput( fd, U, P ); }
214 : :
215 : : //! Public interface to return point history output
216 : : std::vector< std::vector< tk::real > >
217 : 0 : histOutput( const std::vector< HistData >& h,
218 : : const std::vector< std::size_t >& inpoel,
219 : : const tk::UnsMesh::Coords& coord,
220 : : const tk::Fields& U,
221 : : const tk::Fields& P ) const
222 : 0 : { return self->histOutput( h, inpoel, coord, U, P ); }
223 : :
224 : : //! Public interface to returning analytic solution
225 : : tk::InitializeFn::result_type
226 : 0 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
227 : 0 : { return self->analyticSolution( xi, yi, zi, t ); }
228 : :
229 : : //! Public interface to returning the analytic solution for conserved vars
230 : : tk::InitializeFn::result_type
231 : : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
232 : : { return self->solution( xi, yi, zi, t ); }
233 : :
234 : : //! Public interface to returning the specific total energy
235 : : tk::real
236 : : sp_totalenergy( std::size_t e, const tk::Fields& unk ) const
237 : : { return self->sp_totalenergy( e, unk ); }
238 : :
239 : : //! Public interface to returning the relevant sound speed in each cell
240 : : void
241 : 28 : soundspeed(
242 : : std::size_t nielem,
243 : : const tk::Fields& U,
244 : : const tk::Fields& P,
245 : : std::vector< tk::real >& ss ) const
246 : 28 : { return self->soundspeed( nielem, U, P, ss ); }
247 : :
248 : : //! Copy assignment
249 : : FVPDE& operator=( const FVPDE& x )
250 : : { FVPDE tmp(x); *this = std::move(tmp); return *this; }
251 : : //! Copy constructor
252 : : FVPDE( const FVPDE& x ) : self( x.self->copy() ) {}
253 : : //! Move assignment
254 : : FVPDE& operator=( FVPDE&& ) noexcept = default;
255 : : //! Move constructor
256 : 81 : FVPDE( FVPDE&& ) noexcept = default;
257 : :
258 : : private:
259 : : //! \brief Concept is a pure virtual base class specifying the requirements
260 : : //! of polymorphic objects deriving from it
261 : : struct Concept {
262 : 81 : Concept() = default;
263 : 0 : Concept( const Concept& ) = default;
264 : 81 : virtual ~Concept() = default;
265 : : virtual Concept* copy() const = 0;
266 : : virtual std::size_t nprim() const = 0;
267 : : virtual std::size_t nmat() const = 0;
268 : : virtual void IcBoxElems( const tk::Fields&,
269 : : std::size_t,
270 : : std::vector< std::unordered_set< std::size_t > >& ) const = 0;
271 : : virtual void initialize(
272 : : const tk::Fields&,
273 : : const std::vector< std::size_t >&,
274 : : const tk::UnsMesh::Coords&,
275 : : const std::vector< std::unordered_set< std::size_t > >&,
276 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&,
277 : : tk::Fields&,
278 : : tk::real,
279 : : const std::size_t nielem ) const = 0;
280 : : virtual void lhs( const tk::Fields&, tk::Fields& ) const = 0;
281 : : virtual void updatePrimitives( const tk::Fields&,
282 : : tk::Fields&,
283 : : std::size_t ) const = 0;
284 : : virtual void cleanTraceMaterial( tk::real,
285 : : const tk::Fields&,
286 : : tk::Fields&,
287 : : tk::Fields&,
288 : : std::size_t ) const = 0;
289 : : virtual void reconstruct( const tk::Fields&,
290 : : const inciter::FaceData&,
291 : : const std::map< std::size_t,
292 : : std::vector< std::size_t > >&,
293 : : const std::vector< std::size_t >&,
294 : : const tk::UnsMesh::Coords&,
295 : : tk::Fields&,
296 : : tk::Fields& ) const = 0;
297 : : virtual void limit( const tk::Fields&,
298 : : const inciter::FaceData&,
299 : : const std::map< std::size_t,
300 : : std::vector< std::size_t > >&,
301 : : const std::vector< std::size_t >&,
302 : : const tk::UnsMesh::Coords&,
303 : : const std::vector< int >&,
304 : : tk::Fields&,
305 : : tk::Fields& ) const = 0;
306 : : virtual void rhs( tk::real,
307 : : const tk::Fields&,
308 : : const tk::Fields&,
309 : : const inciter::FaceData&,
310 : : const std::vector< std::size_t >&,
311 : : const tk::UnsMesh::Coords&,
312 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&,
313 : : const tk::Fields&,
314 : : const tk::Fields&,
315 : : tk::Fields&,
316 : : std::vector< int >& ) const = 0;
317 : : virtual tk::real dt( const inciter::FaceData&,
318 : : const tk::Fields&,
319 : : const tk::Fields&,
320 : : const tk::Fields&,
321 : : const tk::Fields&,
322 : : const std::size_t,
323 : : const std::vector< int >&,
324 : : std::vector< tk::real >& ) const = 0;
325 : : virtual std::map< std::string, tk::GetVarFn > OutVarFn() const = 0;
326 : : virtual std::vector< std::string > analyticFieldNames() const = 0;
327 : : virtual std::vector< std::string > surfNames() const = 0;
328 : : virtual std::vector< std::string > histNames() const = 0;
329 : : virtual std::vector< std::string > names() const = 0;
330 : : virtual std::vector< std::vector< tk::real > > surfOutput(
331 : : const inciter::FaceData&,
332 : : const tk::Fields&,
333 : : const tk::Fields& ) const = 0;
334 : : virtual std::vector< std::vector< tk::real > > histOutput(
335 : : const std::vector< HistData >&,
336 : : const std::vector< std::size_t >&,
337 : : const tk::UnsMesh::Coords&,
338 : : const tk::Fields&,
339 : : const tk::Fields& ) const = 0;
340 : : virtual tk::InitializeFn::result_type analyticSolution(
341 : : tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
342 : : virtual tk::InitializeFn::result_type solution(
343 : : tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
344 : : virtual tk::real sp_totalenergy(
345 : : std::size_t, const tk::Fields& ) const = 0;
346 : : virtual void soundspeed(
347 : : std::size_t,
348 : : const tk::Fields&,
349 : : const tk::Fields&,
350 : : std::vector< tk::real >& ) const = 0;
351 : : };
352 : :
353 : : //! \brief Model models the Concept above by deriving from it and overriding
354 : : //! the virtual functions required by Concept
355 : : template< typename T >
356 : : struct Model : Concept {
357 : 81 : explicit Model( T x ) : data( std::move(x) ) {}
358 [ - - ]: 0 : Concept* copy() const override { return new Model( *this ); }
359 : 298 : std::size_t nprim() const override
360 : 298 : { return data.nprim(); }
361 : 0 : std::size_t nmat() const override
362 : 0 : { return data.nmat(); }
363 : 298 : void IcBoxElems( const tk::Fields& geoElem,
364 : : std::size_t nielem,
365 : : std::vector< std::unordered_set< std::size_t > >& inbox )
366 : 298 : const override { data.IcBoxElems( geoElem, nielem, inbox ); }
367 : 298 : void initialize(
368 : : const tk::Fields& L,
369 : : const std::vector< std::size_t >& inpoel,
370 : : const tk::UnsMesh::Coords& coord,
371 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
372 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
373 : : elemblkid,
374 : : tk::Fields& unk,
375 : : tk::real t,
376 : : const std::size_t nielem )
377 : 298 : const override { data.initialize( L, inpoel, coord, inbox, elemblkid, unk,
378 : 298 : t, nielem ); }
379 : 298 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const override
380 : 298 : { data.lhs( geoElem, l ); }
381 : 1966 : void updatePrimitives( const tk::Fields& unk,
382 : : tk::Fields& prim,
383 : : std::size_t nielem )
384 : 1966 : const override { data.updatePrimitives( unk, prim, nielem ); }
385 : 1668 : void cleanTraceMaterial( tk::real t,
386 : : const tk::Fields& geoElem,
387 : : tk::Fields& unk,
388 : : tk::Fields& prim,
389 : : std::size_t nielem )
390 : 1668 : const override { data.cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }
391 : 1668 : void reconstruct( const tk::Fields& geoElem,
392 : : const inciter::FaceData& fd,
393 : : const std::map< std::size_t,
394 : : std::vector< std::size_t > >& esup,
395 : : const std::vector< std::size_t >& inpoel,
396 : : const tk::UnsMesh::Coords& coord,
397 : : tk::Fields& U,
398 : : tk::Fields& P ) const override
399 : : {
400 : 1668 : data.reconstruct( geoElem, fd, esup, inpoel, coord, U, P );
401 : 1668 : }
402 : 1668 : void limit( const tk::Fields& geoFace,
403 : : const inciter::FaceData& fd,
404 : : const std::map< std::size_t, std::vector< std::size_t > >&
405 : : esup,
406 : : const std::vector< std::size_t >& inpoel,
407 : : const tk::UnsMesh::Coords& coord,
408 : : const std::vector< int >& srcFlag,
409 : : tk::Fields& U,
410 : : tk::Fields& P ) const override
411 : : {
412 : 1668 : data.limit( geoFace, fd, esup, inpoel, coord, srcFlag, U, P );
413 : 1668 : }
414 : 1668 : void rhs(
415 : : tk::real t,
416 : : const tk::Fields& geoFace,
417 : : const tk::Fields& geoElem,
418 : : const inciter::FaceData& fd,
419 : : const std::vector< std::size_t >& inpoel,
420 : : const tk::UnsMesh::Coords& coord,
421 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
422 : : elemblkid,
423 : : const tk::Fields& U,
424 : : const tk::Fields& P,
425 : : tk::Fields& R,
426 : : std::vector< int >& srcFlag ) const override
427 : : {
428 : 1668 : data.rhs( t, geoFace, geoElem, fd, inpoel, coord, elemblkid, U, P, R,
429 : : srcFlag );
430 : 1668 : }
431 : 250 : tk::real dt( const inciter::FaceData& fd,
432 : : const tk::Fields& geoFace,
433 : : const tk::Fields& geoElem,
434 : : const tk::Fields& U,
435 : : const tk::Fields& P,
436 : : const std::size_t nielem,
437 : : const std::vector< int >& srcFlag,
438 : : std::vector< tk::real >& local_dte ) const override
439 : : { return data.dt( fd, geoFace, geoElem, U, P, nielem, srcFlag,
440 : 250 : local_dte ); }
441 : 56 : std::map< std::string, tk::GetVarFn > OutVarFn() const override
442 : 56 : { return data.OutVarFn(); }
443 : 0 : std::vector< std::string > analyticFieldNames() const override
444 : 0 : { return data.analyticFieldNames(); }
445 : 28 : std::vector< std::string > surfNames() const override
446 : 28 : { return data.surfNames(); }
447 : 0 : std::vector< std::string > histNames() const override
448 : 0 : { return data.histNames(); }
449 : 21 : std::vector< std::string > names() const override
450 : 21 : { return data.names(); }
451 : 28 : std::vector< std::vector< tk::real > > surfOutput(
452 : : const inciter::FaceData& fd,
453 : : const tk::Fields& U,
454 : : const tk::Fields& P ) const override
455 : 28 : { return data.surfOutput( fd, U, P ); }
456 : 0 : std::vector< std::vector< tk::real > > histOutput(
457 : : const std::vector< HistData >& h,
458 : : const std::vector< std::size_t >& inpoel,
459 : : const tk::UnsMesh::Coords& coord,
460 : : const tk::Fields& U,
461 : : const tk::Fields& P ) const override
462 : 0 : { return data.histOutput( h, inpoel, coord, U, P ); }
463 : : tk::InitializeFn::result_type
464 : 0 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
465 : 0 : const override { return data.analyticSolution( xi, yi, zi, t ); }
466 : : tk::InitializeFn::result_type
467 : 0 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
468 : 0 : const override { return data.solution( xi, yi, zi, t ); }
469 : 0 : tk::real sp_totalenergy( std::size_t e, const tk::Fields& unk )
470 : 0 : const override { return data.sp_totalenergy( e, unk ); }
471 : 28 : void soundspeed(
472 : : std::size_t nielem,
473 : : const tk::Fields& U,
474 : : const tk::Fields& P,
475 : : std::vector< tk::real >& ss )
476 : 28 : const override { return data.soundspeed( nielem, U, P, ss ); }
477 : : T data;
478 : : };
479 : :
480 : : std::unique_ptr< Concept > self; //!< Base pointer used polymorphically
481 : : };
482 : :
483 : : } // inciter::
484 : :
485 : : #endif // FVPDE_h
|
