PLUTO
Macros | Functions
ct_emf_average.c File Reference

Collects different EMF averaging schemes. More...

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

Go to the source code of this file.

Macros

#define EX(k, j, i)   (vz[k][j][i]*By[k][j][i] - vy[k][j][i]*Bz[k][j][i])
 
#define EY(k, j, i)   (vx[k][j][i]*Bz[k][j][i] - vz[k][j][i]*Bx[k][j][i])
 
#define EZ(k, j, i)   (vy[k][j][i]*Bx[k][j][i] - vx[k][j][i]*By[k][j][i])
 
#define dEx_dyp(k, j, i)   (emf->exj[k][j][i] - EX(k,j,i))
 
#define dEx_dzp(k, j, i)   (emf->exk[k][j][i] - EX(k,j,i))
 
#define dEy_dxp(k, j, i)   (emf->eyi[k][j][i] - EY(k,j,i))
 
#define dEy_dzp(k, j, i)   (emf->eyk[k][j][i] - EY(k,j,i))
 
#define dEz_dxp(k, j, i)   (emf->ezi[k][j][i] - EZ(k,j,i))
 
#define dEz_dyp(k, j, i)   (emf->ezj[k][j][i] - EZ(k,j,i))
 
#define dEx_dym(k, j, i)   (EX(k,j,i) - emf->exj[k][j-1][i])
 
#define dEx_dzm(k, j, i)   (EX(k,j,i) - emf->exk[k-1][j][i])
 
#define dEy_dxm(k, j, i)   (EY(k,j,i) - emf->eyi[k][j][i-1])
 
#define dEy_dzm(k, j, i)   (EY(k,j,i) - emf->eyk[k-1][j][i])
 
#define dEz_dxm(k, j, i)   (EZ(k,j,i) - emf->ezi[k][j][i-1])
 
#define dEz_dym(k, j, i)   (EZ(k,j,i) - emf->ezj[k][j-1][i])
 
#define dPP(a, x, y, i, j, k)
 
#define dPM(a, x, y, i, j, k)
 
#define dMM(a, x, y, i, j, k)
 
#define dMP(a, x, y, i, j, k)
 
#define dP(a, x, i, j, k)   (a[k][j][i] + 0.5*(d##a##_##d##x[k][j][i]))
 
#define dM(a, x, i, j, k)   (a[k][j][i] - 0.5*(d##a##_##d##x[k][j][i]))
 

Functions

void CT_EMF_ArithmeticAverage (const EMF *Z1, const double w)
 Compute arithmetic average of EMF at cell edges. More...
 
void CT_EMF_IntegrateToCorner (const Data *d, const EMF *emf, Grid *grid)
 
void CT_EMF_HLL_Solver (const Data *d, const EMF *emf, Grid *grid)
 Solve 2D Riemann problem for induction equation. More...
 
void CT_EMF_CMUSCL_Average (const Data *d, const EMF *emf, Grid *grid)
 – More...
 

Detailed Description

Collects different EMF averaging schemes.

Definition in file ct_emf_average.c.

Macro Definition Documentation

#define dEx_dym (   k,
  j,
  i 
)    (EX(k,j,i) - emf->exj[k][j-1][i])
#define dEx_dyp (   k,
  j,
  i 
)    (emf->exj[k][j][i] - EX(k,j,i))
#define dEx_dzm (   k,
  j,
  i 
)    (EX(k,j,i) - emf->exk[k-1][j][i])
#define dEx_dzp (   k,
  j,
  i 
)    (emf->exk[k][j][i] - EX(k,j,i))
#define dEy_dxm (   k,
  j,
  i 
)    (EY(k,j,i) - emf->eyi[k][j][i-1])
#define dEy_dxp (   k,
  j,
  i 
)    (emf->eyi[k][j][i] - EY(k,j,i))
#define dEy_dzm (   k,
  j,
  i 
)    (EY(k,j,i) - emf->eyk[k-1][j][i])
#define dEy_dzp (   k,
  j,
  i 
)    (emf->eyk[k][j][i] - EY(k,j,i))
#define dEz_dxm (   k,
  j,
  i 
)    (EZ(k,j,i) - emf->ezi[k][j][i-1])
#define dEz_dxp (   k,
  j,
  i 
)    (emf->ezi[k][j][i] - EZ(k,j,i))
#define dEz_dym (   k,
  j,
  i 
)    (EZ(k,j,i) - emf->ezj[k][j-1][i])
#define dEz_dyp (   k,
  j,
  i 
)    (emf->ezj[k][j][i] - EZ(k,j,i))
#define dM (   a,
  x,
  i,
  j,
  k 
)    (a[k][j][i] - 0.5*(d##a##_##d##x[k][j][i]))
#define dMM (   a,
  x,
  y,
  i,
  j,
  k 
)
Value:
(a[k][j][i] - 0.5*(d##a##_##d##x[k][j][i] + \
d##a##_##d##y[k][j][i]))
a
static double a
Definition: init.c:135
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
#define dMP (   a,
  x,
  y,
  i,
  j,
  k 
)
Value:
(a[k][j][i] - 0.5*(d##a##_##d##x[k][j][i] - \
d##a##_##d##y[k][j][i]))
a
static double a
Definition: init.c:135
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
#define dP (   a,
  x,
  i,
  j,
  k 
)    (a[k][j][i] + 0.5*(d##a##_##d##x[k][j][i]))
#define dPM (   a,
  x,
  y,
  i,
  j,
  k 
)
Value:
(a[k][j][i] + 0.5*(d##a##_##d##x[k][j][i] - \
d##a##_##d##y[k][j][i]))
a
static double a
Definition: init.c:135
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
#define dPP (   a,
  x,
  y,
  i,
  j,
  k 
)
Value:
(a[k][j][i] + 0.5*(d##a##_##d##x[k][j][i] + \
d##a##_##d##y[k][j][i]))
a
static double a
Definition: init.c:135
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
#define EX (   k,
  j,
  i 
)    (vz[k][j][i]*By[k][j][i] - vy[k][j][i]*Bz[k][j][i])

Definition at line 8 of file ct_emf_average.c.

#define EY (   k,
  j,
  i 
)    (vx[k][j][i]*Bz[k][j][i] - vz[k][j][i]*Bx[k][j][i])

Definition at line 9 of file ct_emf_average.c.

#define EZ (   k,
  j,
  i 
)    (vy[k][j][i]*Bx[k][j][i] - vx[k][j][i]*By[k][j][i])

Definition at line 10 of file ct_emf_average.c.

Function Documentation

void CT_EMF_ArithmeticAverage ( const EMF *  Z1,
const double  w 
)

Compute arithmetic average of EMF at cell edges.

Combine the four electric field values computed at zone faces as upwind Godunov fluxes into an edge-centered value.
The face-centered EMF should have been stored by previous calls to CT_StoreEMF() during the one-dimensional sweeps.
This function employs a simple arithmetic averaging of the face-centered electric field.

References:

  • "A Staggered Mesh Algorithm Using High Order Godunov Fluxes to Ensure Solenoidal Magnetic Fields in Magnetohydrodynamic Simulations"
    Balsara & Spicer, JCP (1999) 149, 270.
Parameters
[in]Z1pointer to EMF structure
[in]wweighting factor
Returns
This function has no return value.
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Aug 16, 31, 2012

Definition at line 13 of file ct_emf_average.c.

40 {
41  int i, j, k;
42 
43  for (k = Z1->kbeg; k <= Z1->kend; k++){
44  for (j = Z1->jbeg; j <= Z1->jend; j++){
45  for (i = Z1->ibeg; i <= Z1->iend; i++){
46  #if DIMENSIONS == 3
47  Z1->ex[k][j][i] = w*( Z1->exk[k][j][i] + Z1->exk[k][j + 1][i]
48  + Z1->exj[k][j][i] + Z1->exj[k + 1][j][i]);
49  Z1->ey[k][j][i] = w*( Z1->eyi[k][j][i] + Z1->eyi[k + 1][j][i]
50  + Z1->eyk[k][j][i] + Z1->eyk[k][j][i + 1]);
51  #endif
52  Z1->ez[k][j][i] = w*( Z1->ezi[k][j][i] + Z1->ezi[k][j + 1][i]
53  + Z1->ezj[k][j][i] + Z1->ezj[k][j][i + 1]);
54  }}}
55 }
ElectroMotiveForce::ex
double *** ex
Definition: ct.h:103
ElectroMotiveForce::ezi
double *** ezi
Ez available at x-faces (i+1/2);.
Definition: ct.h:79
ElectroMotiveForce::ez
double *** ez
Definition: ct.h:105
ElectroMotiveForce::exk
double *** exk
Ex available at z-faces (k+1/2);.
Definition: ct.h:76
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
ElectroMotiveForce::eyi
double *** eyi
Ey available at x-faces (i+1/2);.
Definition: ct.h:77
ElectroMotiveForce::eyk
double *** eyk
Ey available at z-faces (k+1/2);.
Definition: ct.h:78
i
int i
Definition: analysis.c:2
ElectroMotiveForce::ezj
double *** ezj
Ez available at y-faces (j+1/2);.
Definition: ct.h:80
ElectroMotiveForce::ey
double *** ey
Definition: ct.h:104
ElectroMotiveForce::exj
double *** exj
Ex available at y-faces (j+1/2);.
Definition: ct.h:75

Here is the caller graph for this function:

void CT_EMF_CMUSCL_Average ( const Data *  d,
const EMF *  emf,
Grid *  grid 
)

–

Used in the predictor scheme of CTU scheme in conjunction with UCT_HLL average. No riemann solver actually used, but only a simple pointwise average.

References:

  • "A High order Godunov scheme with constrained transport and adaptie mesh refinement for astrophysical magnetohydrodynamics"
    Fromang, Hennebelle and Teyssier, A&A (2006) 457, 371
Returns
This function has no return value.
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 391 of file ct_emf_average.c.

410 {
411  int i,j,k;
412  int ip, jp, kp;
413  double bx2, by2, bz2;
414  double vx4, vy4, vz4;
415  double ***vx, ***vy, ***vz;
416  double ***bx, ***by, ***bz;
417 
418  D_EXPAND(vx = d->Vc[VX1]; bx = d->Vs[BX1s]; ,
419  vy = d->Vc[VX2]; by = d->Vs[BX2s]; ,
420  vz = d->Vc[VX3]; bz = d->Vs[BX3s];)
421 
422  for (k = emf->kbeg; k <= emf->kend; k++){ kp = k + 1;
423  for (j = emf->jbeg; j <= emf->jend; j++){ jp = j + 1;
424  for (i = emf->ibeg; i <= emf->iend; i++){ ip = i + 1;
425 
426  /* -------------------------------------------
427  EMF: Z component
428  ------------------------------------------- */
429 
430  vx4 = 0.25*(vx[k][j][i] + vx[k][j][ip] + vx[k][jp][i] + vx[k][jp][ip]);
431  vy4 = 0.25*(vy[k][j][i] + vy[k][j][ip] + vy[k][jp][i] + vy[k][jp][ip]);
432 
433  bx2 = 0.5*(bx[k][j][i] + bx[k][jp][i]);
434  by2 = 0.5*(by[k][j][i] + by[k][j][ip]);
435 
436  emf->ez[k][j][i] = vy4*bx2 - vx4*by2;
437 
438  /* -------------------------------------------
439  EMF: X component
440  ------------------------------------------- */
441 
442  #if DIMENSIONS == 3
443  vy4 = 0.25*(vy[k][j][i] + vy[k][jp][i] + vy[kp][j][i] + vy[kp][jp][i]);
444  vz4 = 0.25*(vz[k][j][i] + vz[k][jp][i] + vz[kp][j][i] + vz[kp][jp][i]);
445 
446  by2 = 0.5*(by[k][j][i] + by[kp][j][i]);
447  bz2 = 0.5*(bz[k][j][i] + bz[k][jp][i]);
448 
449  emf->ex[k][j][i] = vz4*by2 - vy4*bz2;
450 
451  /* -------------------------------------------
452  EMF: Y component
453  ------------------------------------------- */
454 
455  vz4 = 0.25*(vz[k][j][i] + vz[kp][j][i] + vz[k][j][ip] + vz[kp][j][ip]);
456  vx4 = 0.25*(vx[k][j][i] + vx[kp][j][i] + vx[k][j][ip] + vx[kp][j][ip]);
457 
458  bz2 = 0.5*(bz[k][j][i] + bz[k][j][ip]);
459  bx2 = 0.5*(bx[k][j][i] + bx[kp][j][i]);
460 
461  emf->ey[k][j][i] = vx4*bz2 - vz4*bx2;
462 
463  #endif
464  }}}
465 
466 }
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
VX2
#define VX2
Definition: mod_defs.h:29
ElectroMotiveForce::ex
double *** ex
Definition: ct.h:103
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
BX3s
#define BX3s
Definition: ct.h:29
VX1
#define VX1
Definition: mod_defs.h:28
ElectroMotiveForce::ez
double *** ez
Definition: ct.h:105
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
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
VX3
#define VX3
Definition: mod_defs.h:30
BX2s
#define BX2s
Definition: ct.h:28
ElectroMotiveForce::ey
double *** ey
Definition: ct.h:104

Here is the call graph for this function:

Here is the caller graph for this function:

void CT_EMF_HLL_Solver ( const Data *  d,
const EMF *  emf,
Grid *  grid 
)

Solve 2D Riemann problem for induction equation.

Solve 2-D Riemann problem using the 2D HLL Riemann flux formula of Londrillo & Del Zanna (2004), JCP, eq. 56. Here N, W, E, S refer to the following configuration: 

*                    |
*                    N
*                 NW | NE
*                    | 
*               --W--+--E--
*                    | 
*                 SW | SE
*                    S
*                    |
*
Parameters
[in]dpointer to PLUTO Data structure
[in]gridpointer to Grid structure;
Returns
This function has no return value.
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 203 of file ct_emf_average.c.

250 {
251  int i,j,k;
252  int ip, jp, kp;
253  double a_xp, a_yp, a_zp;
254  double a_xm, a_ym, a_zm;
255 
256  double eSE, eSW, eNW, eNE;
257  double bS, bN, bW, bE;
258 
259  double *x, *y, *z;
260  double *xr, *yr, *zr;
261  double B0[3];
262 
263  double ***bx, ***by, ***bz;
264  double ***vx, ***vy, ***vz;
265  double ***dvx_dx, ***dvx_dy, ***dvx_dz;
266  double ***dvy_dx, ***dvy_dy, ***dvy_dz;
267  double ***dvz_dx, ***dvz_dy, ***dvz_dz;
268 
269  double ***dbx_dy, ***dbx_dz;
270  double ***dby_dx, ***dby_dz;
271  double ***dbz_dx, ***dbz_dy;
272 
273  D_EXPAND(vx = d->Vc[VX1]; bx = d->Vs[BX1s]; ,
274  vy = d->Vc[VX2]; by = d->Vs[BX2s]; ,
275  vz = d->Vc[VX3]; bz = d->Vs[BX3s];)
276 
277  dvx_dx = emf->dvx_dx; dvx_dy = emf->dvx_dy; dvx_dz = emf->dvx_dz;
278  dvy_dx = emf->dvy_dx; dvy_dy = emf->dvy_dy; dvy_dz = emf->dvy_dz;
279  dvz_dx = emf->dvz_dx; dvz_dy = emf->dvz_dy; dvz_dz = emf->dvz_dz;
280 
281  dbx_dy = emf->dbx_dy; dbx_dz = emf->dbx_dz;
282  dby_dx = emf->dby_dx; dby_dz = emf->dby_dz;
283  dbz_dx = emf->dbz_dx; dbz_dy = emf->dbz_dy;
284 
285  x = grid[IDIR].x; y = grid[JDIR].x; z = grid[KDIR].x;
286  xr = grid[IDIR].xr; yr = grid[JDIR].xr; zr = grid[KDIR].xr;
287 
288  for (k = emf->kbeg; k <= emf->kend; k++){ kp = k + 1;
289  for (j = emf->jbeg; j <= emf->jend; j++){ jp = j + 1;
290  for (i = emf->ibeg; i <= emf->iend; i++){ ip = i + 1;
291 
292  /* -------------------------------------------
293  EMF: Z component at (i+1/2, j+1/2, k)
294  ------------------------------------------- */
295 
296  a_xp = MAX(emf->SxR[k][j][i], emf->SxR[k][jp][i]);
297  a_xm = MAX(emf->SxL[k][j][i], emf->SxL[k][jp][i]);
298  a_yp = MAX(emf->SyR[k][j][i], emf->SyR[k][j][ip]);
299  a_ym = MAX(emf->SyL[k][j][i], emf->SyL[k][j][ip]);
300 
301  bS = dP(bx, y, i , j , k); bW = dP(by, x, i , j , k);
302  bN = dM(bx, y, i , jp, k); bE = dM(by, x, ip, j , k);
303 
304  #if BACKGROUND_FIELD == YES
305  BackgroundField (xr[i], yr[j], z[k], B0);
306  bS += B0[0]; bW += B0[1];
307  bN += B0[0]; bE += B0[1];
308  #endif
309 
310  eSW = dPP(vy,x,y,i ,j ,k)*bS - dPP(vx,x,y,i ,j ,k)*bW;
311  eSE = dMP(vy,x,y,ip,j ,k)*bS - dMP(vx,x,y,ip,j ,k)*bE;
312  eNE = dMM(vy,x,y,ip,jp,k)*bN - dMM(vx,x,y,ip,jp,k)*bE;
313  eNW = dPM(vy,x,y,i ,jp,k)*bN - dPM(vx,x,y,i ,jp,k)*bW;
314 
315  emf->ez[k][j][i] = a_xp*a_yp*eSW + a_xm*a_yp*eSE
316  + a_xm*a_ym*eNE + a_xp*a_ym*eNW;
317  emf->ez[k][j][i] /= (a_xp + a_xm)*(a_yp + a_ym);
318  emf->ez[k][j][i] -= a_yp*a_ym*(bN - bS)/(a_yp + a_ym);
319  emf->ez[k][j][i] += a_xp*a_xm*(bE - bW)/(a_xp + a_xm);
320 
321  /* -------------------------------------------
322  EMF: X component at (i, j+1/2, k+1/2)
323  ------------------------------------------- */
324 
325  #if DIMENSIONS == 3
326  a_xp = MAX(emf->SyR[k][j][i], emf->SyR[kp][j][i]);
327  a_xm = MAX(emf->SyL[k][j][i], emf->SyL[kp][j][i]);
328  a_yp = MAX(emf->SzR[k][j][i], emf->SzR[k][jp][i]);
329  a_ym = MAX(emf->SzL[k][j][i], emf->SzL[k][jp][i]);
330 
331  bS = dP(by, z, i, j, k ); bW = dP(bz, y, i, j , k);
332  bN = dM(by, z, i, j, kp); bE = dM(bz, y, i, jp, k);
333 
334  #if BACKGROUND_FIELD == YES
335  BackgroundField (x[i], yr[j], zr[k], B0);
336  bS += B0[1]; bW += B0[2];
337  bN += B0[1]; bE += B0[2];
338  #endif
339 
340  eSW = dPP(vz,y,z,i,j ,k )*bS - dPP(vy,y,z,i ,j ,k )*bW;
341  eSE = dMP(vz,y,z,i,jp,k )*bS - dMP(vy,y,z,i ,jp,k )*bE;
342  eNE = dMM(vz,y,z,i,jp,kp)*bN - dMM(vy,y,z,i ,jp,kp)*bE;
343  eNW = dPM(vz,y,z,i,j ,kp)*bN - dPM(vy,y,z,i ,j ,kp)*bW;
344 
345  emf->ex[k][j][i] = a_xp*a_yp*eSW + a_xm*a_yp*eSE
346  + a_xm*a_ym*eNE + a_xp*a_ym*eNW;
347  emf->ex[k][j][i] /= (a_xp + a_xm)*(a_yp + a_ym);
348  emf->ex[k][j][i] -= a_yp*a_ym*(bN - bS)/(a_yp + a_ym);
349  emf->ex[k][j][i] += a_xp*a_xm*(bE - bW)/(a_xp + a_xm);
350 
351  /* -------------------------------------------
352  EMF: Y component at (i+1/2, j, k+1/2)
353  ------------------------------------------- */
354 
355  a_xp = MAX(emf->SzR[k][j][i], emf->SzR[k][j][ip]);
356  a_xm = MAX(emf->SzL[k][j][i], emf->SzL[k][j][ip]);
357  a_yp = MAX(emf->SxR[k][j][i], emf->SxR[kp][j][i]);
358  a_ym = MAX(emf->SxL[k][j][i], emf->SxL[kp][j][i]);
359 
360  bS = dP(bz, x, i , j, k); bW = dP(bx, z, i, j, k);
361  bN = dM(bz, x, ip, j, k); bE = dM(bx, z, i, j, kp);
362 
363  #if BACKGROUND_FIELD == YES
364  BackgroundField (xr[i], y[j], zr[k], B0);
365  bS += B0[2]; bW += B0[0];
366  bN += B0[2]; bE += B0[0];
367  #endif
368 
369  eSW = dPP(vx,z,x,i, j,k )*bS - dPP(vz,z,x,i ,j ,k )*bW;
370  eSE = dMP(vx,z,x,i, j,kp)*bS - dMP(vz,z,x,i ,j ,kp)*bE;
371  eNE = dMM(vx,z,x,ip,j,kp)*bN - dMM(vz,z,x,ip,j ,kp)*bE;
372  eNW = dPM(vx,z,x,ip,j,k )*bN - dPM(vz,z,x,ip,j ,k )*bW;
373 
374  emf->ey[k][j][i] = a_xp*a_yp*eSW + a_xm*a_yp*eSE
375  + a_xm*a_ym*eNE + a_xp*a_ym*eNW;
376  emf->ey[k][j][i] /= (a_xp + a_xm)*(a_yp + a_ym);
377  emf->ey[k][j][i] -= a_yp*a_ym*(bN - bS)/(a_yp + a_ym);
378  emf->ey[k][j][i] += a_xp*a_xm*(bE - bW)/(a_xp + a_xm);
379  #endif /* DIMENSIONS == 3 */
380  }}}
381 
382 }
MAX
#define MAX(a, b)
Definition: macros.h:101
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
dP
#define dP(a, x, i, j, k)
VX2
#define VX2
Definition: mod_defs.h:29
ElectroMotiveForce::SxL
double *** SxL
Definition: ct.h:93
ElectroMotiveForce::SzR
double *** SzR
Definition: ct.h:98
dPP
#define dPP(a, x, y, i, j, k)
ElectroMotiveForce::ex
double *** ex
Definition: ct.h:103
BackgroundField
void BackgroundField(double x1, double x2, double x3, double *B0)
Definition: init.c:79
ElectroMotiveForce::SxR
double *** SxR
Definition: ct.h:94
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
BX3s
#define BX3s
Definition: ct.h:29
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
ElectroMotiveForce::ez
double *** ez
Definition: ct.h:105
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
ElectroMotiveForce::SyR
double *** SyR
Definition: ct.h:96
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
ElectroMotiveForce::SzL
double *** SzL
Definition: ct.h:97
dM
#define dM(a, x, i, j, k)
VX3
#define VX3
Definition: mod_defs.h:30
BX2s
#define BX2s
Definition: ct.h:28
dPM
#define dPM(a, x, y, i, j, k)
ElectroMotiveForce::ey
double *** ey
Definition: ct.h:104
JDIR
#define JDIR
Definition: pluto.h:194
dMP
#define dMP(a, x, y, i, j, k)
ElectroMotiveForce::SyL
double *** SyL
Definition: ct.h:95
dMM
#define dMM(a, x, y, i, j, k)

Here is the call graph for this function:

Here is the caller graph for this function:

void CT_EMF_IntegrateToCorner ( const Data *  d,
const EMF *  emf,
Grid *  grid 
)

Add derivatives to the 4-point arithmetic average of magnetic fields. Obtain the electric field at corners.

References:

  • "An unsplit Godunov method for ideal MHD via constrained transport"
    Gardiner & Stone, JCP (2005) 205, 509. See Eq. (41), (45) and (50).
Parameters
[in]dpointer to PLUTO Data structure
[in]emfpointer to EMF structure
[in]gridpointer to Grid structure
Returns
This function has no return value.
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Aug 16, 31, 2012

Definition at line 58 of file ct_emf_average.c.

98 {
99  int i, j, k;
100  int iu, ju, ku;
101  signed char sx, sy, sz;
102  double ***vx, ***vy, ***vz;
103  double ***Bx, ***By, ***Bz;
104 
105  #ifdef CTU
106  if (g_intStage == 1) return; /* -- not needed in predictor step of CTU -- */
107  #endif
108 
109  EXPAND(vx = d->Vc[VX1]; Bx = d->Vc[BX1]; ,
110  vy = d->Vc[VX2]; By = d->Vc[BX2]; ,
111  vz = d->Vc[VX3]; Bz = d->Vc[BX3];)
112 
113  for (k = emf->kbeg; k <= emf->kend + KOFFSET; k++){
114  for (j = emf->jbeg; j <= emf->jend + 1 ; j++){
115  for (i = emf->ibeg; i <= emf->iend + 1 ; i++){
116 
117  D_EXPAND(sx = emf->svx[k][j][i]; ,
118  sy = emf->svy[k][j][i]; ,
119  sz = emf->svz[k][j][i];)
120 
121  D_EXPAND(iu = sx > 0 ? i:i+1; , /* -- upwind index -- */
122  ju = sy > 0 ? j:j+1; ,
123  ku = sz > 0 ? k:k+1;)
124 
125  /* ----------------------------------------
126  Span X - Faces: dEz/dy, dEy/dz
127  ---------------------------------------- */
128 
129  if (sx == 0) {
130  emf->ez[k][j][i] += 0.5*(dEz_dyp(k,j,i) + dEz_dyp(k,j,i+1));
131  emf->ez[k][j-1][i] -= 0.5*(dEz_dym(k,j,i) + dEz_dym(k,j,i+1));
132  #if DIMENSIONS == 3
133  emf->ey[k][j][i] += 0.5*(dEy_dzp(k,j,i) + dEy_dzp(k,j,i+1));
134  emf->ey[k-1][j][i] -= 0.5*(dEy_dzm(k,j,i) + dEy_dzm(k,j,i+1));
135  #endif
136  }else{
137  emf->ez[k][j][i] += dEz_dyp(k,j,iu);
138  emf->ez[k][j-1][i] -= dEz_dym(k,j,iu);
139  #if DIMENSIONS == 3
140  emf->ey[k][j][i] += dEy_dzp(k,j,iu);
141  emf->ey[k-1][j][i] -= dEy_dzm(k,j,iu);
142  #endif
143  }
144 
145  /* ----------------------------------------
146  Span Y - Faces: dEz/dx, dEx/dz
147  ---------------------------------------- */
148 
149  if (sy == 0) {
150  emf->ez[k][j][i] += 0.5*(dEz_dxp(k,j,i) + dEz_dxp(k,j+1,i));
151  emf->ez[k][j][i-1] -= 0.5*(dEz_dxm(k,j,i) + dEz_dxm(k,j+1,i));
152  #if DIMENSIONS == 3
153  emf->ex[k][j][i] += 0.5*(dEx_dzp(k,j,i) + dEx_dzp(k,j+1,i));
154  emf->ex[k-1][j][i] -= 0.5*(dEx_dzm(k,j,i) + dEx_dzm(k,j+1,i));
155  #endif
156  }else{
157  emf->ez[k][j][i] += dEz_dxp(k,ju,i);
158  emf->ez[k][j][i-1] -= dEz_dxm(k,ju,i);
159  #if DIMENSIONS == 3
160  emf->ex[k][j][i] += dEx_dzp(k,ju,i);
161  emf->ex[k-1][j][i] -= dEx_dzm(k,ju,i);
162  #endif
163  }
164 
165  /* ----------------------------------------
166  Span Z - Faces: dEx/dy, dEy/dx
167  ---------------------------------------- */
168 
169  #if DIMENSIONS == 3
170 
171  if (sz == 0) {
172  emf->ex[k][j][i] += 0.5*(dEx_dyp(k,j,i) + dEx_dyp(k+1,j,i));
173  emf->ex[k][j-1][i] -= 0.5*(dEx_dym(k,j,i) + dEx_dym(k+1,j,i));
174  emf->ey[k][j][i] += 0.5*(dEy_dxp(k,j,i) + dEy_dxp(k+1,j,i));
175  emf->ey[k][j][i-1] -= 0.5*(dEy_dxm(k,j,i) + dEy_dxm(k+1,j,i));
176  }else{
177  emf->ex[k][j][i] += dEx_dyp(ku,j,i);
178  emf->ex[k][j-1][i] -= dEx_dym(ku,j,i);
179  emf->ey[k][j][i] += dEy_dxp(ku,j,i);
180  emf->ey[k][j][i-1] -= dEy_dxm(ku,j,i);
181  }
182 
183  #endif
184  }}}
185 }
dEy_dzm
#define dEy_dzm(k, j, i)
VX2
#define VX2
Definition: mod_defs.h:29
ElectroMotiveForce::ex
double *** ex
Definition: ct.h:103
ElectroMotiveForce::svy
signed char *** svy
Definition: ct.h:83
g_intStage
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
Definition: globals.h:98
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
dEy_dxp
#define dEy_dxp(k, j, i)
VX1
#define VX1
Definition: mod_defs.h:28
dEy_dxm
#define dEy_dxm(k, j, i)
dEy_dzp
#define dEy_dzp(k, j, i)
ElectroMotiveForce::ez
double *** ez
Definition: ct.h:105
Bx
static double Bx
Definition: hlld.c:62
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
dEz_dyp
#define dEz_dyp(k, j, i)
ElectroMotiveForce::svz
signed char *** svz
Definition: ct.h:83
BX3
#define BX3
Definition: mod_defs.h:27
ElectroMotiveForce::svx
signed char *** svx
Definition: ct.h:83
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
dEx_dzp
#define dEx_dzp(k, j, i)
dEx_dyp
#define dEx_dyp(k, j, i)
i
int i
Definition: analysis.c:2
dEx_dym
#define dEx_dym(k, j, i)
dEx_dzm
#define dEx_dzm(k, j, i)
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
dEz_dym
#define dEz_dym(k, j, i)
dEz_dxp
#define dEz_dxp(k, j, i)
BX2
#define BX2
Definition: mod_defs.h:26
ElectroMotiveForce::ey
double *** ey
Definition: ct.h:104
dEz_dxm
#define dEz_dxm(k, j, i)

Here is the call graph for this function:

Here is the caller graph for this function: