PLUTO
Functions
weno3_states.c File Reference
#include "pluto.h"
Include dependency graph for weno3_states.c:

Go to the source code of this file.

Functions

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

Function Documentation

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.

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.

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 6 of file weno3_states.c.

21 {
22  int k, nv, i;
23  double dmm, dx2;
24  double **v, **vp, **vm, *dvp, *dvm, *dx;
25  double **Lv, **Rv, *lambda;
26  double *R, *L, *P, *M;
27  double b0, S0, wp, tau;
28  double b1, S1, wm;
29  double dwp[NVAR], dwp_lim[NVAR];
30  double dwm[NVAR], dwm_lim[NVAR];
31  double dvpR, dvmR;
32  static double **Rg, **Lg, **Pg, **Mg; /* -- interpolation coeffs -- */
33  static double **dv;
34 
35 /* -----------------------------------------------------
36  0. Allocate memory and set pointer shortcuts
37  ----------------------------------------------------- */
38 
39  if (dv == NULL) {
40  dv = ARRAY_2D(NMAX_POINT, NVAR, double);
41  Rg = ARRAY_2D(DIMENSIONS, NMAX_POINT, double);
42  Lg = ARRAY_2D(DIMENSIONS, NMAX_POINT, double);
43  Pg = ARRAY_2D(DIMENSIONS, NMAX_POINT, double);
44  Mg = ARRAY_2D(DIMENSIONS, NMAX_POINT, double);
45  WENO3_COEFF(Rg, Lg, Pg, Mg, grid);
46  }
47 
48  v = state->v;
49  vp = state->vp;
50  vm = state->vm;
51 
52 #if RECONSTRUCT_4VEL
53  ConvertTo4vel (state->v, beg-1, end+1);
54 #endif
55 
56  R = Rg[g_dir]; L = Lg[g_dir];
57  P = Pg[g_dir]; M = Mg[g_dir];
58  dx = grid[g_dir].dx;
59 
60 /* ----------------------------------------------------
61  1. compute slopes and left and right interface values
62  ---------------------------------------------------- */
63 
64  #if CHAR_LIMITING == NO /* ----------------------------------------
65  a: Limiter on primitive variables
66  ---------------------------------------- */
67  for (i = beg-1; i <= end; i++){
68  for (nv = NVAR; nv--; ){
69  dv[i][nv] = v[i+1][nv] - v[i][nv];
70  }}
71 
72  for (i = beg; i <= end; i++){
73  dvp = dv[i];
74  dvm = dv[i-1];
75 
76  #if SHOCK_FLATTENING == MULTID
77  if (state->flag[i] & FLAG_MINMOD){
78  for (nv = NVAR; nv--; ) {
79  dmm = MINMOD(dvp[nv], dvm[nv]);
80  vp[i][nv] = v[i][nv] + 0.5*dmm;
81  vm[i][nv] = v[i][nv] - 0.5*dmm;
82  }
83  continue;
84  }
85  #endif
86 
87  dx2 = dx[i]*dx[i];
88  for (nv = 0; nv < NVAR; nv++){
89  b0 = dvp[nv]*dvp[nv] + dx2;
90  b1 = dvm[nv]*dvm[nv] + dx2;
91 
92  tau = dvp[nv] - dvm[nv];
93  tau = tau*tau;
94 
95  S0 = 1.0 + tau/b0;
96  S1 = 1.0 + tau/b1;
97 /*
98 S0 = 1.0/(b0 + 1.e-6)/(b0 + 1.e-6);
99 S1 = 1.0/(b1 + 1.e-6)/(b1 + 1.e-6);
100  vp[nv] = v[i][nv] + (S0*dvp[nv] + 0.5*S1*dvm[nv])/(2.0*S0 + S1);
101  vm[nv] = v[i][nv] - (S1*dvm[nv] + 0.5*S0*dvp[nv])/(2.0*S1 + S0);
102 */
103 
104  vp[i][nv] = v[i][nv] + (S0*R[i]*dvp[nv] + P[i]*S1*R[i-1]*dvm[nv])
105  /(S0 + P[i]*S1);
106  vm[i][nv] = v[i][nv] - (M[i]*S0*L[i]*dvp[nv] + S1*L[i-1]*dvm[nv])
107  /(M[i]*S0 + S1);
108 
109  }
110  }
111 
112  #else /* --------------------------------------------
113  b: Limiter on characteristic variables
114  -------------------------------------------- */
115 
116  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
117  i = beg-1;
118  for (nv = NVAR; nv--; ) dv[i][nv] = v[i+1][nv] - v[i][nv];
119 
120  for (i = beg; i <= end; i++){
121 
122  NVAR_LOOP(nv) dv[i][nv] = v[i+1][nv] - v[i][nv];
123 
124  Lv = state->Lp[i];
125  Rv = state->Rp[i];
126  lambda = state->lambda[i];
127  dvp = dv[i];
128  dvm = dv[i-1];
129 
130  /* -------------------------------
131  project undivided differences
132  onto characteristic space
133  ------------------------------- */
134 
135  PrimEigenvectors (v[i], state->a2[i], state->h[i], lambda, Lv, Rv);
136  PrimToChar(Lv, dvp, dwp);
137  PrimToChar(Lv, dvm, dwm);
138 
139  #if SHOCK_FLATTENING == MULTID
140  if (state->flag[i] & FLAG_MINMOD){
141  for (nv = NVAR; nv--; ) {
142  dmm = MINMOD(dvp[nv], dvm[nv]);
143  vp[i][nv] = v[i][nv] + 0.5*dmm;
144  vm[i][nv] = v[i][nv] - 0.5*dmm;
145  }
146  continue;
147  }
148  #endif
149 
150  dx2 = dx[i]*dx[i];
151  for (k = NVAR; k--; ){
152  b0 = dwp[k]*dwp[k] + dx2;
153  b1 = dwm[k]*dwm[k] + dx2;
154 
155  tau = dwp[k] - dwm[k];
156  tau = tau*tau;
157 
158  S0 = 1.0 + tau/b0;
159  S1 = 1.0 + tau/b1;
160 
161  dwp_lim[k] = (S0*dwp[k] + 0.5*S1*dwm[k])/(2.0*S0 + S1);
162  dwm_lim[k] = (S1*dwm[k] + 0.5*S0*dwp[k])/(2.0*S1 + S0);
163 /*
164  dwp_lim[k] = (S0*R[i]*dwp[k] + P[i]*S1*R[i-1]*dwm[k])/(S0 + P[i]*S1);
165  dwm_lim[k] = (M[i]*S0*L[i]*dwp[nv] + S1*L[i-1]*dwm[nv])/(M[i]*S0 + S1);
166 */
167 
168  }
169  for (nv = NFLX; nv--; ){
170  dvpR = dvmR = 0.0;
171  #ifdef STAGGERED_MHD
172  if (nv == BXn) continue;
173  #endif
174  for (k = NFLX; k--; ){
175  dvpR += dwp_lim[k]*Rv[nv][k];
176  dvmR += dwm_lim[k]*Rv[nv][k];
177  }
178  vp[i][nv] = v[i][nv] + dvpR;
179  vm[i][nv] = v[i][nv] - dvmR;
180  }
181 
182  /* --------------------------------------------------
183  Passive scalars are treated in a separate loop
184  since the characteristic structure is simpler,
185  that is, L=R=diag(1) and the characteristic
186  variable coincides with the primitive one.
187  -------------------------------------------------- */
188 
189  #if NFLX != NVAR
190  for (nv = NFLX; nv < NVAR; nv++ ){
191  vp[i][nv] = v[i][nv] + dwp_lim[nv];
192  vm[i][nv] = v[i][nv] - dwm_lim[nv];
193  }
194  #endif
195 
196  }
197 
198  #endif /* CHAR_LIMITING == YES */
199 
200 /* --------------------------------------------------
201  2. Ensure positivity of density and pressure
202  -------------------------------------------------- */
203 
204  for (i = beg; i <= end; i++){
205  dvp = dv[i];
206  dvm = dv[i-1];
207  if (vp[i][RHO] < 0.0 || vm[i][RHO] < 0.0){
208  for (nv = NFLX; nv--; ){
209  dmm = MINMOD(dvp[RHO], dvm[RHO]);
210  vp[i][RHO] = v[i][RHO] + 0.5*dmm;
211  vm[i][RHO] = v[i][RHO] - 0.5*dmm;
212  }
213  }
214 
215  #if HAVE_ENERGY
216  if (vp[i][PRS] < 0.0 || vm[i][PRS] < 0.0){
217  for (nv = NFLX; nv--; ){
218  dmm = MINMOD(dvp[PRS], dvm[PRS]);
219  vp[i][PRS] = v[i][PRS] + 0.5*dmm;
220  vm[i][PRS] = v[i][PRS] - 0.5*dmm;
221  }
222  }
223  #endif
224  #if ENTROPY_SWITCH
225  if (vp[i][ENTR] < 0.0 || vm[i][ENTR] < 0.0){
226  dmm = MINMOD(dvp[ENTR], dvm[ENTR]);
227  vp[i][ENTR] = v[i][ENTR] + 0.5*dmm;
228  vm[i][ENTR] = v[i][ENTR] - 0.5*dmm;
229  }
230  #endif
231 
232  /* -- Relativistic Limiter -- */
233 
234  #if PHYSICS == RHD || PHYSICS == RMHD
235  VelocityLimiter (v[i], vp[i], vm[i]);
236  #endif
237 
238  }
239 
240 /* -------------------------------------------
241  3. Shock flattening
242  ------------------------------------------- */
243 
244 #if SHOCK_FLATTENING == ONED
245  Flatten (state, beg, end, grid);
246 #endif
247 
248 /* -------------------------------------------
249  4. Assign face-centered magnetic field
250  ------------------------------------------- */
251 
252  #ifdef STAGGERED_MHD
253  for (i = beg - 1; i <= end; i++) {
254  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
255  }
256  #endif
257 
258 /* --------------------------------------------------------
259  5. Evolve L/R states and center value by dt/2
260  -------------------------------------------------------- */
261 
262 #if TIME_STEPPING == CHARACTERISTIC_TRACING
263  CharTracingStep(state, beg, end, grid);
264 #elif TIME_STEPPING == HANCOCK
265  HancockStep(state, beg, end, grid);
266 #endif
267 
268 /* --------------------------------------------------------
269  6. Convert back to 3-velocity
270  -------------------------------------------------------- */
271 
272 #if RECONSTRUCT_4VEL
273  ConvertTo3vel (state->v, beg-1, end+1);
274  ConvertTo3vel (state->vp, beg, end);
275  ConvertTo3vel (state->vm, beg, end);
276 #endif
277 
278 /* ----------------------------------------------
279  7. Obtain L/R states in conservative variables
280  ---------------------------------------------- */
281 
282  PrimToCons (state->vp, state->up, beg, end);
283  PrimToCons (state->vm, state->um, beg, end);
284 }
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
WENO3_COEFF
static void WENO3_COEFF(double **, double **, double **, double **, Grid *)
Definition: weno3_states.c:287
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
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
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2
HancockStep
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 WENO3_COEFF ( double **  R,
double **  L,
double **  P,
double **  M,
Grid *  grid 
)
static

Definition at line 287 of file weno3_states.c.

294 {
295  int d, i, beg, end;
296  double *dV;
297 
298  for (d = 0; d < DIMENSIONS; d++){
299  dV = grid[d].dV;
300  beg = 0; end = grid[d].np_tot - 2;
301  for (i = beg; i <= end; i++){
302  R[d][i] = dV[i + 1]/(dV[i] + dV[i + 1]);
303  L[d][i] = 1.0 - R[d][i];
304  }
305  beg = 1; end = grid[d].np_tot - 2;
306  for (i = beg; i <= end; i++){
307  P[d][i] = dV[i + 1]/(dV[i] + dV[i - 1]);
308  M[d][i] = dV[i - 1]/(dV[i] + dV[i + 1]);
309  }
310  }
311 }
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
GRID::beg
int beg
Global start index for the local array.
Definition: structs.h:115
i
int i
Definition: analysis.c:2
GRID::np_tot
int np_tot
Total number of points in the local domain (boundaries included).
Definition: structs.h:100
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the caller graph for this function: