Peano 4
Loading...
Searching...
No Matches
SelfSimilarInfallDG.cpp
Go to the documentation of this file.
4#include "EulerKernels.h"
5#include "spherical-accretion/GravityModel.h"
6
7
8tarch::logging::Log benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallDG::_log( "benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallDG" );
9
10
11
13
14
18
19
22 const tarch::la::Vector<3,double>& cellSize,
23 const tarch::la::Vector<3,int>& index,
24 double density
25) {
26 double overDensity = density - BaseDensity;
27 double radius = tarch::la::norm2( x );
28
30 ::exahype2::dg::getQuadratureWeight(cellSize,index,QuadratureWeights1d)>0.0,
31 cellSize, index, QuadratureWeights1d[0]
32 );
33 double mass = overDensity * ::exahype2::dg::getQuadratureWeight(cellSize,index,QuadratureWeights1d);
34
35 _accumulator.addMass(mass,radius);
36}
37
38
40 const double * __restrict__ Q, // Q[5+0]
43 double t
44) {
45 if (
47 and
48 tarch::la::norm2(x) <= InitialTopHatRadius + h(0)
49 ) {
50 return ::exahype2::RefinementCommand::Refine;
51 }
52 else if (
53 DynamicAMR
54 and
55 ( Q[4]>1.1*BaseDensity )
56 ) {
57 return ::exahype2::RefinementCommand::Refine;
58 }
59 else if (
60 DynamicAMR
61 ) {
62 return ::exahype2::RefinementCommand::Erase;
63 }
64 else {
65 return ::exahype2::RefinementCommand::Keep;
66 }
67}
68
69
70
72 double * __restrict__ Q,
76 bool gridIsConstructed
77) {
78 logTraceInWith2Arguments( "initialCondition(...)", x, gridIsConstructed );
79
80 for (int i=0; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
81 Q[i] = 0.0;
82 }
83
84/*
85 applications::exahype2::euler::sphericalaccretion::initialiseOverdensity_topHat(
86 Q,
87 x,
88 InitialTopHatRadius,
89 AdditionalMass,
90 BaseDensity,
91 pInitial,
92 Gamma
93 );
94*/
95
96// applications::exahype2::euler::sphericalaccretion::initialiseOverdensity_hyperbolicSecant(
97
98
100 Q,
101 x,
102 InitialTopHatRadius,
103 AdditionalMass,
104 BaseDensity,
105 pInitial,
106 Gamma
107 );
108
109 /*
110 applications::exahype2::euler::sphericalaccretion::initialiseHomogeneousDensity(
111 Q,
112 BaseDensity,
113 pInitial,
114 Gamma
115 );
116*/
117
118 logTraceOut( "initialCondition(...)" );
119}
120
121
123 const double * __restrict__ Qinside, // Qinside[5+0]
124 double * __restrict__ Qoutside, // Qoutside[5+0]
126 double t,
127 int normal
128) {
129 logTraceInWith3Arguments( "boundaryConditions(...)", x, t, normal );
130
131 ::applications::exahype2::euler::NeumannBoundaryConditions(
132 Qinside,
133 Qoutside,
134 x,
135 0.0, // no h in classical sense
136 t,
137 normal,
138 NumberOfUnknowns, // defined in superclass
139 NumberOfAuxiliaryVariables // defined in superclass
140 );
141
142 logTraceOut( "boundaryConditions(...)" );
143}
144
145
146double ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallDG::maxEigenvalue(
147 const double * __restrict__ Q, // Q[5+0],
149 double t,
150 double dt,
151 int normal
152) {
153 logTraceInWith4Arguments( "maxEigenvalue(...)", x, t, dt, normal );
154
156 Q,
157 x,
158 0.0, // no volume, as point-wise
159 t,
160 dt,
161 normal,
162 NumberOfUnknowns,
163 NumberOfAuxiliaryVariables,
164 Gamma
165 );
166
167 logTraceOutWith1Argument( "maxEigenvalue(...)", result );
168
169 return result;
170}
171
172
173void ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallDG::flux(
174 const double * __restrict__ Q, // Q[5+0],
176 double t,
177 double dt,
178 int normal,
179 double * __restrict__ F // F[5]
180) {
181 logTraceInWith4Arguments( "flux(...)", x, t, dt, normal );
182
184 Q,
185 x,
186 0.0, // no volume, as point-wise
187 t,
188 dt,
189 normal,
190 NumberOfUnknowns,
191 NumberOfAuxiliaryVariables,
192 F,
193 Gamma
194 );
195
196 logTraceOut( "flux(...)" );
197}
198
199
200void ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallDG::sourceTerm(
201 const double * __restrict__ Q, // Q[5+0]
203 double t,
204 double dt,
205 double * __restrict__ S // S[5
206) {
207 logTraceInWith3Arguments( "sourceTerm(...)", x, t, dt );
208
209 for (int i=0; i<NumberOfUnknowns; i++) {
210 S[i] = 0.0;
211 }
212
213 double massInterpolated = _accumulator.getMass_linearInterpolation( tarch::la::norm2( x ) );
214
216 S,
217 x,
218 Q,
219 massInterpolated,
220 aInitial,
221 t
222 );
223
224 logTraceOut( "sourceTerm(...)" );
225}
226
227
228void ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallDG::finishTimeStep() {
229 AbstractSelfSimilarInfallDG::finishTimeStep();
230 if (isLastGridSweepOfTimeStep()) {
231 _accumulator.finishAccumulation();
232 TotalMassInPreviousTimeStep = _accumulator.getTotalMass();
233 }
234}
235
236
237bool ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallDG::cellCanUseStatelessPDETerms(
238 const tarch::la::Vector<Dimensions,double>& cellCentre,
240 double t,
241 double dt
242) const {
243 const double radius = tarch::la::norm2( cellCentre );
244 return radius > _accumulator.getMaxRelevantRadius() + std::sqrt(Dimensions) * cellH(0) * cellH(0);
245
246}
#define assertion3(expr, param0, param1, param2)
#define logTraceOutWith1Argument(methodName, argument0)
Definition Log.h:380
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:373
#define logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
Definition Log.h:372
#define logTraceInWith2Arguments(methodName, argument0, argument1)
Definition Log.h:371
void initialCondition(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, const tarch::la::Vector< Dimensions, double > &h, const tarch::la::Vector< Dimensions, int > &index, bool gridIsConstructed) override
static double TotalMassInPreviousTimeStep
After each time step, we backup the total mass in this static variable.
virtual void boundaryConditions(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) override
::exahype2::RefinementCommand refinementCriterion(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, const tarch::la::Vector< Dimensions, double > &h, double t) override
Evaluate refinement criterion over cell.
void addDensity(const tarch::la::Vector< Dimensions, double > &x, const tarch::la::Vector< 3, double > &cellSize, const tarch::la::Vector< 3, int > &index, double density)
Add mass corresponding to one degree of freedom.
Log Device.
Definition Log.h:516
void addGravitationalSource_AlphaCDM(double *__restrict__ S, const tarch::la::Vector< Dimensions, double > &x, const double *__restrict__ Q, double mass, double aInitial, double t)
void initialiseOverdensity_Gaussian(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double topHatRadius, double additionalMass, double baselineDensity, double baselinePressure, double gamma)
Initialise overdensity according to Gaussian.
static void flux(const double *__restrict__ Q, int normal, double *__restrict__ F) InlineMethod
static double maxEigenvalue(const double *__restrict__ Q, int normal) InlineMethod
double getQuadratureWeight(const tarch::la::Vector< 3, double > &cellSize, const tarch::la::Vector< 3, int > &index, const double *__restrict__ quadratureWeights)
Compute integral over shape function over cell defined by index.
Definition DGUtils.cpp:25
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
Simple vector class.
Definition Vector.h:134