1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559 | // *****************************************************************************
/*!
\file src/PDE/DGPDE.hpp
\copyright 2012-2015 J. Bakosi,
2016-2018 Los Alamos National Security, LLC.,
2019-2021 Triad National Security, LLC.
All rights reserved. See the LICENSE file for details.
\brief Partial differential equation base for discontinuous Galerkin PDEs
\details This file defines a generic partial differential equation (PDE)
class for PDEs that use discontinuous Galerkin spatial discretization.
The class uses runtime polymorphism without client-side inheritance:
inheritance is confined to the internals of the class, invisible to
client-code. The class exclusively deals with ownership enabling client-side
value semantics. Credit goes to Sean Parent at Adobe:
https://github.com/sean-parent/sean-parent.github.com/wiki/
Papers-and-Presentations.
*/
// *****************************************************************************
#ifndef DGPDE_h
#define DGPDE_h
#include <array>
#include <string>
#include <vector>
#include <memory>
#include <unordered_set>
#include <functional>
#include "Types.hpp"
#include "Fields.hpp"
#include "FaceData.hpp"
#include "UnsMesh.hpp"
#include "Inciter/InputDeck/InputDeck.hpp"
#include "FunctionPrototypes.hpp"
#include "History.hpp"
namespace inciter {
extern ctr::InputDeck g_inputdeck;
using ncomp_t = kw::ncomp::info::expect::type;
using bcconf_t = kw::sideset::info::expect::type;
using BCStateFn =
std::vector< std::pair< std::vector< bcconf_t >, tk::StateFn > >;
//! Extract BC configuration ignoring if BC not specified
//! \note A more preferable way of catching errors such as this function
//! hides is during parsing, so that we don't even get here if BCs are
//! not correctly specified. For now we simply ignore if BCs are not
//! specified by allowing empty BC vectors from the user input.
template< class Eq > struct ConfigBC {
std::size_t system; //! Compflow system id
BCStateFn& state; //!< BC state config: sidesets + statefn
const std::vector< tk::StateFn >& fn; //!< BC state functions
std::size_t c; //!< Counts BC types configured
//! Constructor
ConfigBC( std::size_t sys,
BCStateFn& s,
const std::vector< tk::StateFn >& f ) :
system(sys), state(s), fn(f), c(0) {}
//! Function to call for each BC type
template< typename U > void operator()( brigand::type_<U> ) {
std::vector< bcconf_t > cfg;
const auto& v = g_inputdeck.get< tag::param, Eq, tag::bc, U >();
if (v.size() > system) cfg = v[system];<--- Variable 'cfg' is assigned a value that is never used.
Assert( fn.size() > c, "StateFn missing for BC type" );
state.push_back( { cfg, fn[c++] } );
}
};
//! State function for invalid/un-configured boundary conditions
[[noreturn]] tk::StateFn::result_type
invalidBC( ncomp_t, ncomp_t, const std::vector< tk::real >&,
tk::real, tk::real, tk::real, tk::real,
const std::array< tk::real, 3> & );
//! \brief Partial differential equation base for discontinuous Galerkin PDEs
//! \details This class uses runtime polymorphism without client-side
//! inheritance: inheritance is confined to the internals of the this class,
//! invisible to client-code. The class exclusively deals with ownership
//! enabling client-side value semantics. Credit goes to Sean Parent at Adobe:
//! https://github.com/sean-parent/sean-parent.github.com/wiki/
//! Papers-and-Presentations. For example client code that models a DGPDE,
//! see inciter::CompFlow.
class DGPDE {
private:
using ncomp_t = kw::ncomp::info::expect::type;
public:
//! Default constructor taking no arguments for Charm++
explicit DGPDE() = default;
//! \brief Constructor taking an object modeling Concept.
//! \details The object of class T comes pre-constructed.
//! \param[in] x Instantiated object of type T given by the template
//! argument.
template< typename T > explicit DGPDE( T x ) :
self( std::make_unique< Model<T> >( std::move(x) ) ) {}
//! \brief Constructor taking a function pointer to a constructor of an
//! object modeling Concept.
//! \details Passing std::function allows late execution of the constructor,
//! i.e., as late as inside this class' constructor, and thus usage from
//! a factory. Note that there are at least two different ways of using
//! this constructor:
//! - Bind T's constructor arguments and place it in std::function<T()>
//! and passing no arguments as args.... This case then instantiates the
//! model via its constructor and stores it in here.
//! - Bind a single placeholder argument to T's constructor and pass it in
//! as host's args..., which then forwards it to model's constructor. This
//! allows late binding, i.e., binding the argument only here.
//! \see See also the wrapper tk::recordModel() which does the former and
//! tk::recordModelLate() which does the latter, both defined in
//! src/Base/Factory.h.
//! \param[in] x Function pointer to a constructor of an object modeling
//! Concept.
//! \param[in] args Zero or more constructor arguments
template< typename T, typename...Args >
explicit DGPDE( std::function<T(Args...)> x, Args&&... args ) :
self( std::make_unique< Model<T> >(
std::move( x( std::forward<Args>(args)... ) ) ) ) {}
//! Public interface to find number of primitive quantities for the diff eq
std::size_t nprim() const
{ return self->nprim(); }
//! Public interface to find number of materials for the diff eq
std::size_t nmat() const
{ return self->nmat(); }
//! Public interface to find Dofs for each equation in pde system
void numEquationDofs(std::vector< std::size_t >& numEqDof) const
{ return self->numEquationDofs(numEqDof); }
//! Public interface to determine elements that lie inside the IC box
void IcBoxElems( const tk::Fields& geoElem,
std::size_t nielem,
std::vector< std::unordered_set< std::size_t > >& inbox ) const
{ self->IcBoxElems( geoElem, nielem, inbox ); }
//! Public interface to setting the initial conditions for the diff eq
void initialize(
const tk::Fields& L,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
const std::vector< std::unordered_set< std::size_t > >& inbox,
tk::Fields& unk,
tk::real t,
const std::size_t nielem ) const
{ self->initialize( L, inpoel, coord, inbox, unk, t, nielem ); }
//! Public interface to computing the left-hand side matrix for the diff eq
void lhs( const tk::Fields& geoElem, tk::Fields& l ) const
{ self->lhs( geoElem, l ); }
//! Public interface to updating the interface cells for the diff eq
void updateInterfaceCells( tk::Fields& unk,
std::size_t nielem,
std::vector< std::size_t >& ndofel ) const
{ self->updateInterfaceCells( unk, nielem, ndofel ); }
//! Public interface to updating the primitives for the diff eq
void updatePrimitives( const tk::Fields& unk,
const tk::Fields& L,
const tk::Fields& geoElem,
tk::Fields& prim,
std::size_t nielem ) const
{ self->updatePrimitives( unk, L, geoElem, prim, nielem ); }
//! Public interface to cleaning up trace materials for the diff eq
void cleanTraceMaterial( const tk::Fields& geoElem,
tk::Fields& unk,
tk::Fields& prim,
std::size_t nielem ) const
{ self->cleanTraceMaterial( geoElem, unk, prim, nielem ); }
//! Public interface to reconstructing the second-order solution
void reconstruct( tk::real t,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const inciter::FaceData& fd,
const std::map< std::size_t, std::vector< std::size_t > >&
esup,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
tk::Fields& U,
tk::Fields& P ) const
{
self->reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P );
}
//! Public interface to limiting the second-order solution
void limit( tk::real t,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const inciter::FaceData& fd,
const std::map< std::size_t, std::vector< std::size_t > >& esup,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
const std::vector< std::size_t >& ndofel,
const std::vector< std::size_t >& gid,
const std::unordered_map< std::size_t, std::size_t >& bid,
const std::vector< std::vector<tk::real> >& uNodalExtrm,
const std::vector< std::vector<tk::real> >& pNodalExtrm,
tk::Fields& U,
tk::Fields& P,
std::vector< std::size_t >& shockmarker ) const
{
self->limit( t, geoFace, geoElem, fd, esup, inpoel, coord, ndofel, gid,
bid, uNodalExtrm, pNodalExtrm, U, P, shockmarker );
}
//! Public interface to computing the P1 right-hand side vector
void rhs( tk::real t,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const inciter::FaceData& fd,
const std::vector< std::size_t >& inpoel,
const std::vector< std::unordered_set< std::size_t > >& boxelems,
const tk::UnsMesh::Coords& coord,
const tk::Fields& U,
const tk::Fields& P,
const std::vector< std::size_t >& ndofel,
tk::Fields& R ) const
{
self->rhs( t, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
ndofel, R );
}
//! Evaluate the adaptive indicator and mark the ndof for each element
void eval_ndof( std::size_t nunk,
const tk::UnsMesh::Coords& coord,
const std::vector< std::size_t >& inpoel,
const inciter::FaceData& fd,
const tk::Fields& unk,
inciter::ctr::PrefIndicatorType indicator,
std::size_t ndof,
std::size_t ndofmax,
tk::real tolref,
std::vector< std::size_t >& ndofel ) const
{
self->eval_ndof( nunk, coord, inpoel, fd, unk, indicator, ndof, ndofmax,
tolref, ndofel );
}
//! Public interface for computing the minimum time step size
tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
const std::vector< std::size_t >& inpoel,
const inciter::FaceData& fd,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const std::vector< std::size_t >& ndofel,
const tk::Fields& U,
const tk::Fields& P,
const std::size_t nielem ) const
{ return self->dt( coord, inpoel, fd, geoFace, geoElem, ndofel, U, P, nielem ); }
//! Public interface to returning analytic field output labels
std::vector< std::string > analyticFieldNames() const
{ return self->analyticFieldNames(); }
//! Public interface to returning time history field output labels
std::vector< std::string > histNames() const { return self->histNames(); }
//! Public interface to returning variable names
std::vector< std::string > names() const { return self->names(); }
//! Public interface to returning surface field output
std::vector< std::vector< tk::real > >
surfOutput( const std::map< int, std::vector< std::size_t > >& bnd,
tk::Fields& U ) const
{ return self->surfOutput( bnd, U ); }
//! Public interface to return point history output
std::vector< std::vector< tk::real > >
histOutput( const std::vector< HistData >& h,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
const tk::Fields& U,
const tk::Fields& P ) const
{ return self->histOutput( h, inpoel, coord, U, P ); }
//! Public interface to returning analytic solution
tk::InitializeFn::result_type
analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
{ return self->analyticSolution( xi, yi, zi, t ); }
//! Public interface to returning the analytic solution for conserved vars
tk::InitializeFn::result_type
solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
{ return self->solution( xi, yi, zi, t ); }
//! Copy assignment
DGPDE& operator=( const DGPDE& x )
{ DGPDE tmp(x); *this = std::move(tmp); return *this; }
//! Copy constructor
DGPDE( const DGPDE& x ) : self( x.self->copy() ) {}
//! Move assignment
DGPDE& operator=( DGPDE&& ) noexcept = default;
//! Move constructor
DGPDE( DGPDE&& ) noexcept = default;
private:
//! \brief Concept is a pure virtual base class specifying the requirements
//! of polymorphic objects deriving from it
struct Concept {
Concept() = default;
Concept( const Concept& ) = default;
virtual ~Concept() = default;
virtual Concept* copy() const = 0;
virtual std::size_t nprim() const = 0;
virtual std::size_t nmat() const = 0;
virtual void numEquationDofs(std::vector< std::size_t >&) const = 0;
virtual void IcBoxElems( const tk::Fields&,
std::size_t,
std::vector< std::unordered_set< std::size_t > >& ) const = 0;
virtual void initialize(
const tk::Fields&,
const std::vector< std::size_t >&,
const tk::UnsMesh::Coords&,
const std::vector< std::unordered_set< std::size_t > >&,
tk::Fields&,
tk::real,
const std::size_t nielem ) const = 0;
virtual void lhs( const tk::Fields&, tk::Fields& ) const = 0;
virtual void updateInterfaceCells( tk::Fields&,
std::size_t,
std::vector< std::size_t >& ) const = 0;
virtual void updatePrimitives( const tk::Fields&,
const tk::Fields&,
const tk::Fields&,
tk::Fields&,
std::size_t ) const = 0;
virtual void cleanTraceMaterial( const tk::Fields&,
tk::Fields&,
tk::Fields&,
std::size_t ) const = 0;
virtual void reconstruct( tk::real,
const tk::Fields&,
const tk::Fields&,
const inciter::FaceData&,
const std::map< std::size_t,
std::vector< std::size_t > >&,
const std::vector< std::size_t >&,
const tk::UnsMesh::Coords&,
tk::Fields&,
tk::Fields& ) const = 0;
virtual void limit( tk::real,
const tk::Fields&,
const tk::Fields&,
const inciter::FaceData&,
const std::map< std::size_t,
std::vector< std::size_t > >&,
const std::vector< std::size_t >&,
const tk::UnsMesh::Coords&,
const std::vector< std::size_t >&,
const std::vector< std::size_t >&,
const std::unordered_map< std::size_t, std::size_t >&,
const std::vector< std::vector<tk::real> >&,
const std::vector< std::vector<tk::real> >&,
tk::Fields&,
tk::Fields&,
std::vector< std::size_t >& ) const = 0;
virtual void rhs( tk::real,
const tk::Fields&,
const tk::Fields&,
const inciter::FaceData&,
const std::vector< std::size_t >&,
const std::vector< std::unordered_set< std::size_t > >&,
const tk::UnsMesh::Coords&,
const tk::Fields&,
const tk::Fields&,
const std::vector< std::size_t >&,
tk::Fields& ) const = 0;
virtual void eval_ndof( std::size_t,
const tk::UnsMesh::Coords&,
const std::vector< std::size_t >&,
const inciter::FaceData&,
const tk::Fields&,
inciter::ctr::PrefIndicatorType,
std::size_t,
std::size_t,
tk::real,
std::vector< std::size_t >& ) const = 0;
virtual tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
const std::vector< std::size_t >&,
const inciter::FaceData&,
const tk::Fields&,
const tk::Fields&,
const std::vector< std::size_t >&,
const tk::Fields&,
const tk::Fields&,
const std::size_t ) const = 0;
virtual std::vector< std::string > analyticFieldNames() const = 0;
virtual std::vector< std::string > histNames() const = 0;
virtual std::vector< std::string > names() const = 0;
virtual std::vector< std::vector< tk::real > > surfOutput(
const std::map< int, std::vector< std::size_t > >&,
tk::Fields& ) const = 0;
virtual std::vector< std::vector< tk::real > > histOutput(
const std::vector< HistData >&,
const std::vector< std::size_t >&,
const tk::UnsMesh::Coords&,
const tk::Fields&,
const tk::Fields& ) const = 0;
virtual tk::InitializeFn::result_type analyticSolution(
tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
virtual tk::InitializeFn::result_type solution(
tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
};
//! \brief Model models the Concept above by deriving from it and overriding
//! the virtual functions required by Concept
template< typename T >
struct Model : Concept {
explicit Model( T x ) : data( std::move(x) ) {}
Concept* copy() const override { return new Model( *this ); }
std::size_t nprim() const override
{ return data.nprim(); }
std::size_t nmat() const override
{ return data.nmat(); }
void numEquationDofs(std::vector< std::size_t >& numEqDof) const override
{ data.numEquationDofs(numEqDof); }
void IcBoxElems( const tk::Fields& geoElem,
std::size_t nielem,
std::vector< std::unordered_set< std::size_t > >& inbox )
const override { data.IcBoxElems( geoElem, nielem, inbox ); }
void initialize(
const tk::Fields& L,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
const std::vector< std::unordered_set< std::size_t > >& inbox,
tk::Fields& unk,
tk::real t,
const std::size_t nielem )
const override { data.initialize( L, inpoel, coord, inbox, unk, t,
nielem ); }
void lhs( const tk::Fields& geoElem, tk::Fields& l ) const override
{ data.lhs( geoElem, l ); }
void updateInterfaceCells( tk::Fields& unk,
std::size_t nielem,
std::vector< std::size_t >& ndofel )
const override { data.updateInterfaceCells( unk, nielem, ndofel ); }
void updatePrimitives( const tk::Fields& unk,
const tk::Fields& L,
const tk::Fields& geoElem,
tk::Fields& prim,
std::size_t nielem )
const override { data.updatePrimitives( unk, L, geoElem, prim, nielem ); }
void cleanTraceMaterial( const tk::Fields& geoElem,
tk::Fields& unk,
tk::Fields& prim,
std::size_t nielem )
const override { data.cleanTraceMaterial( geoElem, unk, prim, nielem ); }
void reconstruct( tk::real t,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const inciter::FaceData& fd,
const std::map< std::size_t,
std::vector< std::size_t > >& esup,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
tk::Fields& U,
tk::Fields& P ) const override
{
data.reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P );
}
void limit( tk::real t,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const inciter::FaceData& fd,
const std::map< std::size_t, std::vector< std::size_t > >&
esup,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
const std::vector< std::size_t >& ndofel,
const std::vector< std::size_t >& gid,
const std::unordered_map< std::size_t, std::size_t >& bid,
const std::vector< std::vector<tk::real> >& uNodalExtrm,
const std::vector< std::vector<tk::real> >& pNodalExtrm,
tk::Fields& U,
tk::Fields& P,
std::vector< std::size_t >& shockmarker ) const override
{
data.limit( t, geoFace, geoElem, fd, esup, inpoel, coord, ndofel, gid,
bid, uNodalExtrm, pNodalExtrm, U, P, shockmarker );
}
void rhs(
tk::real t,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const inciter::FaceData& fd,
const std::vector< std::size_t >& inpoel,
const std::vector< std::unordered_set< std::size_t > >& boxelems,
const tk::UnsMesh::Coords& coord,
const tk::Fields& U,
const tk::Fields& P,
const std::vector< std::size_t >& ndofel,
tk::Fields& R ) const override
{
data.rhs( t, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
ndofel, R );
}
void eval_ndof( std::size_t nunk,
const tk::UnsMesh::Coords& coord,
const std::vector< std::size_t >& inpoel,
const inciter::FaceData& fd,
const tk::Fields& unk,
inciter::ctr::PrefIndicatorType indicator,
std::size_t ndof,
std::size_t ndofmax,
tk::real tolref,
std::vector< std::size_t >& ndofel ) const override
{ data.eval_ndof( nunk, coord, inpoel, fd, unk, indicator, ndof, ndofmax,
tolref, ndofel ); }
tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
const std::vector< std::size_t >& inpoel,
const inciter::FaceData& fd,
const tk::Fields& geoFace,
const tk::Fields& geoElem,
const std::vector< std::size_t >& ndofel,
const tk::Fields& U,
const tk::Fields& P,
const std::size_t nielem ) const override
{ return data.dt( coord, inpoel, fd, geoFace, geoElem, ndofel, U, P,
nielem ); }
std::vector< std::string > analyticFieldNames() const override
{ return data.analyticFieldNames(); }
std::vector< std::string > histNames() const override
{ return data.histNames(); }
std::vector< std::string > names() const override
{ return data.names(); }
std::vector< std::vector< tk::real > > surfOutput(
const std::map< int, std::vector< std::size_t > >& bnd,
tk::Fields& U ) const override
{ return data.surfOutput( bnd, U ); }
std::vector< std::vector< tk::real > > histOutput(
const std::vector< HistData >& h,
const std::vector< std::size_t >& inpoel,
const tk::UnsMesh::Coords& coord,
const tk::Fields& U,
const tk::Fields& P ) const override
{ return data.histOutput( h, inpoel, coord, U, P ); }
tk::InitializeFn::result_type
analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
const override { return data.analyticSolution( xi, yi, zi, t ); }
tk::InitializeFn::result_type
solution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
const override { return data.solution( xi, yi, zi, t ); }
T data;
};
std::unique_ptr< Concept > self; //!< Base pointer used polymorphically
};
} // inciter::
#endif // DGPDE_h
|