Line data Source code
1 : #include "crpropa/advectionField/AdvectionField.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Common.h"
4 :
5 : #include "gtest/gtest.h"
6 : #include <stdexcept>
7 : #include <cmath>
8 :
9 : namespace crpropa {
10 :
11 2 : TEST(testUniformAdvectionField, SimpleTest) {
12 1 : UniformAdvectionField A(Vector3d(-1, 5, 3));
13 1 : Vector3d a = A.getField(Vector3d(1, 0, 0));
14 1 : double D = A.getDivergence(Vector3d(1, 0, 0));
15 1 : EXPECT_DOUBLE_EQ(a.x, -1);
16 1 : EXPECT_DOUBLE_EQ(a.y, 5);
17 1 : EXPECT_DOUBLE_EQ(a.z, 3);
18 1 : EXPECT_DOUBLE_EQ(D, 0.);
19 1 : }
20 :
21 2 : TEST(testAdvectionFieldList, SimpleTest) {
22 : // Test a list of three advection fields
23 : AdvectionFieldList A;
24 2 : A.addField(new UniformAdvectionField(Vector3d(1, 0, 0)));
25 2 : A.addField(new UniformAdvectionField(Vector3d(0, 2, 0)));
26 2 : A.addField(new UniformAdvectionField(Vector3d(0, 0, 3)));
27 1 : Vector3d a = A.getField(Vector3d(0.));
28 1 : double D = A.getDivergence(Vector3d(1, 2, 3));
29 1 : EXPECT_DOUBLE_EQ(a.x, 1);
30 1 : EXPECT_DOUBLE_EQ(a.y, 2);
31 1 : EXPECT_DOUBLE_EQ(a.z, 3);
32 1 : EXPECT_DOUBLE_EQ(D, 0.);
33 1 : }
34 :
35 2 : TEST(testConstantSphericalAdvectionField, SimpleTest) {
36 :
37 : Vector3d origin(1, 0, 0);
38 : double V_wind(10);
39 :
40 1 : ConstantSphericalAdvectionField A(origin, V_wind);
41 :
42 : // Check the properties of the advection field
43 1 : EXPECT_DOUBLE_EQ(A.getOrigin().x, origin.x);
44 1 : EXPECT_DOUBLE_EQ(A.getOrigin().y, origin.y);
45 1 : EXPECT_DOUBLE_EQ(A.getOrigin().z, origin.z);
46 :
47 1 : EXPECT_DOUBLE_EQ(A.getVWind(), V_wind);
48 :
49 : // Field should be radial with center (1,0,0)
50 : Vector3d Pos(2, 1, 1);
51 1 : Vector3d V = A.getField(Pos);
52 :
53 1 : EXPECT_DOUBLE_EQ(V.x, 10.*pow(3, -0.5));
54 1 : EXPECT_DOUBLE_EQ(V.y, 10.*pow(3, -0.5));
55 1 : EXPECT_DOUBLE_EQ(V.z, 10.*pow(3, -0.5));
56 :
57 : // Divergence should be 2*V/r
58 : Vector3d Pos2(2, 0, 0);
59 1 : double D = A.getDivergence(Pos2);
60 1 : EXPECT_DOUBLE_EQ(D, 2*10);
61 1 : }
62 :
63 2 : TEST(testSphericalAdvectionField, SimpleTest) {
64 :
65 : Vector3d origin(1, 0, 0);
66 : double R_max(10);
67 : double V_max(1000);
68 : double tau(3.);
69 : double alpha(2.);
70 :
71 1 : SphericalAdvectionField A(origin, R_max, V_max, tau, alpha);
72 :
73 : // Check the properties of the advection field
74 1 : EXPECT_DOUBLE_EQ(A.getOrigin().x, origin.x);
75 1 : EXPECT_DOUBLE_EQ(A.getOrigin().y, origin.y);
76 1 : EXPECT_DOUBLE_EQ(A.getOrigin().z, origin.z);
77 :
78 1 : EXPECT_DOUBLE_EQ(A.getRadius(), R_max);
79 1 : EXPECT_DOUBLE_EQ(A.getVMax(), V_max);
80 1 : EXPECT_DOUBLE_EQ(A.getTau(), tau);
81 1 : EXPECT_DOUBLE_EQ(A.getAlpha(), alpha);
82 :
83 : // Field/divergence should be zero for R>10
84 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(12,0,0)).getR(), 0.);
85 1 : EXPECT_DOUBLE_EQ(A.getDivergence(Vector3d(12,0,0)), 0.);
86 :
87 : // Check field and divergence
88 : Vector3d Pos(2, 0, 0);
89 1 : Vector3d a = A.getField(Pos);
90 1 : Vector3d a0 = a.getUnitVector();
91 1 : double d = A.getDivergence(Pos);
92 :
93 1 : EXPECT_DOUBLE_EQ(a0.x, 1.);
94 1 : EXPECT_DOUBLE_EQ(a0.y, 0.);
95 1 : EXPECT_DOUBLE_EQ(a0.z, 0.);
96 1 : EXPECT_DOUBLE_EQ(a.getR(), V_max*(1.-exp(-1./tau)) );
97 :
98 1 : EXPECT_NEAR(d, 1044.624919, 1e-5);
99 :
100 : // Check asymptotic of the Field
101 1 : EXPECT_NEAR(A.getField(Vector3d(11, 0, 0)).getR(), 1000., 1e-4);
102 :
103 : // Divergence should be 2*V_max/r for r>>1
104 1 : EXPECT_NEAR(A.getDivergence(Vector3d(11, 0, 0)), 2*1000./10., 1e-4);
105 1 : }
106 :
107 :
108 2 : TEST(testSphericalAdvectionShock, SimpleTest) {
109 :
110 : Vector3d origin(0, 0, 0);
111 : double R_0(10);
112 : double V_0(1000);
113 : double lambda(0.1);
114 : double R_rot(1.);
115 : double V_rot(100);
116 :
117 :
118 1 : SphericalAdvectionShock A(origin, R_0, V_0, lambda);
119 :
120 : // Check the properties of the advection field
121 1 : EXPECT_DOUBLE_EQ(A.getOrigin().x, origin.x);
122 1 : EXPECT_DOUBLE_EQ(A.getOrigin().y, origin.y);
123 1 : EXPECT_DOUBLE_EQ(A.getOrigin().z, origin.z);
124 :
125 1 : EXPECT_DOUBLE_EQ(A.getR0(), R_0);
126 1 : EXPECT_DOUBLE_EQ(A.getV0(), V_0);
127 1 : EXPECT_DOUBLE_EQ(A.getLambda(), lambda);
128 :
129 : // Default values for R_Rot is R_0
130 : // Default value for Azimuthal speed is 0
131 1 : EXPECT_DOUBLE_EQ(A.getRRot(), R_0);
132 1 : EXPECT_DOUBLE_EQ(A.getAzimuthalSpeed(), 0.);
133 :
134 : // Field should drop to 0.625*V_0 at the shock
135 : // That's a difference to the analytic shock model where it should
136 : // drop to v_r(R_0)=0.25*V_0.
137 1 : EXPECT_DOUBLE_EQ(A.getField(Vector3d(R_0,0,0)).getR(), 0.625*V_0);
138 :
139 : // Field divergence should be zero for R>>R_0
140 1 : EXPECT_NEAR(A.getDivergence(Vector3d(15,0,0)), 0., 1e-10);
141 :
142 : // Check field and divergence
143 : Vector3d Pos(2, 0, 0);
144 1 : Vector3d a = A.getField(Pos);
145 1 : Vector3d a0 = a.getUnitVector();
146 1 : double d = A.getDivergence(Pos);
147 :
148 1 : EXPECT_DOUBLE_EQ(a0.x, 1.);
149 1 : EXPECT_DOUBLE_EQ(a0.y, 0.);
150 1 : EXPECT_DOUBLE_EQ(a0.z, 0.);
151 1 : EXPECT_DOUBLE_EQ(a.getR(), V_0);
152 :
153 : //Set explicitely the azimuthal rotation speed
154 1 : A.setRRot(R_rot);
155 1 : A.setAzimuthalSpeed(V_rot);
156 :
157 1 : EXPECT_DOUBLE_EQ(A.getRRot(), R_rot);
158 1 : EXPECT_DOUBLE_EQ(A.getAzimuthalSpeed(), V_rot);
159 :
160 : Vector3d pos(1., 0., 0.);
161 1 : Vector3d f = A.getField(pos);
162 1 : EXPECT_DOUBLE_EQ(f.x, 1000.);
163 1 : EXPECT_DOUBLE_EQ(f.y, 100.);
164 1 : EXPECT_DOUBLE_EQ(f.z, 0.);
165 :
166 :
167 1 : }
168 :
169 : } //namespace crpropa
|
