PLUTO
Macros | Functions
ppm_states.c File Reference

Piecewise parabolic reconstruction. More...

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

Go to the source code of this file.

Macros

#define PARABOLIC_LIM   1
 

Functions

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

Detailed Description

Piecewise parabolic reconstruction.

Compute interface states using piecewise parabolic reconstruction inside each zone. Reconstruction is performed in primitive variable (when CHAR_LIMITING == NO) or characteristic variables (when CHAR_LIMITING == YES). The reconstruction process follows the following steps:

  • compute unlimited 4-th or 5-th order interface values from within the zone ( ${\rm vp} = v_{i,+},\; {\rm vm} = v_{i,-} $).
  • interface states must lie between adjacent cell averages. This is achieved by working on the increments,

    \[ \left\{\begin{array}{lll} \min(v_{i+1},v_i) \le v_{i,+} \le \max(v_{i+1},v_i) & \qquad \Longrightarrow\qquad & \delta v_{i,+} \to {\rm minmod}(\delta v_{i,+}, \Delta v_{i+\HALF}) \\ \noalign{\medskip} \min(v_{i-1},v_i) \le v_{i,-} \le \max(v_{i-1},v_i) & \qquad \Longrightarrow\qquad & \delta v_{i,-} \to {\rm minmod}(\delta v_{i,-}, \Delta v_{i+\HALF}) \end{array}\right. \]

    where $\delta v_{i,\pm} = v_{i,\pm} - v_i$. This step is equivalent to Eq. [45] of Mignone (JCP, 2014). The previous operation replaces, starting with PLUTO 4.1, the conventional van-Leer limiting used in the original PPM formulation.
  • the parabola does not take any extrema inside the zone, see Eq. [46] in Mignone (JCP, 2014).

When interpolation is carried out in characteristic variables, we first compute vp and vm (unlimited) in primitive variables and then project the increments in characteristic space:

\[ \delta w_{i,\pm} = \vec{l}_i\cdot\delta v_{i,\pm} \]

At this point, Eq. [45] (written in terms of increments) is imposed on characteristic variables while Eq. [46] can be imposed in characteristic variables (PRIMITIVE_LIM == 0), primitive (PRIMITIVE_LIM == 1) or both (PRIMITIVE_LIM == 2).

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

References

  • "High-order conservative reconstruction schemes for finite volume methods in cylindrical and spherical coordinates" A. Mignone, JCP (2014), 270, 784.
  • "Accurate monotonicity- and extrema-preserving methods through adaptive nonlinear hybridizations" Rider, Greenough and Kamm, JCP (2007) 225, 1827
  • "A limiter for PPM that preserved accuracy at smooth extrema" Colella and Sekora, JCP (2008) 227, 7069

Definition in file ppm_states.c.

Macro Definition Documentation

#define PARABOLIC_LIM   1

Definition at line 257 of file ppm_states.c.

Function Documentation

void States ( const State_1D *  state,
int  beg,
int  end,
Grid *  grid 
)
Parameters
[in]statepointer to State_1D structure
[in]beginitial index of computation
[in]endfinal index of computation
[in]gridpointer to an array of Grid structures

Definition at line 68 of file ppm_states.c.

77 {
78  int i, nv;
79  double dv, **v;
80  double dvp, cp, *wp, *hp, **vp;
81  double dvm, cm, *wm, *hm, **vm;
82  PPM_Coeffs ppm_coeffs;
83  PLM_Coeffs plm_coeffs;
84 
85 /* ---------------------------------------------------------
86  0. Set pointers, compute geometrical coefficients
87  --------------------------------------------------------- */
88 
89  v = state->v;
90  vm = state->vm;
91  vp = state->vp;
92 
93  PPM_CoefficientsGet(&ppm_coeffs, g_dir);
94  #if SHOCK_FLATTENING == MULTID
95  PLM_CoefficientsGet(&plm_coeffs, g_dir);
96  #endif
97 
98  hp = ppm_coeffs.hp;
99  hm = ppm_coeffs.hm;
100 
101 /* ---------------------------------------
102  1b. Need to reconstruct 4-velocity ?
103  --------------------------------------- */
104 
105 #if RECONSTRUCT_4VEL
106  ConvertTo4vel (state->v, beg-2, end+2);
107 #endif
108 
109 /* ---------------------------------------------------------
110  2a. Define unlimited left and right interface values and
111  make sure they lie between adjacent cell averages.
112  --------------------------------------------------------- */
113 
114  #if PPM_ORDER == 3 || PPM_ORDER == 5
115  for (i = beg; i <= end; i++) {
116  wp = ppm_coeffs.wp[i];
117  wm = ppm_coeffs.wm[i];
118  VAR_LOOP(nv){
119  #if PPM_ORDER == 3
120  vp[i][nv] = wp[-1]*v[i-1][nv] + wp[0]*v[i][nv] + wp[1]*v[i+1][nv];
121  vm[i][nv] = wm[-1]*v[i-1][nv] + wm[0]*v[i][nv] + wm[1]*v[i+1][nv];
122  #elif PPM_ORDER == 5
123  vp[i][nv] = wp[-2]*v[i-2][nv] + wp[-1]*v[i-1][nv] +
124  wp[ 0]*v[i][nv] + wp[ 1]*v[i+1][nv] + wp[ 2]*v[i+2][nv];
125 
126  vm[i][nv] = wm[-2]*v[i-2][nv] + wm[-1]*v[i-1][nv] + wm[0]*v[i][nv]
127  + wm[ 1]*v[i+1][nv] + wm[ 2]*v[i+2][nv];
128  #endif
129  dvp = vp[i][nv] - v[i][nv];
130  dvm = vm[i][nv] - v[i][nv];
131 
132  dv = v[i+1][nv] - v[i][nv];
133  vp[i][nv] = v[i][nv] + MINMOD(dvp, dv);
134 
135  dv = v[i][nv] - v[i-1][nv];
136  vm[i][nv] = v[i][nv] + MINMOD(dvm, -dv);
137  }
138  }
139  #elif PPM_ORDER == 4 /* -- set a unique interface value -- */
140  for (i = beg-1; i <= end; i++) {
141  wp = ppm_coeffs.wp[i];
142  VAR_LOOP(nv){
143  vp[i][nv] = wp[-1]*v[i-1][nv] + wp[0]*v[i][nv]
144  + wp[ 1]*v[i+1][nv] + wp[2]*v[i+2][nv];
145 
146  dv = v[i+1][nv] - v[i][nv];
147  dvp = vp[i][nv] - v[i][nv];
148  vp[i][nv] = v[i][nv] + MINMOD(dvp, dv);
149  }
150  }
151  for (i = beg; i <= end; i++) VAR_LOOP(nv) vm[i][nv] = vp[i-1][nv];
152  #endif
153 
154 /* ---------------------------------------------------------
155  2b. Apply parabolic limiter: no new extrema should appear
156  in the parabola defined by vp, vm and v.
157  --------------------------------------------------------- */
158 
159  for (i = beg; i <= end; i++) {
160 #if SHOCK_FLATTENING == MULTID
161  if (state->flag[i] & FLAG_MINMOD) {
162  wp = plm_coeffs.wp;
163  wm = plm_coeffs.wm;
164  VAR_LOOP(nv) {
165  dvp = (v[i+1][nv] - v[i][nv])*wp[i];
166  dvm = (v[i][nv] - v[i-1][nv])*wm[i];
167  dv = MINMOD(dvp, dvm);
168  vp[i][nv] = v[i][nv] + dv*plm_coeffs.dp[i];
169  vm[i][nv] = v[i][nv] - dv*plm_coeffs.dm[i];
170  }
171  #if PHYSICS == RHD || PHYSICS == RMHD
172  VelocityLimiter (v[i], vp[i], vm[i]);
173  #endif
174  continue;
175  }
176 #endif
177 
178  #if PPM_ORDER == 0
179  cm = cp = 2.0;
180  #else
181  cm = (hm[i] + 1.0)/(hp[i] - 1.0);
182  cp = (hp[i] + 1.0)/(hm[i] - 1.0);
183  #endif
184 
185  for (nv = 0; nv < NVAR; nv++){
186  dvp = vp[i][nv] - v[i][nv];
187  dvm = vm[i][nv] - v[i][nv];
188 
189  if (dvp*dvm >= 0.0) dvp = dvm = 0.0;
190  else{
191  if (fabs(dvp) >= cm*fabs(dvm)) dvp = -cm*dvm;
192  else if (fabs(dvm) >= cp*fabs(dvp)) dvm = -cp*dvp;
193  }
194  vp[i][nv] = v[i][nv] + dvp;
195  vm[i][nv] = v[i][nv] + dvm;
196  }
197  #if PHYSICS == RHD || PHYSICS == RMHD
198  VelocityLimiter (v[i], vp[i], vm[i]);
199  #endif
200  }
201 
202 /* --------------------------------------------------------
203  3a. Apply 1D shock flattening
204  -------------------------------------------------------- */
205 
206  #if SHOCK_FLATTENING == YES
207  Flatten (state, beg, end, grid);
208  #endif
209 
210 /* -------------------------------------------
211  4. Assign face-centered magnetic field
212  ------------------------------------------- */
213 
214  #ifdef STAGGERED_MHD
215  for (i = beg-1; i <= end; i++) {
216  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
217  }
218  #endif
219 
220 /* --------------------------------------------------------
221  5. Evolve L/R states and center value by dt/2
222  -------------------------------------------------------- */
223 
224 #if TIME_STEPPING == CHARACTERISTIC_TRACING
225  CharTracingStep (state, beg, end, grid);
226 #endif
227 
228 /* --------------------------------------------------------
229  6. Convert back to 3-velocity
230  -------------------------------------------------------- */
231 
232 #if RECONSTRUCT_4VEL
233  ConvertTo3vel (state->v, beg-2, end+2);
234  ConvertTo3vel (state->vp, beg, end);
235  ConvertTo3vel (state->vm, beg, end);
236 #endif
237 
238 /* ----------------------------------------------
239  7. Obtain L/R states in conservative variables
240  ---------------------------------------------- */
241 
242  PrimToCons (state->vp, state->up, beg, end);
243  PrimToCons (state->vm, state->um, beg, end);
244 }
STATE_1D::v
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
VelocityLimiter
void VelocityLimiter(double *, double *, double *)
Definition: vel_limiter.c:16
ConvertTo4vel
void ConvertTo4vel(double **, int, int)
Definition: four_vel.c:4
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
Flatten
void Flatten(const State_1D *, int, int, Grid *)
Definition: flatten.c:4
CharTracingStep
void CharTracingStep(const State_1D *, int, int, Grid *)
Definition: char_tracing.c:198
PPM_CoefficientsGet
void PPM_CoefficientsGet(PPM_Coeffs *ppm_coeffs, int dir)
Definition: ppm_coeffs.c:591
PPM_COEFFS
Definition: ppm_coeffs.h:29
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
FLAG_MINMOD
#define FLAG_MINMOD
Reconstruct using MINMOD limiter.
Definition: pluto.h:179
BXn
int BXn
Definition: globals.h:75
STATE_1D::vp
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
Definition: structs.h:142
PPM_COEFFS::hm
double * hm
Definition: ppm_coeffs.h:33
PPM_COEFFS::wp
double ** wp
Definition: ppm_coeffs.h:30
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
GRID::beg
int beg
Global start index for the local array.
Definition: structs.h:115
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
ConvertTo3vel
void ConvertTo3vel(double **, int, int)
Definition: four_vel.c:25
MINMOD
#define MINMOD(a, b)
Definition: macros.h:112
STATE_1D::vm
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
PLM_COEFFS
Definition: plm_coeffs.h:38
PLM_CoefficientsGet
void PLM_CoefficientsGet(PLM_Coeffs *plm_coeffs, int dir)
Definition: plm_coeffs.c:79
PLM_COEFFS::dm
double * dm
Definition: plm_coeffs.h:44
PLM_COEFFS::dp
double * dp
Definition: plm_coeffs.h:43
PrimToCons
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
Definition: mappers.c:26
i
int i
Definition: analysis.c:2
PPM_COEFFS::hp
double * hp
Definition: ppm_coeffs.h:32
STATE_1D::bn
double * bn
Face magentic field, bn = bx(i+1/2)
Definition: structs.h:165
PPM_COEFFS::wm
double ** wm
Definition: ppm_coeffs.h:31
STATE_1D::um
double ** um
same as vm, in conservative vars
Definition: structs.h:146
PLM_COEFFS::wm
double * wm
Definition: plm_coeffs.h:42
STATE_1D::vL
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
STATE_1D::up
double ** up
same as vp, in conservative vars
Definition: structs.h:147
NVAR
#define NVAR
Definition: pluto.h:609
PLM_COEFFS::wp
double * wp
Definition: plm_coeffs.h:41

Here is the call graph for this function: