PLUTO
plm_states.c File Reference

Piecewise linear reconstruction. More...

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

Go to the source code of this file.

Functions

static void MonotonicityTest (double **, double **, double **, int, int)
 
static void FourthOrderLinear (const State_1D *, int, int, Grid *)
 
void States (const State_1D *state, int beg, int end, Grid *grid)
 

Detailed Description

Piecewise linear reconstruction.

Compute interface states using piecewise linear reconstruction inside each zone. Reconstruction is performed in primitive variable when CHAR_LIMITING == NO ) or characteristic variables when CHAR_LIMITING == YES .

The convention used throughout is that vL and vR are left and right states with respect to the interface, while vm and vp refer to the cell center:

                    vL(i)-> <-vR(i)
     |----------*----------|----------*----------|
      <-vm(i)  (i)  vp(i)->         (i+1)

The default setting (LIMITER == DEFAULT) applies a different limiter to each variable:

  • MC for density
  • VanLeer for velcoity and magnetic field
  • MinMod for pressure.

Otherwise the same limiter can be imposed to all variables from definitions.h.

The approach followed here is taken from Mignone (JCP, 2014) "High-order conservative reconstruction schemes for finite volume methods in cylindrical and spherical geometries" where interface states are constructed as

\[ V_i^{\pm} = V_i \;\pm\; \overline{\Delta V}_i \; d^\pm_i\,, \qquad \overline{\Delta V}_i = {\rm Lim}\left(\Delta V_i^+,\,\Delta V_i^-,\, c^+_i, c^-_i\right)\,, \qquad \Delta V^\pm_i = \pm\left(V_{i\pm1}-V_i\right)w^\pm_i \]

where

\[ w^\pm_i = \left|\frac{\Delta\xi_i}{\bar{\xi}_{i\pm1} - \bar{\xi}_i}\right| \,, \qquad d^\pm_i = \left|\frac{\xi_{i\pm\HALF} - \bar{\xi}_i}{\Delta\xi_i}\right|\,, \qquad c^\pm_i = \left|\frac{\bar{\xi}_{i\pm1} - \bar{\xi}_i} {\xi_{i\pm\HALF}- \bar{\xi}_i}\right| \]

are interpolation coefficients returned by the PLM_Coeffs structure (see the PLM_CoefficientsGet() function). The slope limiter $\overline{\Delta V}_i=$ dv_lim[i][nv] is a function of the forward and backward (undivided) derivatives and geometrical coefficients in case of non-uniform or non-Cartesian grids. Different limiters are implemented through small macros defined in plm_coeffs.h. Note that $ \bar{\xi} $ defines the centroid of volume which differs from the cell center in cylindrical and spherical geometires.

A stencil of 3 zones is required for all limiters except for the FOURTH_ORDER_LIM which requires 5 zones.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
June 16, 2015

References

  • "High-order conservative reconstruction schemes for finite volume methods in cylindrical and spherical coordinates" A. Mignone, JCP (2014), 270, 784.

Definition in file plm_states.c.

Function Documentation

void FourthOrderLinear ( const State_1D state,
int  beg,
int  end,
Grid grid 
)
static

Definition at line 301 of file plm_states.c.

318 {
319  int i, nv;
320  static double **s;
321  static double **dv, **dvf, **dvc, **dvlim;
322  double scrh, dvp, dvm, dvl;
323  double **v, **vp, **vm;
324 
325  v = state->v;
326  vp = state->vp;
327  vm = state->vm;
328 
329  if (s == NULL){
330  s = ARRAY_2D(NMAX_POINT, NVAR, double);
331  dv = ARRAY_2D(NMAX_POINT, NVAR, double);
332  dvf = ARRAY_2D(NMAX_POINT, NVAR, double);
333  dvc = ARRAY_2D(NMAX_POINT, NVAR, double);
334  dvlim = ARRAY_2D(NMAX_POINT, NVAR, double);
335  }
336 
337  #if TIME_STEPPING == HANCOCK && PHYSICS != RMHD
338  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
339  #endif
340 /*
341  #if GEOMETRY != CARTESIAN
342  print1 ("! FourthOrderLinear: only Cartesian geometry supported\n");
343  QUIT_PLUTO(1);
344  #endif
345 */
346 /* -----------------------------------------------------------
347  compute undivided differences
348  ----------------------------------------------------------- */
349 
350  for (i = beg-2; i <= end+1; i++){
351  for (nv = 0; nv < NVAR; nv++) dv[i][nv] = v[i+1][nv] - v[i][nv];
352  }
353 
354  for (i = beg - 1; i <= end + 1; i++){
355  for (nv = 0; nv < NVAR; nv++){
356  dvp = dv[i][nv]; dvm = dv[i-1][nv];
357  dvc[i][nv] = 0.5*(dvp + dvm);
358  s[i][nv] = (dvp > 0.0 ? 0.5:-0.5)
359  + (dvm > 0.0 ? 0.5:-0.5);
360  dvlim[i][nv] = 2.0*MIN(fabs(dvp), fabs(dvm));
361  dvf[i][nv] = MIN(fabs(dvc[i][nv]), dvlim[i][nv])*s[i][nv];
362  }}
363 
364  for (i = beg; i <= end; i++){
365  for (nv = 0; nv < NVAR; nv++){
366  if (dv[i][nv]*dv[i-1][nv] > 0.0) {
367  scrh = 4.0/3.0*dvc[i][nv] - (dvf[i+1][nv] + dvf[i-1][nv])/6.0;
368  dvlim[i][nv] = MIN(fabs(scrh), dvlim[i][nv])*s[i][nv];
369  }else{
370  dvlim[i][nv] = 0.0;
371  }
372  }
373 
374  for (nv = NVAR; nv--; ) {
375  vp[i][nv] = v[i][nv] + 0.5*dvlim[i][nv];
376  vm[i][nv] = v[i][nv] - 0.5*dvlim[i][nv];
377  }
378  #if (PHYSICS == RHD || PHYSICS == RMHD)
379  VelocityLimiter(v[i], vp[i], vm[i]);
380  #endif
381  }
382 
383 /* -------------------------------------------
384  Shock flattening
385  ------------------------------------------- */
386 
387  #if SHOCK_FLATTENING == MULTID && CHAR_LIMITING == NO
388  Flatten (state, beg, end, grid);
389  #endif
390 
391 /* -------------------------------------------
392  Shock flattening
393  ------------------------------------------- */
394 
395  #if SHOCK_FLATTENING == ONED
396  Flatten (state, beg, end, grid);
397  #endif
398 
399 /* -------------------------------------------
400  Assign face-centered magnetic field
401  ------------------------------------------- */
402 
403  #ifdef STAGGERED_MHD
404  for (i = beg - 1; i <= end; i++) {
405  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
406  }
407  #endif
408 
409 /* --------------------------------------------------------
410  evolve center values by dt/2
411  -------------------------------------------------------- */
412 
413 
414  #if TIME_STEPPING == HANCOCK
415  HancockStep(state, beg, end, grid);
416  #endif
417 
418 /* -------------------------------------------
419  compute states in conservative variables
420  ------------------------------------------- */
421 
422  PrimToCons (state->vp, state->up, beg, end);
423  PrimToCons (state->vm, state->um, beg, end);
424 }
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
void VelocityLimiter(double *, double *, double *)
Definition: vel_limiter.c:16
int end
Global end index for the local array.
Definition: structs.h:116
void Flatten(const State_1D *, int, int, Grid *)
Definition: flatten.c:4
tuple scrh
Definition: configure.py:200
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
double * h
Enthalpy.
Definition: structs.h:159
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
Definition: eos.c:16
int BXn
Definition: globals.h:75
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
Definition: structs.h:142
#define MIN(a, b)
Definition: macros.h:104
#define CELL_CENTER
Definition: pluto.h:205
int beg
Global start index for the local array.
Definition: structs.h:115
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
#define s
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
Definition: mappers.c:26
int i
Definition: analysis.c:2
double * bn
Face magentic field, bn = bx(i+1/2)
Definition: structs.h:165
double ** um
same as vm, in conservative vars
Definition: structs.h:146
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
double * a2
Sound speed squared.
Definition: structs.h:158
double ** up
same as vp, in conservative vars
Definition: structs.h:147
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
#define NVAR
Definition: pluto.h:609
void HancockStep(const State_1D *, int, int, Grid *)
Definition: hancock.c:136

Here is the call graph for this function:

Here is the caller graph for this function:

void MonotonicityTest ( double **  v,
double **  vp,
double **  vm,
int  beg,
int  end 
)
static

Definition at line 679 of file plm_states.c.

