71 #if PHYSICS == MHD || PHYSICS == RMHD
88 double dx1, x1, x1p, x1m;
89 double dx2, x2, x2p, x2m;
90 double dx3, x3, x3p, x3m;
91 double bx, by, bz, br, bphi, bth;
92 double us_p[256], us_m[256];
113 x1p = grid[
IDIR].
x[
i] + dx1; x1m = grid[
IDIR].
x[
i] - dx1;
114 x2p = grid[
JDIR].
x[
j] + dx2; x2m = grid[
JDIR].
x[
j] - dx2;
115 x3p = grid[
KDIR].
x[
k] + dx3; x3m = grid[
KDIR].
x[
k] - dx3;
125 #if GEOMETRY == CARTESIAN
129 Init (us_p, x1f, x2p, x3);
130 Init (us_m, x1f, x2m, x3);
132 bx = (us_p[
AX3] - us_m[
AX3])/dx2;
135 Init (us_p, x1f, x2, x3p);
136 Init (us_m, x1f, x2, x3m);
138 bx -= (us_p[
AX2] - us_m[
AX2])/dx3;
143 Init (us_p, x1p, x2f, x3);
144 Init (us_m, x1m, x2f, x3);
146 by = -(us_p[
AX3] - us_m[
AX3])/dx1;
149 Init (us_p, x1, x2f, x3p);
150 Init (us_m, x1, x2f, x3m);
152 by += (us_p[
AX1] - us_m[
AX1])/dx3;
157 Init (us_p, x1p, x2, x3f);
158 Init (us_m, x1m, x2, x3f);
160 bz = (us_p[
AX2] - us_m[
AX2])/dx1;
162 Init (us_p, x1, x2p, x3f);
163 Init (us_m, x1, x2m, x3f);
165 bz -= (us_p[
AX1] - us_m[
AX1])/dx2;
167 b[0] = bx; b[1] = by; b[2] = bz;
169 #elif GEOMETRY == CYLINDRICAL
173 Init (us_p, x1f, x2p, x3);
174 Init (us_m, x1f, x2m, x3);
176 br = - (us_p[
AX3] - us_m[
AX3])/dx2;
180 Init (us_p, x1p, x2f, x3);
181 Init (us_m, x1m, x2f, x3);
183 bz = (x1p*us_p[
AX3] - x1m*us_m[
AX3])/(x1*dx1);
190 Init (us_p, x1, x2p, x3);
191 Init (us_m, x1, x2m, x3);
193 bphi = (us_p[
AX1] - us_m[
AX1])/dx2;
195 Init (us_p, x1p, x2, x3);
196 Init (us_m, x1m, x2, x3);
198 bphi -= (us_p[
AX2] - us_m[
AX2])/dx1;
201 b[0] = br; b[1] = bz; b[2] = bphi;
203 #elif GEOMETRY == POLAR
207 Init (us_p, x1f, x2p, x3);
208 Init (us_m, x1f, x2m, x3);
210 br = (us_p[
AX3] - us_m[
AX3])/(x1f*dx2);
212 Init (us_p, x1f, x2, x3p);
213 Init (us_m, x1f, x2, x3m);
215 br -= (us_p[
AX2] - us_m[
AX2])/dx3;
219 Init (us_p, x1p, x2f, x3);
220 Init (us_m, x1m, x2f, x3);
222 bphi = -(us_p[
AX3] - us_m[
AX3])/dx1;
225 Init (us_p, x1, x2f, x3p);
226 Init (us_m, x1, x2f, x3m);
228 bphi += (us_p[
AX1] - us_m[
AX1])/dx3;
233 Init (us_p, x1p, x2, x3f);
234 Init (us_m, x1m, x2, x3f);
236 bz = (x1p*us_p[
AX2] - x1m*us_m[
AX2])/(x1*dx1);
238 Init (us_p, x1, x2p, x3f);
239 Init (us_m, x1, x2m, x3f);
241 bz -= (us_p[
AX1] - us_m[
AX1])/(x1*dx2);
243 b[0] = br; b[1] = bphi; b[2] = bz;
245 #elif GEOMETRY == SPHERICAL
249 Init (us_p, x1f, x2p, x3);
250 Init (us_m, x1f, x2m, x3);
252 br = (sin(x2p)*us_p[
AX3] - sin(x2m)*us_m[
AX3])/(x1f*(cos(x2m) - cos(x2p)));
255 Init (us_p, x1f, x2, x3p);
256 Init (us_m, x1f, x2, x3m);
258 br -= 1.0/(x1f*sin(x2)*dx3)*(us_p[
AX2] - us_m[
AX2]);
263 Init (us_p, x1p, x2f, x3);
264 Init (us_m, x1m, x2f, x3);
266 bth = - (x1p*us_p[
AX3] - x1m*us_m[
AX3])/(x1*dx1);
269 Init (us_p, x1, x2f, x3p);
270 Init (us_m, x1, x2f, x3m);
272 bth += (us_p[
AX1] - us_m[
AX1])/(x1*sin(x2f)*dx3);
279 Init (us_p, x1p, x2, x3f);
280 Init (us_m, x1m, x2, x3f);
282 bphi = (x1p*us_p[
AX2] - x1m*us_m[
AX2])/(x1*dx1);
284 Init (us_p, x1, x2p, x3f);
285 Init (us_m, x1, x2m, x3f);
287 bphi -= (us_p[
AX1] - us_m[
AX1])/(x1*dx2);
289 b[0] = br; b[1] = bth; b[2] = bphi;
void Init(double *, double, double, double)
void VectorPotentialDiff(double *, int, int, int, Grid *)