Peano 4
Loading...
Searching...
No Matches
SelfSimilarInfallFV.cpp
Go to the documentation of this file.
3
4tarch::logging::Log benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallFV::_log( "benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallFV" );
5
6
7
9
10
14
15
17 const tarch::la::Vector<Dimensions,double>& volumeCentre,
19 double density
20) {
21 double overDensity = density - BaseDensity;
22 double radius = tarch::la::norm2( volumeCentre );
23
24 double mass = overDensity * tarch::la::volume( volumeH );
25
26 logDebug( "touchCellFirstTime(...)", "dump data of " << volumeCentre << ": " << density << " x " << volumeH << "->" << radius << " x " << mass );
27
28 _accumulator.addMass(mass,radius);
29}
30
31
33 const double * __restrict__ Q, // Q[5+0],
34 const tarch::la::Vector<Dimensions,double>& volumeCentre,
36 double t
37) {
38 if (
40 and
41 tarch::la::norm2(volumeCentre) <= InitialTopHatRadius + volumeH(0) * NumberOfFiniteVolumesPerAxisPerPatch
42 ) {
43 return ::exahype2::RefinementCommand::Refine;
44 }
45 else if (
46 DynamicAMR
47 and
48 ( Q[4]>1.1*BaseDensity )
49 ) {
50 return ::exahype2::RefinementCommand::Refine;
51 }
52 else if (
53 DynamicAMR
54 ) {
55 return ::exahype2::RefinementCommand::Erase;
56 }
57 else {
58 return ::exahype2::RefinementCommand::Keep;
59 }
60}
61
62
64 double * __restrict__ Q,
67 bool gridIsConstructed
68) {
69 logTraceInWith3Arguments( "initialCondition(...)", volumeX, volumeH, gridIsConstructed );
70
71 for (int i=0; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
72 Q[i] = 0.0;
73 }
74
76 Q,
77 volumeX,
78 InitialTopHatRadius,
79 AdditionalMass,
80 BaseDensity,
81 pInitial,
82 Gamma
83 );
84
85 logDebug( "initialCondition(...)", volumeX << " (" << tarch::la::norm2(volumeX) << ") -> " << Q[0] );
86
87/*
88 applications::exahype2::euler::sphericalaccretion::initialiseHomogeneousDensity(
89 Q,
90 BaseDensity,
91 pInitial,
92 Gamma
93 );
94*/
95
96 logTraceOut( "initialCondition(...)" );
97}
98
99
101 const double * __restrict__ Qinside, // Qinside[5+0]
102 double * __restrict__ Qoutside, // Qoutside[5+0]
103 const tarch::la::Vector<Dimensions,double>& faceCentre,
105 double t,
106 int normal
107) {
108 logTraceInWith4Arguments( "boundaryConditions(...)", faceCentre, volumeH, t, normal );
109
110 ::applications::exahype2::euler::NeumannBoundaryConditions(
111 Qinside,
112 Qoutside,
113 faceCentre,
114 volumeH,
115 t,
116 normal,
117 NumberOfUnknowns, // defined in superclass
118 NumberOfAuxiliaryVariables // defined in superclass
119 );
120
121 logTraceOut( "boundaryConditions(...)" );
122}
123
124
125double ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallFV::maxEigenvalue(
126 const double * __restrict__ Q,
127 const tarch::la::Vector<Dimensions,double>& faceCentre,
129 double t,
130 double dt,
131 int normal
132) {
133 logTraceInWith4Arguments( "maxEigenvalue(...)", faceCentre, volumeH, t, normal );
134
136 Q,
137 faceCentre,
138 volumeH,
139 t,
140 dt,
141 normal,
142 NumberOfUnknowns,
143 NumberOfAuxiliaryVariables,
144 Gamma
145 );
146
147 logTraceOutWith1Argument( "maxEigenvalue(...)", result );
148
149 return result;
150}
151
152
153void ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallFV::flux(
154 const double * __restrict__ Q,
155 const tarch::la::Vector<Dimensions,double>& faceCentre,
157 double t,
158 double dt,
159 int normal,
160 double * __restrict__ F // F[5]
161) {
162 logTraceInWith4Arguments( "flux(...)", faceCentre, volumeH, t, normal );
163
165 Q,
166 faceCentre,
167 volumeH,
168 t,
169 dt,
170 normal,
171 NumberOfUnknowns,
172 NumberOfAuxiliaryVariables,
173 F,
174 Gamma
175 );
176
177 logTraceOut( "flux(...)" );
178}
179
180
181void ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallFV::sourceTerm(
182 const double * __restrict__ Q, // Q[5+0]
185 double t,
186 double dt,
187 double * __restrict__ S // S[5
188) {
189 logTraceInWith4Arguments( "sourceTerm(...)", volumeX, volumeH, t, dt );
190
191 for (int i=0; i<NumberOfUnknowns; i++) {
192 S[i] = 0.0;
193 }
194
195 double massInterpolated = _accumulator.getMass_linearInterpolation( tarch::la::norm2( volumeX ) );
196
198 S,
199 volumeX,
200 Q,
201 massInterpolated,
202 aInitial,
203 t
204 );
205
206 logTraceOut( "sourceTerm(...)" );
207}
208
209
210void ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallFV::finishTimeStep() {
211 AbstractSelfSimilarInfallFV::finishTimeStep();
212 if (isLastGridSweepOfTimeStep()) {
213 _accumulator.finishAccumulation();
214 TotalMassInPreviousTimeStep = _accumulator.getTotalMass();
215 }
216}
217
218
219bool ::benchmarks::exahype2::euler::sphericalaccretionupscaling::SelfSimilarInfallFV::patchCanUseStatelessPDETerms(
220 const tarch::la::Vector<Dimensions,double>& patchCentre,
222 double t,
223 double dt
224) const {
225 const double radius = tarch::la::norm2( patchCentre );
226 return radius > _accumulator.getMaxRelevantRadius() + std::sqrt(Dimensions) * patchH(0) * patchH(0);
227}
228
#define logDebug(methodName, logMacroMessageStream)
Definition Log.h:50
#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
void addDensity(const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double density)
Add density of one voxel to global bookkeeping.
static double TotalMassInPreviousTimeStep
After each time step, we backup the total mass in this static variable.
void initialCondition(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, bool gridIsConstructed) override
We start with a top hat overdensity.
::exahype2::RefinementCommand refinementCriterion(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t) override
virtual void boundaryConditions(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal) override
Apply Neumann boundary conditions.
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_topHat(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double topHatRadius, double additionalMass, double baselineDensity, double baselinePressure, double gamma)
Set the initial overdensity profile as top hat.
static void flux(const double *__restrict__ Q, int normal, double *__restrict__ F) InlineMethod
static double maxEigenvalue(const double *__restrict__ Q, int normal) InlineMethod
Scalar volume(const Vector< Size, Scalar > &vector)
Computes the volume of the tetrahedron spanned by the Cartesian unit vectors scaled by the correspond...
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