PLUTO
Functions
limo3_states.c File Reference

LimO3 reconstruction. More...

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

Go to the source code of this file.

Functions

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

Detailed Description

LimO3 reconstruction.

Provide a three-point stencil, third-order reconstruction algorithm based on the limiter function of Cada & Torrilhon

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

References

  • "Compact third-order limiter functions for finite volume methods", Cada & Torrilhon, JCP (2009) 228, 4118.

Definition in file limo3_states.c.

Function Documentation

double LimO3Func ( double  dvp,
double  dvm,
double  dx 
)
static

Implement the 3rd-order limiter function, Eq. [3.22]

Definition at line 246 of file limo3_states.c.

253 {
254  double r = 0.1;
255  double a,b,c,q, th, lim;
256  double eta, psi, eps = 1.e-12;
257 
258  th = dvm/(dvp + 1.e-16);
259 
260  q = (2.0 + th)/3.0;
261 
262  a = MIN(1.5,2.0*th);
263  a = MIN(q,a);
264  b = MAX(-0.5*th,a);
265  c = MIN(q,b);
266  psi = MAX(0.0,c);
267 
268  eta = r*dx;
269  eta = (dvm*dvm + dvp*dvp)/(eta*eta);
270  if ( eta <= 1.0 - eps) {
271  lim = q;
272  }else if (eta >= 1.0 + eps){
273  lim = psi;
274  }else{
275  psi = (1.0 - (eta - 1.0)/eps)*q
276  + (1.0 + (eta - 1.0)/eps)*psi;
277  lim = 0.5*psi;
278  }
279  return (lim);
280 }
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
eta
static double *** eta[3]
Definition: res_functions.c:94
GRID::dx
double * dx
Definition: structs.h:83
MIN
#define MIN(a, b)
Definition: macros.h:104
eps
#define eps
menu.c
tuple c
Definition: menu.py:375
q
static Runtime q
Definition: runtime_setup.c:490

Here is the caller graph for this function:

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

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 21 of file limo3_states.c.

26 {
27  int k, nv, i;
28  double dmm;
29  double **v, **vp, **vm, *dvp, *dvm, *dx;
30  double **L, **R, *lambda;
31  double dwp[NVAR], dwp_lim[NVAR];
32  double dwm[NVAR], dwm_lim[NVAR];
33  double dvpR, dvmR;
34  static double **dv;
35 
36 /* ----------------------------------------------------
37  0. Allocate memory, set pointer shortcuts
38  ---------------------------------------------------- */
39 
40  if (dv == NULL) {
41  dv = ARRAY_2D(NMAX_POINT, NVAR, double);
42  }
43 
44  v = state->v;
45  vp = state->vp;
46  vm = state->vm;
47 
48 #if RECONSTRUCT_4VEL
49  ConvertTo4vel (state->v, beg-1, end+1);
50 #endif
51 
52  dx = grid[g_dir].dx;
53 
54 /* ----------------------------------------------------
55  1. compute slopes and left and right interface values
56  ---------------------------------------------------- */
57 
58  #if CHAR_LIMITING == NO /* ----------------------------------------
59  Limiter on primitive variables
60  ---------------------------------------- */
61  for (i = beg-1; i <= end; i++){
62  for (nv = NVAR; nv--; ){
63  dv[i][nv] = v[i+1][nv] - v[i][nv];
64  }}
65 
66  for (i = beg; i <= end; i++){
67  dvp = dv[i];
68  dvm = dv[i-1];
69 
70  #if SHOCK_FLATTENING == MULTID
71  if (state->flag[i] & FLAG_MINMOD){
72  for (nv = NVAR; nv--; ) {
73  dmm = MINMOD(dvp[nv], dvm[nv]);
74  vp[i][nv] = v[i][nv] + 0.5*dmm;
75  vm[i][nv] = v[i][nv] - 0.5*dmm;
76  }
77  continue;
78  }
79  #endif
80 
81  for (nv = NVAR; nv--; ){
82  vp[i][nv] = v[i][nv] + 0.5*dvp[nv]*LimO3Func(dvp[nv], dvm[nv], dx[i]);
83  vm[i][nv] = v[i][nv] - 0.5*dvm[nv]*LimO3Func(dvm[nv], dvp[nv], dx[i]);
84  }
85  }
86 
87  #else /* --------------------------------------------
88  Limiter on characteristic variables
89  -------------------------------------------- */
90 
91  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
92  i = beg-1;
93  NVAR_LOOP(nv) dv[i][nv] = v[i+1][nv] - v[i][nv];
94 
95  for (i = beg; i <= end; i++){
96 
97  NVAR_LOOP(nv) dv[i][nv] = v[i+1][nv] - v[i][nv];
98  L = state->Lp[i];
99  R = state->Rp[i];
100  lambda = state->lambda[i];
101  dvp = dv[i];
102  dvm = dv[i-1];
103 
104  /* -------------------------------
105  project undivided differences
106  onto characteristic space
107  ------------------------------- */
108 
109  PrimEigenvectors (v[i], state->a2[i], state->h[i], lambda, L, R);
110  PrimToChar(L, dvp, dwp);
111  PrimToChar(L, dvm, dwm);
112 
113  #if SHOCK_FLATTENING == MULTID
114  if (state->flag[i] & FLAG_MINMOD){
115  for (nv = NVAR; nv--; ) {
116  dmm = MINMOD(dvp[nv], dvm[nv]);
117  vp[i][nv] = v[i][nv] + 0.5*dmm;
118  vm[i][nv] = v[i][nv] - 0.5*dmm;
119  }
120  continue;
121  }
122  #endif
123 
124  /* -----------------------------
125  limit undivided differences
126  ----------------------------- */
127 
128  for (k = NFLX; k--; ){
129  dwp_lim[k] = dwp[k]*LimO3Func(dwp[k], dwm[k], dx[i]);
130  dwm_lim[k] = dwm[k]*LimO3Func(dwm[k], dwp[k], dx[i]);
131  }
132  for (nv = NFLX; nv--; ){
133  dvpR = dvmR = 0.0;
134  #ifdef STAGGERED_MHD
135  if (nv == BXn) continue;
136  #endif
137  for (k = NFLX; k--; ){
138  dvpR += dwp_lim[k]*R[nv][k];
139  dvmR += dwm_lim[k]*R[nv][k];
140  }
141  vp[i][nv] = v[i][nv] + 0.5*dvpR;
142  vm[i][nv] = v[i][nv] - 0.5*dvmR;
143  }
144 
145  /* --------------------------------------
146  Compute limited slopes for tracers
147  exploiting the simple characteristic
148  structure, L=R=diag(1).
149  -------------------------------------- */
150 
151  #if NFLX != NVAR
152  for (nv = NFLX; nv < NVAR; nv++ ){
153  dvpR = dvp[nv]*LimO3Func(dvp[nv], dvm[nv], dx[i]);
154  dvmR = dvm[nv]*LimO3Func(dvm[nv], dvp[nv], dx[i]);
155  vp[i][nv] = v[i][nv] + 0.5*dvpR;
156  vm[i][nv] = v[i][nv] - 0.5*dvmR;
157  }
158  #endif
159  }
160 
161  #endif /* CHAR_LIMITING == YES */
162 
163 /* --------------------------------------------------
164  2. Since the third-order limiter is not TVD, we
165  need to ensure positivity of density and pressure
166  -------------------------------------------------- */
167 
168  for (i = beg; i <= end; i++){
169  dvp = dv[i];
170  dvm = dv[i-1];
171  if (vp[i][RHO] < 0.0 || vm[i][RHO] < 0.0){
172  dmm = MINMOD(dvp[RHO], dvm[RHO]);
173  vp[i][RHO] = v[i][RHO] + 0.5*dmm;
174  vm[i][RHO] = v[i][RHO] - 0.5*dmm;
175  }
176 
177  #if HAVE_ENERGY
178  if (vp[i][PRS] < 0.0 || vm[i][PRS] < 0.0){
179  dmm = MINMOD(dvp[PRS], dvm[PRS]);
180  vp[i][PRS] = v[i][PRS] + 0.5*dmm;
181  vm[i][PRS] = v[i][PRS] - 0.5*dmm;
182  }
183  #endif
184  #if ENTROPY_SWITCH
185  if (vp[i][ENTR] < 0.0 || vm[i][ENTR] < 0.0){
186  dmm = MINMOD(dvp[ENTR], dvm[ENTR]);
187  vp[i][ENTR] = v[i][ENTR] + 0.5*dmm;
188  vm[i][ENTR] = v[i][ENTR] - 0.5*dmm;
189  }
190  #endif
191 
192  /* -- relativistic limiter --*/
193 
194  #if PHYSICS == RHD || PHYSICS == RMHD
195  VelocityLimiter (v[i], vp[i], vm[i]);
196  #endif
197  }
198 
199 /* -------------------------------------------
200  3. Shock flattening
201  ------------------------------------------- */
202 
203  #if SHOCK_FLATTENING == ONED
204  Flatten (state, beg, end, grid);
205  #endif
206 
207 /* -------------------------------------------
208  4. Assign face-centered magnetic field
209  ------------------------------------------- */
210 
211  #ifdef STAGGERED_MHD
212  for (i = beg - 1; i <= end; i++) {
213  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
214  }
215  #endif
216 
217 /* --------------------------------------------------------
218  5. Evolve L/R states and center value by dt/2
219  -------------------------------------------------------- */
220 
221 #if TIME_STEPPING == CHARACTERISTIC_TRACING
222  CharTracingStep(state, beg, end, grid);
223 #elif TIME_STEPPING == HANCOCK
224  HancockStep(state, beg, end, grid);
225 #endif
226 
227 /* --------------------------------------------------------
228  6. Convert back to 3-velocity
229  -------------------------------------------------------- */
230 
231 #if RECONSTRUCT_4VEL
232  ConvertTo3vel (state->v, beg-1, end+1);
233  ConvertTo3vel (state->vp, beg, end);
234  ConvertTo3vel (state->vm, beg, end);
235 #endif
236 
237 /* ----------------------------------------------
238  7. Obtain L/R states in conservative variables
239  ---------------------------------------------- */
240 
241  PrimToCons (state->vp, state->up, beg, end);
242  PrimToCons (state->vm, state->um, beg, end);
243 }
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
PrimEigenvectors
void PrimEigenvectors(double *q, double cs2, double h, double *lambda, double **LL, double **RR)
Definition: eigenv.c:91
RHO
#define RHO
Definition: mod_defs.h:19
CharTracingStep
void CharTracingStep(const State_1D *, int, int, Grid *)
Definition: char_tracing.c:198
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
GRID::dx
double * dx
Definition: structs.h:83
STATE_1D::h
double * h
Enthalpy.
Definition: structs.h:159
FLAG_MINMOD
#define FLAG_MINMOD
Reconstruct using MINMOD limiter.
Definition: pluto.h:179
SoundSpeed2
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
Definition: eos.c:16
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
CELL_CENTER
#define CELL_CENTER
Definition: pluto.h:205
NFLX
#define NFLX
Definition: mod_defs.h:32
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
NVAR_LOOP
#define NVAR_LOOP(n)
Definition: pluto.h:618
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
k
int k
Definition: analysis.c:2
ConvertTo3vel
void ConvertTo3vel(double **, int, int)
Definition: four_vel.c:25
MINMOD
#define MINMOD(a, b)
Definition: macros.h:112
PrimToChar
void PrimToChar(double **L, double *v, double *w)
Definition: eigenv.c:593
STATE_1D::vm
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
LimO3Func
static double LimO3Func(double, double, double)
Definition: limo3_states.c:246
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
PrimToCons
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
Definition: mappers.c:26
i
int i
Definition: analysis.c:2
STATE_1D::bn
double * bn
Face magentic field, bn = bx(i+1/2)
Definition: structs.h:165
STATE_1D::um
double ** um
same as vm, in conservative vars
Definition: structs.h:146
STATE_1D::vL
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
STATE_1D::a2
double * a2
Sound speed squared.
Definition: structs.h:158
STATE_1D::up
double ** up
same as vp, in conservative vars
Definition: structs.h:147
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
NVAR
#define NVAR
Definition: pluto.h:609
HancockStep
void HancockStep(const State_1D *, int, int, Grid *)
Definition: hancock.c:136

Here is the call graph for this function: