Line data Source code
1 : #include "crpropa/advectionField/AdvectionField.h"
2 : #include "crpropa/advectionField/TimeDependentAdvectionField.h"
3 : #include "crpropa/Units.h"
4 : #include "crpropa/Common.h"
5 :
6 : #include "gtest/gtest.h"
7 : #include <stdexcept>
8 : #include <cmath>
9 :
10 : namespace crpropa {
11 :
12 2 : TEST(testUniformAdvectionField, SimpleTest) {
13 1 : UniformAdvectionField A(Vector3d(-1, 5, 3));
14 1 : Vector3d a = A.getField(Vector3d(1, 0, 0));
15 1 : double D = A.getDivergence(Vector3d(1, 0, 0));
16 1 : EXPECT_DOUBLE_EQ(a.x, -1);
17 1 : EXPECT_DOUBLE_EQ(a.y, 5);
18 1 : EXPECT_DOUBLE_EQ(a.z, 3);
19 1 : EXPECT_DOUBLE_EQ(D, 0.);
20 1 : }
21 :
22 2 : TEST(testAdvectionFieldList, SimpleTest) {
23 : // Test a list of three advection fields
24 : AdvectionFieldList A;
25 1 : A.addField(new UniformAdvectionField(Vector3d(1, 0, 0)));
26 1 : A.addField(new UniformAdvectionField(Vector3d(0, 2, 0)));
27 1 : A.addField(new UniformAdvectionField(Vector3d(0, 0, 3)));
28 1 : Vector3d a = A.getField(Vector3d(0.));
29 1 : double D = A.getDivergence(Vector3d(1, 2, 3));
30 1 : EXPECT_DOUBLE_EQ(a.x, 1);
31 1 : EXPECT_DOUBLE_EQ(a.y, 2);
32 1 : EXPECT_DOUBLE_EQ(a.z, 3);
33 1 : EXPECT_DOUBLE_EQ(D, 0.);
34 1 : }
35 :
36 2 : TEST(testConstantSphericalAdvectionField, SimpleTest) {
37 :
38 : Vector3d origin(1, 0, 0);
39 : double V_wind(10);
40 :
41 1 : ConstantSphericalAdvectionField A(origin, V_wind);
42 :
43 : // Check the properties of the advection field
44 1 : EXPECT_DOUBLE_EQ(A.getOrigin().x, origin.x);
45 1 : EXPECT_DOUBLE_EQ(A.getOrigin().y, origin.y);
46 1 : EXPECT_DOUBLE_EQ(A.getOrigin().z, origin.z);
47 :
48 1 : EXPECT_DOUBLE_EQ(A.getVWind(), V_wind);
49 :
50 : // Field should be radial with center (1,0,0)
51 : Vector3d Pos(2, 1, 1);
52 1 : Vector3d V = A.getField(Pos);
53 :
54 1 : EXPECT_DOUBLE_EQ(V.x, 10.*pow(3, -0.5));
55 1 : EXPECT_DOUBLE_EQ(V.y, 10.*pow(3, -0.5));
56 1 : EXPECT_DOUBLE_EQ(V.z, 10.*pow(3, -0.5));
57 :
58 : // Divergence should be 2*V/r
59 : Vector3d Pos2(2, 0, 0);
60 1 : double D = A.getDivergence(Pos2);
61 1 : EXPECT_DOUBLE_EQ(D, 2*10);
62 1 : }
63 :
64 2 : TEST(testSphericalAdvectionField, SimpleTest) {
65 :
66 : Vector3d origin(1, 0, 0);
67 : double R_max(10);
68 : double V_max(1000);
69 : double tau(3.);
70 : double alpha(2.);
71 :
72 1 : SphericalAdvectionField A(origin, R_max, V_max, tau, alpha);
73 :
74 : // Check the properties of the advection field
75 1 : EXPECT_DOUBLE_EQ(A.getOrigin().x, origin.x);
76 1 : EXPECT_DOUBLE_EQ(A.getOrigin().y, origin.y);
77 1 : EXPECT_DOUBLE_EQ(A.getOrigin().z, origin.z);
78 :
79 1 : EXPECT_DOUBLE_EQ(A.getRadius(), R_max);
80 1 : EXPECT_DOUBLE_EQ(A.getVMax(), V_max);
81 1 : EXPECT_DOUBLE_EQ(A.getTau(), tau);
82 1 : EXPECT_DOUBLE_EQ(A.getAlpha(), alpha);
83 :
84 : // Field/divergence should be zero for R>10
85 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(12,0,0)).getR(), 0.);
86 1 : EXPECT_DOUBLE_EQ(A.getDivergence(Vector3d(12,0,0)), 0.);
87 :
88 : // Check field and divergence
89 : Vector3d Pos(2, 0, 0);
90 1 : Vector3d a = A.getField(Pos);
91 1 : Vector3d a0 = a.getUnitVector();
92 1 : double d = A.getDivergence(Pos);
93 :
94 1 : EXPECT_DOUBLE_EQ(a0.x, 1.);
95 1 : EXPECT_DOUBLE_EQ(a0.y, 0.);
96 1 : EXPECT_DOUBLE_EQ(a0.z, 0.);
97 1 : EXPECT_DOUBLE_EQ(a.getR(), V_max*(1.-exp(-1./tau)) );
98 :
99 1 : EXPECT_NEAR(d, 1044.624919, 1e-5);
100 :
101 : // Check asymptotic of the Field
102 1 : EXPECT_NEAR(A.getField(Vector3d(11, 0, 0)).getR(), 1000., 1e-4);
103 :
104 : // Divergence should be 2*V_max/r for r>>1
105 1 : EXPECT_NEAR(A.getDivergence(Vector3d(11, 0, 0)), 2*1000./10., 1e-4);
106 1 : }
107 :
108 :
109 2 : TEST(testSphericalAdvectionShock, SimpleTest) {
110 :
111 : Vector3d origin(0, 0, 0);
112 : double R_0(10);
113 : double V_0(1000);
114 : double lambda(0.1);
115 : double R_rot(1.);
116 : double V_rot(100);
117 :
118 :
119 1 : SphericalAdvectionShock A(origin, R_0, V_0, lambda);
120 :
121 : // Check the properties of the advection field
122 1 : EXPECT_DOUBLE_EQ(A.getOrigin().x, origin.x);
123 1 : EXPECT_DOUBLE_EQ(A.getOrigin().y, origin.y);
124 1 : EXPECT_DOUBLE_EQ(A.getOrigin().z, origin.z);
125 :
126 1 : EXPECT_DOUBLE_EQ(A.getR0(), R_0);
127 1 : EXPECT_DOUBLE_EQ(A.getV0(), V_0);
128 1 : EXPECT_DOUBLE_EQ(A.getLambda(), lambda);
129 :
130 : // Default values for R_Rot is R_0
131 : // Default value for Azimuthal speed is 0
132 1 : EXPECT_DOUBLE_EQ(A.getRRot(), R_0);
133 1 : EXPECT_DOUBLE_EQ(A.getAzimuthalSpeed(), 0.);
134 :
135 : // Field should drop to 0.625*V_0 at the shock
136 : // That's a difference to the analytic shock model where it should
137 : // drop to v_r(R_0)=0.25*V_0.
138 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(R_0,0,0)).getR(), 0.625*V_0);
139 :
140 : // Field divergence should be zero for R>>R_0
141 1 : EXPECT_NEAR(A.getDivergence(Vector3d(15,0,0)), 0., 1e-10);
142 :
143 : // Check field and divergence
144 : Vector3d Pos(2, 0, 0);
145 1 : Vector3d a = A.getField(Pos);
146 1 : Vector3d a0 = a.getUnitVector();
147 1 : double d = A.getDivergence(Pos);
148 :
149 1 : EXPECT_DOUBLE_EQ(a0.x, 1.);
150 1 : EXPECT_DOUBLE_EQ(a0.y, 0.);
151 1 : EXPECT_DOUBLE_EQ(a0.z, 0.);
152 1 : EXPECT_DOUBLE_EQ(a.getR(), V_0);
153 :
154 : //Set explicitely the azimuthal rotation speed
155 1 : A.setRRot(R_rot);
156 1 : A.setAzimuthalSpeed(V_rot);
157 :
158 1 : EXPECT_DOUBLE_EQ(A.getRRot(), R_rot);
159 1 : EXPECT_DOUBLE_EQ(A.getAzimuthalSpeed(), V_rot);
160 :
161 : Vector3d pos(1., 0., 0.);
162 1 : Vector3d f = A.getField(pos);
163 1 : EXPECT_DOUBLE_EQ(f.x, 1000.);
164 1 : EXPECT_DOUBLE_EQ(f.y, 100.);
165 1 : EXPECT_DOUBLE_EQ(f.z, 0.);
166 :
167 :
168 1 : }
169 :
170 1 : TEST(testOneDimensionalTimeDependentShock, SimpleTest) {
171 :
172 : double V_sh(1);
173 : double V_1(3./4.);
174 : double V_0(0.);
175 : double L_sh(.01);
176 : double X_sh0(0);
177 : double T_sh0(0);
178 :
179 1 : OneDimensionalTimeDependentShock A(V_sh, V_1, V_0, L_sh);
180 :
181 : // Check the properties of the advection field
182 1 : EXPECT_DOUBLE_EQ(A.getVshock(), V_sh);
183 1 : EXPECT_DOUBLE_EQ(A.getV1(), V_1);
184 1 : EXPECT_DOUBLE_EQ(A.getV0(), V_0);
185 1 : EXPECT_DOUBLE_EQ(A.getShockWidth(), L_sh);
186 1 : EXPECT_DOUBLE_EQ(A.getShockPosition(0.), X_sh0);
187 1 : EXPECT_DOUBLE_EQ(A.getShockTime(), T_sh0);
188 :
189 : // Field should be zero for t=0
190 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(1,0,0)).getR(), 0.);
191 1 : EXPECT_DOUBLE_EQ(A.getDivergence(Vector3d(1,0,0)), 0.);
192 :
193 : // Shock position at t=1
194 1 : double xsh = A.getShockPosition(10.);
195 1 : EXPECT_DOUBLE_EQ(xsh, V_sh * 10.);
196 :
197 : // Preshock and postshock speeds:
198 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(xsh+5.,0,0), 10.).getR(), V_0);
199 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(xsh-5.,0,0), 10.).getR(), V_1);
200 :
201 1 : }
202 :
203 1 : TEST(testSedovTaylorBlastWave, SimpleTest) {
204 :
205 : double E0(1);
206 : double rho0(1);
207 : double L_sh(0.01);
208 :
209 1 : SedovTaylorBlastWave A(E0, rho0, L_sh);
210 :
211 : // Check the properties of the advection field
212 1 : EXPECT_DOUBLE_EQ(A.getEnergy(), E0);
213 1 : EXPECT_DOUBLE_EQ(A.getDensity(), rho0);
214 1 : EXPECT_DOUBLE_EQ(A.getShockWidth(), L_sh);
215 :
216 : // Field should be zero for t=0
217 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(1,0,0)).getR(), 0.);
218 1 : EXPECT_DOUBLE_EQ(A.getDivergence(Vector3d(1,0,0)), 0.);
219 :
220 : // Shock position at t=1
221 1 : double R = A.getShockRadius(1.);
222 1 : EXPECT_DOUBLE_EQ(R, pow(E0 / rho0, 1./5.));
223 :
224 : // Shock speed at t=1
225 1 : double V = A.getShockSpeed(1.);
226 1 : EXPECT_DOUBLE_EQ(V, 2. / 5. * pow(E0 / rho0, 1. / 5.));
227 :
228 : // Check field
229 : Vector3d Pos(2, 0, 0);
230 1 : Vector3d a = A.getField(Pos);
231 :
232 1 : EXPECT_DOUBLE_EQ(a.getR(), 0.5 * 2. * pow(2., 8) * 2. / 5. * pow(E0 / rho0, 1. / 5.) * (1 - tanh( (2 - 1) * pow(E0 / rho0, 1. / 5.) / L_sh) ));
233 :
234 : // Check asymptotic of the Field
235 1 : EXPECT_NEAR(A.getField(Vector3d(11, 0, 0)).getR(), 0, 1e-4);
236 1 : EXPECT_NEAR(A.getDivergence(Vector3d(11, 0, 0)), 0, 1e-4);
237 :
238 1 : }
239 :
240 : } //namespace crpropa
|
