Peano 4
Loading...
Searching...
No Matches
GravityModel.cpp
Go to the documentation of this file.
2
3
5 double* __restrict__ Q, double baselineDensity, double baselinePressure, double gamma
6) {
7 Q[0] = baselineDensity;
8 Q[1] = 0.0;
9 Q[2] = 0.0;
10 Q[3] = 0.0;
11 Q[4] = 0.0;
12}
13
14
16 double* __restrict__ Q,
18 double topHatRadius,
19 double additionalMass,
20 double baselineDensity,
21 double baselinePressure,
22 double gamma
23) {
24 double radius = tarch::la::norm2(x);
25
26 double additionalDensity = (radius < topHatRadius)
27 ? additionalMass * 3.0 / 4.0 / tarch::la::PI / std::pow(topHatRadius, 3.0)
28 : 0.0;
29
30
31 Q[0] = baselineDensity + additionalDensity;
32 Q[1] = 0.0;
33 Q[2] = 0.0;
34 Q[3] = 0.0;
35 Q[4] = 0.5 * baselinePressure / (gamma - 1);
36}
37
38
40 double* __restrict__ Q,
42 double topHatRadius,
43 double additionalMass,
44 double baselineDensity,
45 double baselinePressure,
46 double gamma
47) {
48 double standardDeviation = topHatRadius;
49 double gaussianScaling = 1.0;
50
51 double radius = tarch::la::norm2(x);
52
53 double exponent = -0.5 * (radius / standardDeviation) * (radius / standardDeviation);
54
55 double gaussian = 1.0 / (standardDeviation * std::sqrt(2.0 * tarch::la::PI)) * std::exp(exponent);
56
57 double additionalDensity = additionalMass * gaussian * gaussianScaling;
58
59 Q[0] = baselineDensity + additionalDensity;
60 Q[1] = 0.0;
61 Q[2] = 0.0;
62 Q[3] = 0.0;
63 Q[4] = 0.5 * baselinePressure / (gamma - 1);
64}
65
66
68 double* __restrict__ Q,
70 double topHatRadius,
71 double additionalMass,
72 double baselineDensity,
73 double baselinePressure,
74 double gamma
75) {
76 double radius = tarch::la::norm2(x);
77
78 double topHatAdditionalDensity = additionalMass * 3.0 / 4.0 / tarch::la::PI / std::pow(topHatRadius, 3.0);
79
80 double scalingOfTopHat = 1.2;
81 double rNormalised = radius / (scalingOfTopHat * topHatRadius);
82
83 double scaling = tarch::la::smallerEquals(rNormalised, 1.0) ? std::exp(-1.0 / (1 - rNormalised * rNormalised)) : 0.0;
84
85 double additionalDensity = topHatAdditionalDensity * scaling;
86
87 Q[0] = baselineDensity + additionalDensity;
88 Q[1] = 0.0;
89 Q[2] = 0.0;
90 Q[3] = 0.0;
91 Q[4] = 0.5 * baselinePressure / (gamma - 1);
92}
93
94
96 double* __restrict__ Q,
98 double topHatRadius,
99 double additionalMass,
100 double baselineDensity,
101 double baselinePressure,
102 double gamma
103) {
104 double radius = tarch::la::norm2(x);
105
106 double topHatAdditionalDensity = additionalMass * 3.0 / 4.0 / tarch::la::PI / std::pow(topHatRadius, 3.0);
107
108 double scaledRadius = radius / topHatRadius;
109 double volumetricScaling = 1e-1;
110
111 double additionalDensity = topHatAdditionalDensity * 2.0 / (std::exp(scaledRadius) + std::exp(-scaledRadius))
112 * volumetricScaling;
113
114 assertion(additionalDensity >= 0.0);
115
116 Q[0] = baselineDensity + additionalDensity;
117 Q[1] = 0.0;
118 Q[2] = 0.0;
119 Q[3] = 0.0;
120 Q[4] = 0.5 * baselinePressure / (gamma - 1);
121}
122
123
125 double* __restrict__ S,
127 const double* __restrict__ Q,
128 double mass,
129 double aInitial,
130 double t
131) {
132 double radius = tarch::la::norm2(x);
133
134 if (tarch::la::greater(radius, 0.0)) {
135 double a = 1e-3 * 0.0287 * std::pow((-t / 11.8 + 0.1694 * std::pow(aInitial, -0.5)), -2); // when code time ~
136 // 2*(a_i^(-0.5)-1), a~1
137
138 double normalisedDistance = std::max(1.0, radius);
139
140 double forceDensity = a * Q[0] * mass / std::pow(normalisedDistance, 3.0);
141
142 S[1] += -forceDensity * x(0) / radius; // velocities
143 S[2] += -forceDensity * x(1) / radius; // velocities
144 S[3] += -forceDensity * x(2) / radius; // velocities
145 S[4] += -forceDensity * (Q[1] * x(0) + Q[2] * x(1) + Q[3] * x(2)) / radius / Q[0];
146 }
147}
#define assertion(expr)
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_hyperbolicSecant(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double topHatRadius, double additionalMass, double baselineDensity, double baselinePressure, double gamma)
Alternative to Gaussian.
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.
void initialiseHomogeneousDensity(double *__restrict__ Q, double baselineDensity, double baselinePressure, double gamma)
Baseline configuration where no overdensity is applied.
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.
void initialiseOverdensity_bumpFunction(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double topHatRadius, double additionalMass, double baselineDensity, double baselinePressure, double gamma)
Alternative to Gaussian.
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
constexpr double PI
Definition Scalar.h:12
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
Simple vector class.
Definition Vector.h:134