683 {
684  int i, nv;
685  double vp_max, vp_min, vm_max, vm_min;
686 
687  for (i = beg; i <= end; i++){ VAR_LOOP(nv) {
688  vp_max = MAX(v[i][nv], v[i+1][nv]) + 1.e-9;
689  vp_min = MIN(v[i][nv], v[i+1][nv]) - 1.e-9;
690 
691  vm_max = MAX(v[i][nv], v[i-1][nv]) + 1.e-9;
692  vm_min = MIN(v[i][nv], v[i-1][nv]) - 1.e-9;
693 
694  if (vp[i][nv] > vp_max || vp[i][nv] < vp_min){
695  print ("! States: monotonicity violated at + interface, i = %d\n",i);
696  print ("vp = %12.6e; vmax = %12.6e, vmin = %12.6e\n",
697  vp[i][nv], vp_max, vp_min);
698  QUIT_PLUTO(1);
699  }
700 
701  if (vm[i][nv] > vm_max || vm[i][nv] < vm_min){
702  print ("! States: monotonicity violated at - interface, i = %d\n",i);
703  print ("vm = %12.6e; vmax = %12.6e, vmin = %12.6e\n",
704  vm[i][nv], vm_max, vm_min);
705  QUIT_PLUTO(1);
706  }
707  }}
708 }
#define MAX(a, b)
Definition: macros.h:101
int end
Global end index for the local array.
Definition: structs.h:116
#define MIN(a, b)
Definition: macros.h:104
#define VAR_LOOP(n)
Definition: macros.h:226
int beg
Global start index for the local array.
Definition: structs.h:115
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
int i
Definition: analysis.c:2
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

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

Compute states using piecewise linear interpolation.

Parameters
[in]statepointer to a State_1D structure
[in]begstarting point where vp and vm must be computed
[in]endfinal point where vp and vm must be computed
[in]gridpointer to array of Grid structures
Returns
This function has no return value.

Compute 1D left and right interface states using piecewise linear reconstruction and the characteristic decomposition of the quasi-linear form of the equations.

This is done by first extrapolating the cell center value to the interface using piecewise limited linear reconstruction on the characteristic variables.

Left and right states are then evolved for the half time step using characteristic tracing if necessary.

Definition at line 80 of file plm_states.c.

