PLUTO
Functions
ct.c File Reference

Update staggered magnetic field. More...

#include "pluto.h"
Include dependency graph for ct.c:

Go to the source code of this file.

Functions

void CT_Update (const Data *d, Data_Arr Bs, double dt, Grid *grid)
 
void CT_CheckDivB (double ***bf[], Grid *grid)
 

Detailed Description

Update staggered magnetic field.

Update face-centered magnetic field in the constrained transport formulation using a discrete version of Stoke's theorem. The update consists of a single Euler step:

\[ \mathtt{d->Vs} = \mathtt{Bs} + \Delta t R \]

where d->Vs is the main staggered array used by PLUTO, Bs is the magnetic field to be updated and R is the right hand side already computed during the unsplit integrator. d->Vs and Bs may be the same array or may be different.

References

  • "A staggered mesh algorithm using high-order Godunov fluxes to ensure solenoidal magnetic field in MHD simulations"
    Balsara & Spicer, JCP (1999) 149, 270
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
April 1, 2014

Definition in file ct.c.

Function Documentation

void CT_CheckDivB ( double ***  bf[],
Grid *  grid 
)

Check the divergence-free condition of magnetic field in the constrained transport formalism.

Definition at line 201 of file ct.c.

207 {
208  int i,j,k;
209  double divB;
210  double ***bx, ***by, ***bz;
211  double *dx, *dy, *dz;
212  double *Ar, *s, *dmu, *r, *rp, *rm;
213  double dbmax=0.0;
214 
215  D_EXPAND( bx = bf[0]; ,
216  by = bf[1]; ,
217  bz = bf[2]; )
218 
219  rp = grid[IDIR].xr;
220  rm = grid[IDIR].xl;
221 
222  Ar = grid[IDIR].A;
223  s = grid[JDIR].A;
224  dmu = grid[JDIR].dV;
225 
226  r = grid[IDIR].x;
227  dx = grid[IDIR].dx;
228  dy = grid[JDIR].dx;
229  dz = grid[KDIR].dx;
230 
231 
232 /*
233 for (i = IBEG-1; i >= -1; i--){
234 for (j=0; j < NX2_TOT; j++){
235  divB = bx[0][j][i]-bx[0][j][i+NX];
236  if (fabs(divB) > 1.e-12){
237  printf ("!! periodicity not respected (bx), zone = %d, %d\n",i,j);
238  exit(1);
239  }
240 }}
241 for (i = IBEG-1; i >= 0; i--){
242 for (j = -1; j < NX2_TOT; j++){
243  divB = by[0][j][i]-by[0][j][i+NX];
244  if (fabs(divB) > 1.e-12){
245  printf ("!! periodicity not respected (by)\n");
246  exit(1);
247  }
248 }}
249 */
250  for (k = 0; k < NX3_TOT; k++){
251 
252  #if DIMENSIONS < 3
253  dz[k] = 1.0;
254  #endif
255 
256  for (j = 0; j < NX2_TOT; j++){
257  for (i = 0; i < NX1_TOT; i++){
258 
259  #if GEOMETRY == CARTESIAN
260 
261  divB = D_EXPAND( (bx[k][j][i] - bx[k][j][i-1])*dy[j]*dz[k] ,
262  + (by[k][j][i] - by[k][j-1][i])*dx[i]*dz[k] ,
263  + (bz[k][j][i] - bz[k-1][j][i])*dx[i]*dy[j] );
264 
265  #elif GEOMETRY == CYLINDRICAL
266 
267  divB = dy[j]*(bx[k][j][i]*Ar[i] - bx[k][j][i-1]*Ar[i-1])
268  + fabs(r[i])*dx[i]*(by[k][j][i] - by[k][j-1][i]);
269 
270 /*
271  divB = dy[j]*(bx[k][j][i]*fabs(rp[i]) - bx[k][j][i-1]*fabs(rm[i]))
272  + fabs(r[i])*dx[i]*(by[k][j][i] - by[k][j-1][i]);
273 */
274  #elif GEOMETRY == POLAR
275 
276  divB = D_EXPAND(
277  dz[k]*dy[j]*(bx[k][j][i]*Ar[i] - bx[k][j][i-1]*Ar[i-1]) ,
278  + dx[i]*dz[k]*(by[k][j][i] - by[k][j-1][i]) ,
279  + r[i]*dx[i]*dy[j]*(bz[k][j][i] - bz[k-1][j][i])
280  );
281 
282  #elif GEOMETRY == SPHERICAL
283 
284  divB = D_EXPAND(
285  dmu[j] *dz[k]*(bx[k][j][i]*Ar[i] - bx[k][j][i-1]*Ar[i-1]) ,
286  + r[i]*dx[i]*dz[k]*(by[k][j][i]*s[j] - by[k][j-1][i]*s[j-1]) ,
287  + r[i]*dx[i]*dy[j]*(bz[k][j][i] - bz[k-1][j][i])
288  );
289 
290  #endif
291 
292  dbmax = MAX(dbmax,fabs(divB));
293  if (fabs(divB) > 1.e-6) {
294  print ("! CHECK_DIVB: div(B) = %12.6e != 0, rank = %d, ijk = %d %d %d \n",
295  divB, prank, i,j,k);
296  print ("! rank =%d, g_intStage = %d, Beg, End = [%d, %d] [%d, %d] [%d, %d]\n",
297  prank, g_intStage, IBEG, IEND,
298  JBEG, JEND, KBEG, KEND);
299  D_EXPAND( print ("b1: %12.6e %12.6e\n",bx[k][j][i],bx[k][j][i-1]); ,
300  print ("b2: %12.6e %12.6e\n",by[k][j][i],by[k][j-1][i]); ,
301  print ("b3: %12.6e %12.6e\n",bz[k][j][i],bz[k-1][j][i]); )
302 
303  QUIT_PLUTO(1);
304  }
305  }}
306  }
307 
308 }
MAX
#define MAX(a, b)
Definition: macros.h:101
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
g_intStage
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
Definition: globals.h:98
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
GRID::dx
double * dx
Definition: structs.h:83
prank
int prank
Processor rank.
Definition: globals.h:33
KDIR
#define KDIR
Definition: pluto.h:195
GRID::xl
double * xl
Definition: structs.h:82
IDIR
#define IDIR
Definition: pluto.h:193
NX3_TOT
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
j
int j
Definition: analysis.c:2
IEND
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:37
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
divB
double *** divB
Definition: analysis.c:4
GRID::x
double * x
Definition: structs.h:80
s
#define s
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
i
int i
Definition: analysis.c:2
KBEG
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
KEND
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:45
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
JEND
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:41
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the call graph for this function:

Here is the caller graph for this function:

void CT_Update ( const Data *  d,
Data_Arr  Bs,
double  dt,
Grid *  grid 
)

Update staggered magnetic field using discrete version of Stoke's theorem. Only d->Vs is updated, while Bs is the original array:
d->Vs = Bs + dt * R, where R = curl(E) is the electric field.

Parameters
[in,out]dpointer to PLUTO Data structure. d->Vs will be updated.
[in]Bsthe original array (not updated).
[in]dtstep size
[in]gridpointer to an array of Grid structures

Definition at line 29 of file ct.c.

43 {
44  int i, j, k, nv;
45  int ibeg, jbeg, kbeg;
46  int iend, jend, kend;
47  double rhs_x, rhs_y, rhs_z;
48  double *dx, *dy, *dz, *A1, *A2, *dV1, *dV2;
49  double *r, *rp, *th, r_2;
50  double ***ex, ***ey, ***ez;
51  EMF *emf;
52 
53 /* ---- check div.B ---- */
54 
55  #if CHECK_DIVB_CONDITION == YES
56  if (g_intStage == 1) CT_CheckDivB (d->Vs, grid);
57  #endif
58 
59  emf = CT_GetEMF (d, grid);
60 
61  #if UPDATE_VECTOR_POTENTIAL == YES
62  VectorPotentialUpdate (d, emf, NULL, grid);
63  #endif
64 
65  ex = emf->ex;
66  ey = emf->ey;
67  ez = emf->ez;
68 
69 /* ------------------------------------------------------------
70  Correct EMF if the ShearingBox module is enabled.
71  For CTU, this step is carried only during the final stage.
72  ------------------------------------------------------------ */
73 
74  #ifdef SHEARINGBOX
75  #ifdef CTU
76  if (g_intStage == 2) SB_CorrectEMF (emf, Bs, grid);
77  #else
78  SB_CorrectEMF (emf, d->Vs, grid);
79  #endif
80  #endif
81 
82  dx = grid[IDIR].dx;
83  dy = grid[JDIR].dx;
84  dz = grid[KDIR].dx;
85 
86  r = grid[IDIR].x;
87  rp = grid[IDIR].xr;
88  A1 = grid[IDIR].A;
89  dV1 = grid[IDIR].dV;
90 
91  th = grid[JDIR].x;
92  A2 = grid[JDIR].A;
93  dV2 = grid[JDIR].dV;
94 
95 /* --------------------------------------------------------
96  Update Bx1 at (i+1/2, j, k) faces
97  -------------------------------------------------------- */
98 
99  for (k = emf->kbeg + KOFFSET; k <= emf->kend; k++){
100  for (j = emf->jbeg + 1 ; j <= emf->jend; j++){
101  for (i = emf->ibeg ; i <= emf->iend; i++){
102 
103  #if GEOMETRY == CARTESIAN
104 
105  rhs_x = D_EXPAND( 0.0 ,
106  - dt/dy[j]*(ez[k][j][i] - ez[k][j - 1][i]) ,
107  + dt/dz[k]*(ey[k][j][i] - ey[k - 1][j][i]) );
108 
109  #elif GEOMETRY == CYLINDRICAL
110 
111  rhs_x = - dt/dy[j]*(ez[k][j][i] - ez[k][j - 1][i]);
112 
113  #elif GEOMETRY == POLAR
114 
115  rhs_x = D_EXPAND( 0.0 ,
116  - dt/(A1[i]*dy[j])*(ez[k][j][i] - ez[k][j - 1][i]) ,
117  + dt/dz[k] *(ey[k][j][i] - ey[k - 1][j][i]));
118 
119  #elif GEOMETRY == SPHERICAL
120 
121  rhs_x = D_EXPAND(
122  0.0 ,
123  - dt/(rp[i]*dV2[j])*(A2[j]*ez[k][j][i] - A2[j - 1]*ez[k][j - 1][i]) ,
124  + dt*dy[j]/(rp[i]*dV2[j]*dz[k])*(ey[k][j][i] - ey[k - 1][j][i]));
125 
126  #endif
127 
128  d->Vs[BX1s][k][j][i] = Bs[BX1s][k][j][i] + rhs_x;
129 
130  }}}
131 
132 /* --------------------------------------------------------
133  Update Bx2 at (i, j+1/2, k) faces
134  -------------------------------------------------------- */
135 
136  for (k = emf->kbeg + KOFFSET; k <= emf->kend; k++){
137  for (j = emf->jbeg ; j <= emf->jend; j++){
138  for (i = emf->ibeg + 1 ; i <= emf->iend; i++){
139 
140  #if GEOMETRY == CARTESIAN
141 
142  rhs_y = D_EXPAND( dt/dx[i]*(ez[k][j][i] - ez[k][j][i - 1]) ,
143  ,
144  - dt/dz[k]*(ex[k][j][i] - ex[k - 1][j][i]) );
145 
146  #elif GEOMETRY == CYLINDRICAL
147 
148  rhs_y = dt/dV1[i]*(A1[i]*ez[k][j][i] - A1[i - 1]*ez[k][j][i - 1]);
149 
150  #elif GEOMETRY == POLAR
151 
152  rhs_y = D_EXPAND( dt/dx[i]*(ez[k][j][i] - ez[k][j][i - 1]) ,
153  ,
154  - dt/dz[k]*(ex[k][j][i] - ex[k - 1][j][i]));
155 
156  #elif GEOMETRY == SPHERICAL
157 
158  rhs_y = D_EXPAND(
159  + dt/(r[i]*dx[i])*(rp[i]*ez[k][j][i] - rp[i - 1]*ez[k][j][i - 1]) ,
160  , - dt/(r[i]*A2[j]*dz[k])*(ex[k][j][i] - ex[k - 1][j][i]));
161 
162  #endif
163 
164  d->Vs[BX2s][k][j][i] = Bs[BX2s][k][j][i] + rhs_y;
165 
166  }}}
167 
168 /* --------------------------------------------------------
169  Update Bx3 at (i, j, k+1/2) faces
170  -------------------------------------------------------- */
171 
172  #if DIMENSIONS == 3
173  for (k = emf->kbeg ; k <= emf->kend; k++){
174  for (j = emf->jbeg + 1; j <= emf->jend; j++){
175  for (i = emf->ibeg + 1; i <= emf->iend; i++){
176 
177  #if GEOMETRY == CARTESIAN
178 
179  rhs_z = - dt/dx[i]*(ey[k][j][i] - ey[k][j][i - 1])
180  + dt/dy[j]*(ex[k][j][i] - ex[k][j - 1][i]);
181 
182  #elif GEOMETRY == POLAR
183 
184  rhs_z = - dt/dV1[i]*(A1[i]*ey[k][j][i] - A1[i - 1]*ey[k][j][i - 1])
185  + dt/(r[i]*dy[j])*(ex[k][j][i] - ex[k][j - 1][i]);
186 
187  #elif GEOMETRY == SPHERICAL
188 
189  rhs_z = - dt/(r[i]*dx[i])*(rp[i]*ey[k][j][i] - rp[i - 1]*ey[k][j][i - 1])
190  + dt/(r[i]*dy[j])*(ex[k][j][i] - ex[k][j - 1][i]);
191  #endif
192 
193  d->Vs[BX3s][k][j][i] = Bs[BX3s][k][j][i] + rhs_z;
194 
195  }}}
196  #endif
197 
198 }
emf
static EMF emf
Definition: ct_emf.c:27
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
ElectroMotiveForce::ex
double *** ex
Definition: ct.h:103
GRID::xr
double * xr
Definition: structs.h:81
g_intStage
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
Definition: globals.h:98
CT_GetEMF
EMF * CT_GetEMF(const Data *d, Grid *grid)
Definition: ct_emf.c:205
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
BX3s
#define BX3s
Definition: ct.h:29
GRID::dx
double * dx
Definition: structs.h:83
KDIR
#define KDIR
Definition: pluto.h:195
ElectroMotiveForce::ez
double *** ez
Definition: ct.h:105
IDIR
#define IDIR
Definition: pluto.h:193
CT_CheckDivB
void CT_CheckDivB(double ***bf[], Grid *grid)
Definition: ct.c:201
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
GRID::x
double * x
Definition: structs.h:80
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
i
int i
Definition: analysis.c:2
BX2s
#define BX2s
Definition: ct.h:28
VectorPotentialUpdate
void VectorPotentialUpdate(const Data *d, const void *vp, const State_1D *state, const Grid *grid)
ElectroMotiveForce::ey
double *** ey
Definition: ct.h:104
JDIR
#define JDIR
Definition: pluto.h:194
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the call graph for this function:

Here is the caller graph for this function: