sb_flux.c File Reference

Enforce conservation at the X1 boundaries in the shearing-box module. More...

#include "pluto.h"

Include dependency graph for sb_flux.c:

Go to the source code of this file.

Macros
#define	NVLAST   (NVAR-1)
#define	swL   0.0
#define	swR   (1.0 - swL)

Functions
void	SB_SaveFluxes (State_1D *state, Grid *grid)
void	SB_CorrectFluxes (Data_Arr U, double t, double dt, Grid *grid)

Variables
static double ***	FluxL
	Array of fluxes at the left x-boundary. More...
static double ***	FluxR
	Array of fluxes at the right x-boundary. More...

Detailed Description

Enforce conservation at the X1 boundaries in the shearing-box module.

This file provides functions to store and then modify the upwind fluxes computed during the Riemann solver at the leftmost and righmost physical boundaries. These tasks are performed only when either SB_SYMMETRIZE_HYDRO, SB_SYMMETRIZE_EY, SB_SYMMETRIZE_EZ flags are set to YES in Src/MHD/shearingbox.h and are useful to avoid loss of conservation in the hydrodynamical variables (density, momentum, energy) and/or magnetic fields.

This is first achieved by calling the SB_SaveFluxes() function during the time stepping scheme, with the purpose of storing the leftmost and rightmost conservative fluxes in the x-direction into the static arrays FluxL[] and FluxR[].

These fluxes are then subsequently used by SB_CorrectFluxes() which interpolates the fluxes and properly correct leftmost and rightmost cell-centered flow quantities to ensure conservation.

The treatment of staggered magnetic field is done similarly by SB_CorrectEMF().

References

• "??"
  Mignone et al, in preparation

Authors

A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
G. Muscianisi (g.mus.nosp@m.ican.nosp@m.isi@c.nosp@m.inec.nosp@m.a.it)
G. Bodo (bodo@.nosp@m.oato.nosp@m..inaf.nosp@m..it)

Date

Feb 14, 2014

Definition in file sb_flux.c.

Macro Definition Documentation

#define NVLAST   (NVAR-1)

Definition at line 44 of file sb_flux.c.

#define swL   0.0

Sets the weight coefficient in the range [0,1] used during flux symmetrization at the left and right boundaries,

FL -> swL*FL    + swR*I(FR)
FR -> swL*I(FL) + swR*FR 

where swR=1-swL and I() stands for conservative interpolation.

Definition at line 53 of file sb_flux.c.

#define swR   (1.0 - swL)

Definition at line 54 of file sb_flux.c.

Function Documentation

void SB_CorrectFluxes	(	Data_Arr	U,
		double	t,
		double	dt,
		Grid *	grid
	)

Interpolate x-fluxes FLuxL and FluxR and properly correct leftmost and rightmost cells to ensure conservation.

Parameters

[in,out]	U	data array containing cell-centered quantities
[in]	t	the time step at which fluxes have been computed
[in]	dt	the time step being used in the integrator stage
[in]	grid	pointer to array of Grid structures

Note

Modify the values of FluxR (on the left) and FluxL (on the right) to properly account for the shift induced by the moving boundaries. For safety reason (avoid overwriting arrays when a processor owns both the left and the right side) we first copy, then interpolate and finally average:

• On the left boundary:

  FluxR -> FluxLR -> I(FluxLR) -> fL = wL*FluxL + wR*FluxLR

• On the right boundary:

  FluxL -> FluxLR -> I(FluxLR) -> fR = wL*FluxLR + wR*FluxR

  ---

Definition at line 105 of file sb_flux.c.

{
  int    i, j, k, nv;
  double dtdx;
  static double **fL, **fR;
  static double ****Ftmp;
  
  #if SB_SYMMETRIZE_HYDRO == NO
   return;
  #endif

  if (fL == NULL){
    fL   = ARRAY_2D(NVAR, NMAX_POINT, double);
    fR   = ARRAY_2D(NVAR, NMAX_POINT, double);
    Ftmp = ARRAY_4D(NVAR, NX3_TOT, NX2_TOT, 1, double);
  }

  dtdx = dt/grid[IDIR].dx[IBEG];

{
t = g_time;
#if TIME_STEPPING == RK2
 if (g_intStage == 2) t += g_dt;
#elif TIME_STEPPING == RK3
 if (g_intStage == 2) t += 0.5*g_dt;
 if (g_intStage == 3) t += g_dt;
#endif
#ifdef CTU
 if (g_intStage == 2) t += 0.5*g_dt;
#endif
}

/* -----------------------------------------------------------------------
    Exchange left and right fluxes if they were initially computed on
    different processors.
    After this step, both the left- and right- processor- will share
    the same fluxes FluxL[] and FluxR[].
   ----------------------------------------------------------------------- */

  #ifdef PARALLEL
   ExchangeX (FluxL[0][0], FluxR[0][0], NVAR*NX2_TOT*NX3_TOT, grid);
  #endif

/* --------------------------------------------------------------------- */
/*! \note 
    Modify the values of FluxR (on the left) and FluxL (on the right)
    to properly account for the shift induced by the moving boundaries.
    For safety reason (avoid overwriting arrays when a processor owns 
    both the left and the right side) we first copy, then interpolate and
    finally average:

    - On the left boundary:

      FluxR -> FluxLR -> I(FluxLR) -> fL = wL*FluxL + wR*FluxLR

    - On the right boundary:

      FluxL -> FluxLR -> I(FluxLR) -> fR = wL*FluxLR + wR*FluxR
      
   --------------------------------------------------------------------- */

  if (grid[IDIR].lbound != 0){  /* ---- Left x-boundary ---- */

    RBox box;

  /* -- set grid ranges of FluxR and exchange b.c. -- */

    box.ib = 0; box.ie = 0; 
    box.jb = 0; box.je = NX2_TOT-1; 
    box.kb = 0; box.ke = NX3_TOT-1;

  /* -- copy FluxR on temporary storage and shift -- */

    for (nv = 0; nv <= NVLAST; nv++) {
      BOX_LOOP((&box), k, j, i) Ftmp[nv][k][j][i] = FluxR[nv][k][j];
      SB_SetBoundaryVar(Ftmp[nv], &box, X1_BEG, t, grid);
    }

  /* -- symmetrize fluxes -- */
  
    for (k = KBEG; k <= KEND; k++){
      for (nv = 0; nv <= NVLAST; nv++){
      for (j = JBEG; j <= JEND; j++){
        fL[nv][j] = swL*FluxL[nv][k][j] + swR*Ftmp[nv][k][j][0];
      }}
      for (j = JBEG; j <= JEND; j++){
        #if HAVE_ENERGY
         fL[ENG][j] += swR*sb_vy*(fL[MX2][j] + 0.5*sb_vy*fL[RHO][j]);
        #endif
        #ifdef GLM_MHD
         fL[BX2][j] -= swR*sb_vy*fL[PSI_GLM][j]/glm_ch/glm_ch;
/*         fL[BX2][j] -= 0.5*sb_vy*FluxL[PSI_GLM][k][j]/glm_ch/glm_ch; */
        #endif
        fL[MX2][j] += swR*sb_vy*fL[RHO][j];
        
      /* -- update solution vector -- */

        for (nv = 0; nv <= NVLAST; nv++){
          U[k][j][IBEG][nv] += dtdx*(fL[nv][j] - FluxL[nv][k][j]); 
        } 
      }
    }
  }   

  if (grid[IDIR].rbound != 0){  /* ---- Right x-boundary ---- */

  /* -- set grid ranges of FluxL and exchange b.c. -- */

    RBox box;
    box.ib = 0; box.ie = 0;
    box.jb = 0; box.je = NX2_TOT-1;
    box.kb = 0; box.ke = NX3_TOT-1;

  /* -- copy FluxL on temporary storage and shift -- */

    for (nv = 0; nv <= NVLAST; nv++) {
      BOX_LOOP((&box), k, j, i) Ftmp[nv][k][j][i] = FluxL[nv][k][j];
      SB_SetBoundaryVar(Ftmp[nv], &box, X1_END, t, grid);
    }

  /* -- symmetrize fluxes -- */
  
    for (k = KBEG; k <= KEND; k++){
      for (nv = 0; nv <= NVLAST; nv++){
      for (j = JBEG; j <= JEND; j++){
        fR[nv][j] = swL*Ftmp[nv][k][j][0] + swR*FluxR[nv][k][j];
      }}
      for (j = JBEG; j <= JEND; j++){
        #if HAVE_ENERGY
         fR[ENG][j] += swL*sb_vy*(-fR[MX2][j] + 0.5*sb_vy*fR[RHO][j]);
        #endif
        #ifdef GLM_MHD
         fR[BX2][j] += swL*sb_vy*fR[PSI_GLM][j]/glm_ch/glm_ch;
/*         fR[BX2][j] += 0.5*sb_vy*FluxR[PSI_GLM][k][j]/glm_ch/glm_ch; */
        #endif
        fR[MX2][j] -= swL*sb_vy*fR[RHO][j];

        for (nv = 0; nv <= NVLAST; nv++){
          U[k][j][IEND][nv] -= dtdx*(fR[nv][j] - FluxR[nv][k][j]); 
        }
      }
    }  
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void SB_SaveFluxes	(	State_1D *	state,
		Grid *	grid
	)

Store leftmost and rightmost conservative fluxes in the x direction into FluxL and FluxR.

Parameters

[in]	state	pointer to a State_1D structure
[in]	grid	pointer to an array of Grid structures

Returns

This function has no return value.

Definition at line 57 of file sb_flux.c.

{
  int nv;

  #if SB_SYMMETRIZE_HYDRO == NO
   return;
  #endif

  if (FluxL == NULL){
    FluxL = ARRAY_3D(NVAR, NX3_TOT, NX2_TOT, double);
    FluxR = ARRAY_3D(NVAR, NX3_TOT, NX2_TOT, double);
  }

/* ---------------------------------------------------- 
    Save Fluxes on the LEFT physical boundary
   ---------------------------------------------------- */
  
  if (g_dir == IDIR && grid[IDIR].lbound != 0){
    state->flux[IBEG - 1][MX1] += state->press[IBEG - 1];
    state->press[IBEG - 1]      = 0.0;  
    for (nv = 0; nv <= NVLAST; nv++){
      FluxL[nv][g_k][g_j] = state->flux[IBEG - 1][nv];
    } 
  }

/* ---------------------------------------------------- 
    Save Fluxes on the RIGHT pysical boundary
   ---------------------------------------------------- */
  
  if (g_dir == IDIR && grid[IDIR].rbound != 0){
    state->flux[IEND][MX1] += state->press[IEND];    
    state->press[IEND]      = 0.0;  
    for (nv = 0; nv <= NVLAST; nv++){
      FluxR[nv][g_k][g_j] = state->flux[IEND][nv];
    }
  }
}

Here is the caller graph for this function:

Variable Documentation

double*** FluxL

static

Array of fluxes at the left x-boundary.

Definition at line 41 of file sb_flux.c.

double*** FluxR

static

Array of fluxes at the right x-boundary.

Definition at line 42 of file sb_flux.c.