PLUTO
Functions
ct_fill_mag_field.c File Reference

Fill normal staggered component in the ghost zones. More...

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

Go to the source code of this file.

Functions

void FillMagneticField (const Data *d, int side, Grid *grid)
 Assign the normal component of the staggered magnetic field in the ghost zone-faces. More...
 

Detailed Description

Fill normal staggered component in the ghost zones.

Definition in file ct_fill_mag_field.c.

Function Documentation

void FillMagneticField ( const Data *  d,
int  side,
Grid *  grid 
)

Assign the normal component of the staggered magnetic field in the ghost zone-faces.

Using the div.B = 0 condition in staggered MHD, this function computes the staggered component of magnetic field lying on the zone-face parallel to the boundary specified by "side". This is preformed by solving the div.B = 0 condition for one variable only which in 2-D requires the knowledge of the other 3 components while in 3-D required the knowledge of the other 5 staggered components.
Note that this operation is performed in the outermost ghost zones only since the face value at IBEG-1 or IEND is supposed to be part of the solution and is not changed during this function. Therefore, only nghost-1 faces are assigned:

*  +-----+-----+-----+-----+-----+--
*  |     |     |     |     |     |
*  |     X     X     |     |     |
*  |     |     |     |     |     |
*  +-----+-----+-----+-----+-----+--
*                    |
*  <-----------------> BEG
*   Physical boundary     
*
*   X = components assigned in this function.
*
Parameters
[in,out]dpointer to PLUTO Data structure
[in]sidethe side
[in]gridpointer to PLUTO Grid structure
Todo:
replace the loops with more compact macro, such as X1_BEG_LOOP()...
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Aug 16, 2012

Definition at line 8 of file ct_fill_mag_field.c.

50 {
51  int ibeg, jbeg, kbeg;
52  int iend, jend, kend;
53  int i,j,k;
54  int di,dj,dk;
55 
56  real r, dmu;
57  real bxp, byp, bzp;
58  real bxm, bym, bzm;
59  real dBx = 0.0, dBy = 0.0, dBz = 0.0;
60  real dVx, dVy, dVz;
61  real dx, dy, dz;
62  real rp, sp, Axp, Ayp, Azp;
63  real rm, sm, Axm, Aym, Azm;
64 
65  double ***bx, ***by, ***bz; /* -- staggered magnetic field -- */
66  double ***Bx, ***By, ***Bz; /* -- cell centered mag. field -- */
67 
68  D_EXPAND(bx = d->Vs[BX1s]; Bx = d->Vc[BX1]; ,
69  by = d->Vs[BX2s]; By = d->Vc[BX2]; ,
70  bz = d->Vs[BX3s]; Bz = d->Vc[BX3];)
71 
72  ibeg = 0; iend = NX1_TOT-1; di = 1;
73  jbeg = 0; jend = NX2_TOT-1; dj = 1;
74  #if DIMENSIONS == 3
75  kbeg = 0; kend = NX3_TOT-1; dk = 1;
76  #else
77  kbeg = kend = 0, dk = 1;
78  #endif
79 
80  if (side == X1_BEG) {ibeg = IBEG-1; iend = 0; di = -1;}
81  if (side == X1_END) ibeg = IEND+1;
82 
83  if (side == X2_BEG) {jbeg = JBEG-1; jend = 0; dj = -1;}
84  if (side == X2_END) jbeg = JEND+1;
85 
86  if (side == X3_BEG) {kbeg = KBEG-1; kend = 0; dk = -1;}
87  if (side == X3_END) kbeg = KEND+1;
88 
89  for (k = kbeg; dk*k <= dk*kend; k += dk){
90  for (j = jbeg; dj*j <= dj*jend; j += dj){
91  for (i = ibeg; di*i <= di*iend; i += di){
92 
93  r = grid[IDIR].x[i];
94  rp = fabs(grid[IDIR].xr[i]);
95  rm = fabs(grid[IDIR].xl[i]);
96 
97  dx = grid[IDIR].dx[i];
98  dy = grid[JDIR].dx[j];
99  dz = grid[KDIR].dx[k];
100 
101  #if GEOMETRY == SPHERICAL
102  dmu = grid[JDIR].dV[j];
103  sp = grid[JDIR].A[j];
104  sm = grid[JDIR].A[j - 1];
105  #endif
106 
107  #if DIMENSIONS == 2
108  dz = 1.0;
109  #endif
110 
111  D_EXPAND(bxp = bx[k][j][i]; bxm = bx[k][j][i - 1]; ,
112  byp = by[k][j][i]; bym = by[k][j - 1][i]; ,
113  bzp = bz[k][j][i]; bzm = bz[k - 1][j][i];)
114 
115  /* -------------------------------------------------------
116  Divergence is written as
117 
118  dVx*(Axp*bxp - Axm*dxm) +
119  dVy*(Ayp*byp - Aym*dym) +
120  dVz*(Azp*bzp - Azm*dzm) = 0
121 
122  so that the k-th component can be
123  recovered as
124 
125  bkp = bkm*Akm/Akp +
126  sum_(j != k) (Ajp*bjp - Ajm*bjm)*dVj/(dVk*Akp)
127  ------------------------------------------------------- */
128 
129  #if GEOMETRY == CARTESIAN
130 
131  dVx = dy*dz; Axp = Axm = 1.0;
132  dVy = dx*dz; Ayp = Aym = 1.0;
133  dVz = dx*dy; Azp = Azm = 1.0;
134 
135  #elif GEOMETRY == CYLINDRICAL
136 
137  dVx = dy*dz; Axp = rp ; Axm = rm;
138  dVy = r*dx*dz; Ayp = 1.0; Aym = 1.0;
139  dVz = dx*dy; Azp = 1.0; Azm = 1.0;
140 
141  #elif GEOMETRY == POLAR
142 
143  dVx = dy*dz; Axp = rp ; Axm = rm;
144  dVy = dx*dz; Ayp = 1.0; Aym = 1.0;
145  dVz = r*dx*dy; Azp = 1.0; Azm = 1.0;
146 
147  #elif GEOMETRY == SPHERICAL
148 
149  dVx = dmu*dz; Axp = rp*rp; Axm = rm*rm;
150  dVy = r*dx*dz; Ayp = sp ; Aym = sm;
151  dVz = r*dx*dy; Azp = 1.0 ; Azm = 1.0;
152 
153  #endif
154 
155  D_EXPAND(dBx = dVx*(Axp*bxp - Axm*bxm); ,
156  dBy = dVy*(Ayp*byp - Aym*bym); ,
157  dBz = dVz*(Azp*bzp - Azm*bzm); )
158 
159 /* -------------------------------------------------
160  Assign a single face magnetic field
161  ------------------------------------------------- */
162 
163  if (side == X1_BEG){
164 
165  bxm = (bxp*Axp + (dBy + dBz)/dVx)/Axm;
166  bx[k][j][i - 1] = bxm;
167 
168  }else if (side == X1_END){
169 
170  bxp = (bxm*Axm - (dBy + dBz)/dVx)/Axp;
171  bx[k][j][i] = bxp;
172 
173  }else if (side == X2_BEG){
174 
175  bym = (byp*Ayp + (dBx + dBz)/dVy)/Aym;
176  by[k][j - 1][i] = bym;
177 
178  }else if (side == X2_END){
179 
180  byp = (bym*Aym - (dBx + dBz)/dVy)/Ayp;
181  by[k][j][i] = byp;
182 
183  #if DIMENSIONS == 3
184  }else if (side == X3_BEG){
185 
186  bzm = (bzp*Azp + (dBx + dBy)/dVz)/Azm;
187  bz[k - 1][j][i] = bzm;
188 
189  }else if (side == X3_END){
190 
191  bzp = (bzm*Azm - (dBx + dBy)/dVz)/Azp;
192  bz[k][j][i] = bzp;
193 
194  #endif
195  }
196 
197  /* -------------------------------------------------------
198  Now redefine the cell-centered magnetic field
199  ------------------------------------------------------- */
200 /*
201  #if GEOMETRY == CARTESIAN
202  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
203  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
204  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
205  #elif GEOMETRY == CYLINDRICAL
206  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
207  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
208  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
209  #elif GEOMETRY == POLAR
210  D_EXPAND(Bx[k][j][i] = (Axp*bx[k][j][i] + Axm*bx[k][j][i - 1])/(Axp + Axm); ,
211  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
212  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
213  #elif GEOMETRY == SPHERICAL
214  D_EXPAND(Bx[k][j][i] = (Axp*bx[k][j][i] + Axm*bx[k][j][i - 1])/(Axp + Axm); ,
215  By[k][j][i] = (Ayp*by[k][j][i] + Aym*by[k][j - 1][i])/(Ayp + Aym); ,
216  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
217  #endif
218 */
219  }}}
220 
221 }
X3_BEG
#define X3_BEG
Boundary region at X3 beg.
Definition: pluto.h:150
X1_BEG
#define X1_BEG
Boundary region at X1 beg.
Definition: pluto.h:146
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
real
double real
Definition: pluto.h:488
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
BX3s
#define BX3s
Definition: ct.h:29
X1_END
#define X1_END
Boundary region at X1 end.
Definition: pluto.h:147
GRID::dx
double * dx
Definition: structs.h:83
KDIR
#define KDIR
Definition: pluto.h:195
IDIR
#define IDIR
Definition: pluto.h:193
X2_END
#define X2_END
Boundary region at X2 end.
Definition: pluto.h:149
Bx
static double Bx
Definition: hlld.c:62
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
BX1s
#define BX1s
Definition: ct.h:27
GRID::x
double * x
Definition: structs.h:80
GEOMETRY
#define GEOMETRY
Definition: definitions_01.h:4
CARTESIAN
#define CARTESIAN
Definition: pluto.h:33
BX3
#define BX3
Definition: mod_defs.h:27
X3_END
#define X3_END
Boundary region at X3 end.
Definition: pluto.h:151
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
BX1
#define BX1
Definition: mod_defs.h:25
KBEG
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
BX2s
#define BX2s
Definition: ct.h:28
X2_BEG
#define X2_BEG
Boundary region at X2 beg.
Definition: pluto.h:148
BX2
#define BX2
Definition: mod_defs.h:26
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
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: