glm.h File Reference

Header file for GLM Divergence Cleaning. More...

Go to the source code of this file.

Macros
#define	GLM_MHD
#define	GLM_ALPHA   0.1
	Sets the damping rate of monopoles. More...
#define	GLM_EXTENDED   NO
	The GLM_EXTENDED macro may be turned to YES to enable the extended GLM formalism. More...
#define	COMPUTE_DIVB   NO

Functions
void	GLM_Solve (const State_1D *, double **, double **, int, int, Grid *)
void	GLM_SolveNEW (const State_1D *state, int beg, int end, Grid *grid)
void	GLM_Init (const Data *, const Time_Step *, Grid *)
void	GLM_Source (const Data_Arr, double, Grid *)
void	GLM_ExtendedSource (const State_1D *, double, int, int, Grid *)
void	GLM_ComputeDivB (const State_1D *state, Grid *grid)
double ***	GLM_GetDivB (void)

Variables
double	glm_ch
	The propagation speed of divergence error. More...

Detailed Description

Header file for GLM Divergence Cleaning.

Contains function prototypes and global variable declaration for the GLM formulation to control the divergence-free condition of magnetic field.

Authors

A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
P. Tzeferacos (petro.nosp@m.s.tz.nosp@m.efera.nosp@m.cos@.nosp@m.ph.un.nosp@m.ito..nosp@m.it)

Date

July 24, 2015

References

• "A Second-order unsplit Godunov scheme for cell-centered MHD: The CTU-GLM scheme"
  Mignone & Tzeferacos, JCP (2010) 229, 2117
• "High-order conservative finite difference GLM-MHD scheme for cell-centered MHD"
  Mignone, Tzeferacos & Bodo, JCP (2010) 229, 5896

Definition in file glm.h.

Macro Definition Documentation

#define COMPUTE_DIVB   NO

Definition at line 44 of file glm.h.

#define GLM_ALPHA   0.1

Sets the damping rate of monopoles.

Definition at line 28 of file glm.h.

#define GLM_EXTENDED   NO

The GLM_EXTENDED macro may be turned to YES to enable the extended GLM formalism.

Although it breaks conservation of momentum and energy, it has proven to be more robust in treating low-beta plasma.

Definition at line 32 of file glm.h.

#define GLM_MHD

Definition at line 25 of file glm.h.

Function Documentation

void GLM_ComputeDivB	(	const State_1D *	state,
		Grid *	grid
	)

Definition at line 411 of file glm.c.

{
  int    i,j,k;
  static int old_nstep = -1;
  double   *B, *A, *dV;

  if (divB == NULL) divB = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);

  if (old_nstep != g_stepNumber){
    old_nstep = g_stepNumber;
    DOM_LOOP(k,j,i) divB[k][j][i] = 0.0;
  }

  B  = state->bn;
  A  = grid[g_dir].A;
  dV = grid[g_dir].dV;
  if (g_dir == IDIR){
    k = g_k; j = g_j;
    IDOM_LOOP(i) divB[k][j][i] += (A[i]*B[i] - A[i-1]*B[i-1])/dV[i];
  }else if (g_dir == JDIR){
    k = g_k; i = g_i;
    JDOM_LOOP(j) divB[k][j][i] += (A[j]*B[j] - A[j-1]*B[j-1])/dV[j];
  }else if (g_dir == KDIR){
    j = g_j; i = g_i;
    KDOM_LOOP(k) divB[k][j][i] += (A[k]*B[k] - A[k-1]*B[k-1])/dV[k];
  } 
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GLM_ExtendedSource	(	const State_1D *	state,
		double	dt,
		int	beg,
		int	end,
		Grid *	grid
	)

Add source terms to the right hand side of the conservative equations, momentum and energy equations only. This yields the extended GLM equations given by Eq. (24a)–(24c) in

"Hyperbolic Divergence cleaning for the MHD Equations" Dedner et al. (2002), JcP, 175, 645

Definition at line 168 of file glm.c.

{
  int    i;
  double bx, by, bz;
  double ch2, r, s;
  double *Ar, *Ath, **bgf;
  double *rhs, *v, *Bm, *pm;
  static double *divB, *dpsi, *Bp, *pp;
  Grid   *GG;
/* 
  #if PHYSICS == RMHD
   print1 ("! EGLM not working for the RMHD module\n");
   QUIT_PLUTO(1);
  #endif
*/
  if (divB == NULL){
    divB = ARRAY_1D(NMAX_POINT, double);
    dpsi = ARRAY_1D(NMAX_POINT, double);
    Bp   = ARRAY_1D(NMAX_POINT, double);
    pp   = ARRAY_1D(NMAX_POINT, double);
  }

/* ----------------------- ---------------------
    compute magnetic field normal component 
    interface value by arithmetic averaging
   -------------------------------------------- */

  ch2 = glm_ch*glm_ch;
  for (i = beg - 1; i <= end; i++) {
    Bp[i] = state->flux[i][PSI_GLM]/ch2;
    pp[i] = state->flux[i][BXn];
  }
  Bm  = Bp - 1;
  pm  = pp - 1;
  GG  = grid + g_dir;
  
/* --------------------------------------------
    Compute div.B contribution from the normal 
    direction (1) in different geometries 
   -------------------------------------------- */
  
  #if GEOMETRY == CARTESIAN

   for (i = beg; i <= end; i++) {
     divB[i] = (Bp[i] - Bm[i])/GG->dx[i];
     dpsi[i] = (pp[i] - pm[i])/GG->dx[i];
   }

  #elif GEOMETRY == CYLINDRICAL

   if (g_dir == IDIR){   /* -- r -- */
     Ar  = grid[IDIR].A;
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i]*Ar[i] - Bm[i]*Ar[i - 1])/GG->dV[i];
       dpsi[i] = (pp[i] - pm[i])/GG->dx[i];
     }
   }else if (g_dir == JDIR){  /* -- z -- */
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i] - Bm[i])/GG->dx[i];
       dpsi[i] = (pp[i] - pm[i])/GG->dx[i];
     }
   }

  #elif GEOMETRY == POLAR

   if (g_dir == IDIR){  /* -- r -- */
     Ar  = grid[IDIR].A;
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i]*Ar[i] - Bm[i]*Ar[i - 1])/GG->dV[i];
       dpsi[i] = (pp[i] - pm[i])/GG->dx[i];
     }
   }else if (g_dir == JDIR){  /* -- phi -- */
     r = grid[IDIR].x[g_i];
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i] - Bm[i])/(r*GG->dx[i]);
       dpsi[i] = (pp[i] - pm[i])/(r*GG->dx[i]);
     }
   }else if (g_dir == KDIR){  /* -- z -- */
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i] - Bm[i])/GG->dx[i];
       dpsi[i] = (pp[i] - pm[i])/GG->dx[i];
     }
   }

  #elif GEOMETRY == SPHERICAL

   if (g_dir == IDIR){  /* -- r -- */
     Ar  = grid[IDIR].A;
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i]*Ar[i] - Bm[i]*Ar[i - 1])/GG->dV[i];
       dpsi[i] = (pp[i] - pm[i])/GG->dx[i];
     }
   }else if (g_dir == JDIR){  /* -- theta -- */
     Ath = grid[JDIR].A;
     r   = grid[IDIR].x[g_i];
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i]*Ath[i] - Bm[i]*Ath[i - 1]) /
                 (r*GG->dV[i]);
       dpsi[i] = (pp[i] - pm[i])/(r*GG->dx[i]);

     }
   }else if (g_dir == KDIR){  /* -- phi -- */
     r = grid[IDIR].x[g_i];
     s = sin(grid[JDIR].x[g_j]);
     for (i = beg; i <= end; i++) {
       divB[i] = (Bp[i] - Bm[i])/(r*s*GG->dx[i]);
       dpsi[i] = (pp[i] - pm[i])/(r*s*GG->dx[i]);
     }
   }

  #endif

  #if BACKGROUND_FIELD == YES
   bgf = GetBackgroundField (beg - 1, end, CELL_CENTER, grid);
  #endif

/* ---------------------
     Add source terms
   -------------------- */

  for (i = beg; i <= end; i++) {

    v   = state->vh[i];
    rhs = state->rhs[i];

    EXPAND(bx = v[BX1];  ,
           by = v[BX2];  ,
           bz = v[BX3];)

    #if BACKGROUND_FIELD == YES
     EXPAND(bx += bgf[i][BX1];  ,
            by += bgf[i][BX2];  ,
            bz += bgf[i][BX3];)
    #endif

    EXPAND (rhs[MX1] -= dt*divB[i]*bx;  ,
            rhs[MX2] -= dt*divB[i]*by;  ,
            rhs[MX3] -= dt*divB[i]*bz;)

    #if HAVE_ENERGY
     rhs[ENG] -= dt*v[BXn]*dpsi[i];
    #endif
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

double*** GLM_GetDivB

(

void

)

Definition at line 445 of file glm.c.

{
  return divB;
}

void GLM_Init	(	const Data *	d,
		const Time_Step *	Dts,
		Grid *	grid
	)

Initialize the maximum propagation speed glm_ch.

Definition at line 324 of file glm.c.

{
  int i,j,k,nv;
  double dxmin, gmaxc;

  #if PHYSICS == MHD 
   if (glm_ch < 0.0){
     double        *x1, *x2, *x3;
     static double **v, **lambda, *cs2;
     Index indx;

     if (v == NULL) {
       v      = ARRAY_2D(NMAX_POINT, NVAR, double);
       lambda = ARRAY_2D(NMAX_POINT, NVAR, double);
       cs2    = ARRAY_1D(NMAX_POINT, double);
     }
     x1 = grid[IDIR].x;
     x2 = grid[JDIR].x;
     x3 = grid[KDIR].x;

     glm_ch = 0.0;
     
  /* -- X1 - g_dir -- */

     g_dir = IDIR;  SetIndexes (&indx, grid);
     KDOM_LOOP(k) JDOM_LOOP(j){
       g_j = j; g_k = k;
       IDOM_LOOP(i) for (nv = NVAR; nv--;  ) v[i][nv] = d->Vc[nv][k][j][i];
       SoundSpeed2  (v, cs2, NULL, IBEG, IEND, CELL_CENTER, grid);
       Eigenvalues  (v, cs2, lambda, IBEG, IEND);
       IDOM_LOOP(i){
         glm_ch = MAX(glm_ch, fabs(lambda[i][KFASTP]));
         glm_ch = MAX(glm_ch, fabs(lambda[i][KFASTM]));
       }
     }

  /* -- X2 - g_dir -- */

     #if DIMENSIONS >= 2
     g_dir = JDIR;  SetIndexes (&indx, grid);
     KDOM_LOOP(k) IDOM_LOOP(i){
       g_i = i; g_k = k;
       JDOM_LOOP(j) for (nv = NVAR; nv--;  ) v[j][nv] = d->Vc[nv][k][j][i];
       SoundSpeed2  (v, cs2, NULL, JBEG, JEND, CELL_CENTER, grid);
       Eigenvalues  (v, cs2, lambda, JBEG, JEND);
       JDOM_LOOP(j){
         glm_ch = MAX(glm_ch, fabs(lambda[j][KFASTP]));
         glm_ch = MAX(glm_ch, fabs(lambda[j][KFASTM]));
       }
     }
     #endif

  /* -- X3 - g_dir -- */

     #if DIMENSIONS == 3
     g_dir = KDIR;  SetIndexes (&indx, grid);
     JDOM_LOOP(j) IDOM_LOOP(i){
       g_i = i; g_j = j;
       KDOM_LOOP(k) for (nv = NVAR; nv--;  ) v[k][nv] = d->Vc[nv][k][j][i];
       SoundSpeed2  (v, cs2, NULL, KBEG, KEND, CELL_CENTER, grid);
       Eigenvalues  (v, cs2, lambda, KBEG, KEND);
       KDOM_LOOP(k){
         glm_ch = MAX(glm_ch, fabs(lambda[k][KFASTP]));
         glm_ch = MAX(glm_ch, fabs(lambda[k][KFASTM]));
       }
     }
     #endif

   }
  #elif PHYSICS == RMHD 
   if (glm_ch < 0.0) glm_ch = 1.0;
  #endif

  #ifdef PARALLEL
   MPI_Allreduce (&glm_ch, &gmaxc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
   glm_ch = gmaxc;
  #endif
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GLM_Solve	(	const State_1D *	state,
		double **	VL,
		double **	VR,
		int	beg,
		int	end,
		Grid *	grid
	)

Solve the 2x2 linear hyperbolic GLM-MHD system given by the divergence cleaning approach. Build new states VL and VR for Riemann problem. We use Eq. (42) of Dedner et al (2002)

Parameters

[in,out]	state	pointer to a State_1D structure
[out]	VL	left-interface state to be passed to the Riemann solver
[out]	VR	right-interface state to be passed to the Riemann solver
[in]	beg	starting index of computation
[in]	end	final index of computation
[in]	grid	pointer to array of Grid structures

The purpose of this function is two-fold:

1. assign a unique value to the normal component of magnetic field and to te scalar psi before the actual Riemann solver is called.
2. compute GLM fluxes in Bn and psi.

The following MAPLE script has been used

restart;
with(linalg);
A := matrix(2,2,[ 0, 1, c^2, 0]);

R := matrix(2,2,[ 1, 1, c, -c]);
L := inverse(R);
E := matrix(2,2,[c,0,0,-c]);
multiply(R,multiply(E,L));

Definition at line 24 of file glm.c.

{
  int    i, nv, nflag=0;
  double Bm, psim;
  double dB, dpsi;
  double **v;
  static double *wp, *wm, *Bxp, *Bxm;

  if (wp == NULL){
    wp = ARRAY_1D(NMAX_POINT, double);
    wm = ARRAY_1D(NMAX_POINT, double);
    Bxp = ARRAY_1D(NMAX_POINT, double);
    Bxm = ARRAY_1D(NMAX_POINT, double);
  }

  #ifdef CTU 
   #if PARABOLIC_FLUX != EXPLICIT
    if (g_intStage == 1) nflag = 1;
   #endif

  /* ------------------------------------------------------
      The nflag = 1 is not necessary but it improves the 
      results when dimensionally unsplit CTU is used. 
      In practice, it redefines the input normal field 
      Bn to be at time level n rather than n+1/2.
     ------------------------------------------------------- */

   if (nflag){
     v = state->v;
     for (i = beg-1; i <= end+1; i++){
       dB   = v[i+1][BXn] - v[i][BXn];
       dpsi = v[i+1][PSI_GLM] - v[i][PSI_GLM];
       wm[i] =  0.5*(dB - dpsi/glm_ch);
       wp[i] =  0.5*(dB + dpsi/glm_ch);
     }
     for (i = beg; i <= end+1; i++){
       dB  = MC(wm[i], wm[i-1]);
       dB += MC(wp[i], wp[i-1]);
       Bxp[i] = v[i][BXn] + 0.5*dB;
       Bxm[i] = v[i][BXn] - 0.5*dB;
     }
   }  
  #endif

/* -------------------------------------------------
    Solve the Riemann problem for the 2x2 linear 
    GLM system. Use the solution to assign a single 
    value to both left and right (Bn,psi).
   ------------------------------------------------- */

  for (i = beg; i <= end; i++){
    for (nv = NVAR; nv--;   ){
      VL[i][nv] = state->vL[i][nv];
      VR[i][nv] = state->vR[i][nv];
    }

    #ifdef CTU
     if (nflag){
       VL[i][BXn] = Bxp[i];
       VR[i][BXn] = Bxm[i+1];
     }
    #endif

    dB   = VR[i][BXn]     - VL[i][BXn];
    dpsi = VR[i][PSI_GLM] - VL[i][PSI_GLM];

    Bm   = 0.5*(VL[i][BXn]     + VR[i][BXn])     - 0.5*dpsi/glm_ch;
    psim = 0.5*(VL[i][PSI_GLM] + VR[i][PSI_GLM]) - 0.5*glm_ch*dB;

    state->bn[i] = VL[i][BXn] = VR[i][BXn] = Bm;
    VL[i][PSI_GLM] = VR[i][PSI_GLM] = psim;
  }

  #if COMPUTE_DIVB == YES
   if (g_intStage == 1) GLM_ComputeDivB(state, grid);  
  #endif

}

Here is the call graph for this function:

Here is the caller graph for this function:

void GLM_SolveNEW	(	const State_1D *	state,
		int	beg,
		int	end,
		Grid *	grid
	)

void GLM_Source	(	const Data_Arr	Q,
		double	dt,
		Grid *	grid
	)

Include the parabolic source term of the Lagrangian multiplier equation in a split fashion for the mixed GLM formulation. Ref. Mignone & Tzeferacos, JCP (2010) 229, 2117, Equation (27).

Definition at line 139 of file glm.c.

{
  int    i,j,k;
  double cr, cp, scrh;
  double dx, dy, dz, dtdx;

  #ifdef CHOMBO
   dtdx = dt;
  #else
   dx   = grid[IDIR].dx[IBEG];
   dtdx = dt/dx;
  #endif

  cp   = sqrt(dx*glm_ch/GLM_ALPHA);
  scrh = dtdx*glm_ch*GLM_ALPHA; /* CFL*g_inputParam[ALPHA]; */

  scrh = exp(-scrh);
  DOM_LOOP(k,j,i) Q[PSI_GLM][k][j][i] *= scrh;
  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

double glm_ch

The propagation speed of divergence error.

Definition at line 21 of file glm.c.