92 {
93  int nv, i;
94  double **v, **vp, **vm;
95  double dv_lim[NVAR], dvp[NVAR], dvm[NVAR];
96  double cp, cm, wp, wm, dp, dm;
97  PLM_Coeffs plm_coeffs;
98  static double **dv;
99 
100  #if LIMITER == FOURTH_ORDER_LIM
101  FourthOrderLinear(state, beg, end, grid);
102  return;
103  #endif
104 
105 /* -----------------------------------------------------------
106  0. Memory allocation and pointer shortcuts, geometrical
107  coefficients and conversion to 4vel (if required)
108  ----------------------------------------------------------- */
109 
110  if (dv == NULL) {
111  dv = ARRAY_2D(NMAX_POINT, NVAR, double);
112  }
113 
114  v = state->v;
115  vp = state->vp;
116  vm = state->vm;
117 
118 #if UNIFORM_CARTESIAN_GRID == NO
119  PLM_CoefficientsGet (&plm_coeffs, g_dir);
120 #endif
121 
122 #if RECONSTRUCT_4VEL
123  ConvertTo4vel (state->v, beg-1, end+1);
124 #endif
125 
126 /* -------------------------------------------
127  compute undivided differences
128  ------------------------------------------- */
129 
130  for (i = beg-1; i <= end; i++){
131  NVAR_LOOP(nv) dv[i][nv] = v[i+1][nv] - v[i][nv];
132  }
133 
134 /* -------------------------------------------
135  2. Main spatial loop
136  ------------------------------------------- */
137 
138  for (i = beg; i <= end; i++){
139 
140  /* ---------------------------------------------------------
141  2a. compute forward (dvp) and backward (dvm) derivatives
142  --------------------------------------------------------- */
143 
144  #if UNIFORM_CARTESIAN_GRID == YES
145  cp = cm = 2.0;
146  wp = wm = 1.0;
147  dp = dm = 0.5;
148  NVAR_LOOP(nv) {
149  dvp[nv] = dv[i][nv];
150  dvm[nv] = dv[i-1][nv];
151  }
152  #else
153  cp = plm_coeffs.cp[i]; cm = plm_coeffs.cm[i];
154  wp = plm_coeffs.wp[i]; wm = plm_coeffs.wm[i];
155  dp = plm_coeffs.dp[i]; dm = plm_coeffs.dm[i];
156  NVAR_LOOP(nv) {
157  dvp[nv] = dv[i][nv]*wp;
158  dvm[nv] = dv[i-1][nv]*wm;
159  }
160  #endif
161 
162  /* ---------------------------------------------------------------
163  2b. if shock-flattening is enabled, revert to minmod limiter
164  --------------------------------------------------------------- */
165 
166 #if SHOCK_FLATTENING == MULTID
167  if (state->flag[i] & FLAG_FLAT) {
168  NVAR_LOOP(nv)vp[i][nv] = vm[i][nv] = v[i][nv];
169  continue;
170  }else if (state->flag[i] & FLAG_MINMOD) {
171  NVAR_LOOP(nv){
172  SET_MM_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
173  vp[i][nv] = v[i][nv] + dv_lim[nv]*dp;
174  vm[i][nv] = v[i][nv] - dv_lim[nv]*dm;
175  }
176  #if (PHYSICS == RHD || PHYSICS == RMHD)
177  VelocityLimiter (v[i], vp[i], vm[i]);
178  #endif
179  continue;
180  }
181 #endif
182 
183  /* ---------------------------------------------------------
184  2c. DEFAULT limiting: combination of different limiters
185  (this has some hystorical reasons)
186  --------------------------------------------------------- */
187 
188  #if LIMITER == DEFAULT
189  SET_MC_LIMITER(dv_lim[RHO], dvp[RHO], dvm[RHO], cp, cm);
190  #if PHYSICS != ADVECTION
191  EXPAND(SET_VL_LIMITER(dv_lim[VX1], dvp[VX1], dvm[VX1], cp, cm); ,
192  SET_VL_LIMITER(dv_lim[VX2], dvp[VX2], dvm[VX2], cp, cm); ,
193  SET_VL_LIMITER(dv_lim[VX3], dvp[VX3], dvm[VX3], cp, cm);)
194  #endif
195 
196  #if PHYSICS == MHD || PHYSICS == RMHD
197  EXPAND(SET_VL_LIMITER(dv_lim[BX1], dvp[BX1], dvm[BX1], cp, cm); ,
198  SET_VL_LIMITER(dv_lim[BX2], dvp[BX2], dvm[BX2], cp, cm); ,
199  SET_VL_LIMITER(dv_lim[BX3], dvp[BX3], dvm[BX3], cp, cm);)
200  #ifdef GLM_MHD
201  SET_MC_LIMITER(dv_lim[PSI_GLM], dvp[PSI_GLM], dvm[PSI_GLM], cp, cm);
202  #endif
203  #endif
204 
205  #if HAVE_ENERGY
206  SET_MM_LIMITER(dv_lim[PRS], dvp[PRS], dvm[PRS], cp, cm);
207  #endif
208 
209  #if NFLX != NVAR /* -- scalars: MC lim -- */
210  for (nv = NFLX; nv < NVAR; nv++){
211  SET_MC_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
212  }
213  #endif
214  #endif /* LIMITER == DEFAULT */
215 
216  /* ----------------------------------------
217  2d. construct (+) and (-) states
218  ---------------------------------------- */
219 
220  for (nv = 0; nv < NVAR; nv++){
221  #if LIMITER != DEFAULT /* -- same limiter for all variables -- */
222  SET_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
223  #endif
224 
225  vp[i][nv] = v[i][nv] + dv_lim[nv]*dp;
226  vm[i][nv] = v[i][nv] - dv_lim[nv]*dm;
227  }
228 
229  /* --------------------------------------
230  2e. Relativistic Limiter
231  -------------------------------------- */
232 
233  #if (PHYSICS == RHD || PHYSICS == RMHD)
234  VelocityLimiter (v[i], vp[i], vm[i]);
235  #endif
236  } /* -- end loop on zones -- */
237 
238 /* ----------------------------------------------------
239  3a. Check monotonicity
240  ---------------------------------------------------- */
241 
242 #if CHECK_MONOTONICITY == YES
243  MonotonicityTest(state->v, state->vp, state->vm, beg, end);
244 #endif
245 
246 /* -------------------------------------------
247  3b. Shock flattening
248  ------------------------------------------- */
249 
250 #if SHOCK_FLATTENING == ONED
251  Flatten (state, beg, end, grid);
252 #endif
253 
254 /* -------------------------------------------
255  4. Assign face-centered magnetic field
256  ------------------------------------------- */
257 
258 #ifdef STAGGERED_MHD
259  for (i = beg - 1; i <= end; i++) {
260  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
261  }
262 #endif
263 
264 /* --------------------------------------------------------
265  5. Evolve L/R states and center value by dt/2
266  -------------------------------------------------------- */
267 
268 #if TIME_STEPPING == CHARACTERISTIC_TRACING
269  CharTracingStep(state, beg, end, grid);
270 #elif TIME_STEPPING == HANCOCK
271  HancockStep(state, beg, end, grid);
272  /* The conservative hancock scheme is used only for RMHD with linear
273  interpolation in primitive variables.
274  In this case, the conversion to 3-vel is already carried out in
275  the predictor step and we need to skip it here. */
276  #if PRIMITIVE_HANCOCK == NO
277  PrimToCons (state->vp, state->up, beg, end);
278  PrimToCons (state->vm, state->um, beg, end);
279  return;
280  #endif
281 #endif
282 
283 /* --------------------------------------------------------
284  6. Convert back to 3-velocity
285  -------------------------------------------------------- */
286 
287 #if RECONSTRUCT_4VEL
288  ConvertTo3vel (state->v, beg-1, end+1);
289  ConvertTo3vel (state->vp, beg, end);
290  ConvertTo3vel (state->vm, beg, end);
291 #endif
292 
293 /* ----------------------------------------------
294  7. Obtain L/R states in conservative variables
295  ---------------------------------------------- */
296 
297  PrimToCons (state->vp, state->up, beg, end);
298  PrimToCons (state->vm, state->um, beg, end);
299 }
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
void VelocityLimiter(double *, double *, double *)
Definition: vel_limiter.c:16
void ConvertTo4vel(double **, int, int)
Definition: four_vel.c:4
int end
Global end index for the local array.
Definition: structs.h:116
#define VX2
Definition: mod_defs.h:29
#define SET_VL_LIMITER(dv, dvp, dvm, cp, cm)
Definition: plm_coeffs.h:114
void Flatten(const State_1D *, int, int, Grid *)
Definition: flatten.c:4
#define FLAG_FLAT
Reconstruct using FLAT limiter.
Definition: pluto.h:180
#define RHO
Definition: mod_defs.h:19
void CharTracingStep(const State_1D *, int, int, Grid *)
Definition: char_tracing.c:198
double * cm
Definition: plm_coeffs.h:40
#define PSI_GLM
Definition: mod_defs.h:34
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
#define VX1
Definition: mod_defs.h:28
#define FLAG_MINMOD
Reconstruct using MINMOD limiter.
Definition: pluto.h:179
int BXn
Definition: globals.h:75
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
Definition: structs.h:142
#define NFLX
Definition: mod_defs.h:32
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
int beg
Global start index for the local array.
Definition: structs.h:115
#define NVAR_LOOP(n)
Definition: pluto.h:618
#define SET_MC_LIMITER(dv, dvp, dvm, cp, cm)
Definition: plm_coeffs.h:119
#define PHYSICS
Definition: definitions_01.h:1
if(divB==NULL)
Definition: analysis.c:10
unsigned char * flag
Definition: structs.h:168
double * cp
Definition: plm_coeffs.h:39
void ConvertTo3vel(double **, int, int)
Definition: four_vel.c:25
#define SET_MM_LIMITER(dv, dvp, dvm, cp, cm)
Definition: plm_coeffs.h:76
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
#define MHD
Definition: pluto.h:111
#define BX3
Definition: mod_defs.h:27
void PLM_CoefficientsGet(PLM_Coeffs *plm_coeffs, int dir)
Definition: plm_coeffs.c:79
double * dm
Definition: plm_coeffs.h:44
double * dp
Definition: plm_coeffs.h:43
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
Definition: mappers.c:26
int i
Definition: analysis.c:2
double * bn
Face magentic field, bn = bx(i+1/2)
Definition: structs.h:165
static void MonotonicityTest(double **, double **, double **, int, int)
Definition: plm_states.c:679
#define BX1
Definition: mod_defs.h:25
double ** um
same as vm, in conservative vars
Definition: structs.h:146
#define VX3
Definition: mod_defs.h:30
double * wm
Definition: plm_coeffs.h:42
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
double ** up
same as vp, in conservative vars
Definition: structs.h:147
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
#define BX2
Definition: mod_defs.h:26
#define NVAR
Definition: pluto.h:609
#define RMHD
Definition: pluto.h:112
static void FourthOrderLinear(const State_1D *, int, int, Grid *)
Definition: plm_states.c:301
#define GLM_MHD
Definition: glm.h:25
double * wp
Definition: plm_coeffs.h:41
void HancockStep(const State_1D *, int, int, Grid *)
Definition: hancock.c:136

Here is the call graph for this function: