PLUTO
char_tracing.c
Go to the documentation of this file.
1 #include "pluto.h"
2 
3 #if RECONSTRUCTION == PARABOLIC || RECONSTRUCTION == WENO3
4 
5 /* ----------------------------------------------
6  Set REF_STATE to 1,2,3 to use cell centered
7  value (1), interpolated states (2) or
8  fastest wave (3, the usual PPM prescription)
9  ---------------------------------------------- */
10 
11 #if CHAR_LIMITING == YES
12  #define REF_STATE 3
13 #else
14  #define REF_STATE 3
15 #endif
16 
17 /* ****************************************************************** */
18 void CharTracingStep(const State_1D *state, int beg, int end, Grid *grid)
19 /*
20  *
21  *
22  ******************************************************************** */
23 {
24  int i, nv, k;
25  double dx, dtdx;
26  double dwh, d2w, dwp[NVAR], dwm[NVAR], dvp[NVAR], dvm[NVAR];
27  double nup=1.0, num=-1.0, nu[NVAR];
28  double Spp, Smm;
29  double *vc, *vp, *vm, **L, **R, *lambda;
30  static double **src;
31 
32  if (src == NULL){
33  src = ARRAY_2D(NMAX_POINT, NVAR, double);
34  }
35 
36  #if CHAR_LIMITING == NO
37  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
38  #endif
39 
40  PrimSource (state, beg, end, state->a2, state->h, src, grid);
41  for (i = beg; i <= end; i++){
42 
43  dx = grid[g_dir].dx[beg];
44  dtdx = g_dt/dx;
45 
46  vc = state->v[i];
47  vp = state->vp[i];
48  vm = state->vm[i];
49  L = state->Lp[i];
50  R = state->Rp[i];
51  lambda = state->lambda[i];
52 
53  /* --------------------------------------------------------------
54  1. Compute eigenvectors and eigenvalues if not yet defined
55  -------------------------------------------------------------- */
56 
57  #if CHAR_LIMITING == NO
58  PrimEigenvectors (vc, state->a2[i], state->h[i], lambda, L, R);
59  #if NVAR != NFLX
60  for (k = NFLX; k < NVAR; k++) lambda[k] = vc[V1];
61  #endif
62  #endif
63  for (k = 0; k < NVAR; k++) nu[k] = dtdx*lambda[k];
64 
65  /* --------------------------------------------------------------
66  2. Obtain characteristic variable increments dwp and dwm.
67  -------------------------------------------------------------- */
68 
69  for (nv = NVAR; nv--; ) {
70  dvp[nv] = vp[nv] - vc[nv];
71  dvm[nv] = vm[nv] - vc[nv];
72  }
73  PrimToChar(L, dvp, dwp);
74  PrimToChar(L, dvm, dwm);
75 
76  /* ----------------------------------------------------------------
77  3. Initialize vp and vm to the reference state.
78  Since this is somewhat arbitrary we set it using REF_STATE:
79 
80  REF_STATE==1: use cell-center value;
81  REF_STATE==2: interpolated value at base time level
82  REF_STATE==3: traditional PPM reference state (fastest
83  wave), minimize the size of the term
84  subject to characteristic limiting.
85 
86  Passive scalars use always REF_STATE == 2.
87  ---------------------------------------------------------------- */
88 
89  #if REF_STATE == 1
90  for (nv = NFLX; nv--; ) vp[nv] = vm[nv] = vc[nv];
91  #elif REF_STATE == 3
92  nup = MAX(nu[1], 0.0); num = MIN(nu[0], 0.0);
93  for (nv = NFLX; nv--; ){
94  dwh = vp[nv] - vm[nv];
95  d2w = vp[nv] + vm[nv] - 2.0*vc[nv];
96  vp[nv] -= 0.5*nup*(dwh + d2w*(3.0 - 2.0*nup));
97  vm[nv] -= 0.5*num*(dwh - d2w*(3.0 + 2.0*num));
98  }
99  #endif
100 
101  /* ------------------------------------------------------------------
102  4. Compute L/R states in primitive variables:
103 
104  - evolve characteristic variable increments by dt/2;
105  - discard contributions from waves not reaching the interface;
106  - project characteristic differences dwp and dwm onto right
107  eigenvectors
108  ------------------------------------------------------------------ */
109 
110  for (k = 0; k < NFLX; k++){
111  dwh = dwp[k] - dwm[k];
112  d2w = dwp[k] + dwm[k];
113  if (nu[k] >= 0.0) {
114  #if REF_STATE == 1
115  Spp = dwp[k] - 0.5*nu[k]*(dwh + d2w*(3.0 - 2.0*nu[k]));
116  #elif REF_STATE == 2
117  Spp = - 0.5*nu[k]*(dwh + d2w*(3.0 - 2.0*nu[k]));
118  #elif REF_STATE == 3
119  Spp = -0.5*(nu[k]-nup)*(dwh + 3.0*d2w) + (nu[k]*nu[k] - nup*nup)*d2w;
120  #endif
121  for (nv = NFLX; nv--; ) vp[nv] += Spp*R[nv][k];
122  } else {
123  #if REF_STATE == 1
124  Smm = dwm[k] - 0.5*nu[k]*(dwh - d2w*(3.0 + 2.0*nu[k]));
125  #elif REF_STATE == 2
126  Smm = - 0.5*nu[k]*(dwh - d2w*(3.0 + 2.0*nu[k]));
127  #elif REF_STATE == 3
128  Smm = -0.5*(nu[k]-num)*(dwh - 3.0*d2w) + (nu[k]*nu[k] - num*num)*d2w;
129  #endif
130  for (nv = NFLX; nv--; ) vm[nv] += Smm*R[nv][k];
131  }
132  }
133 
134  /* -------------------------------------------------------
135  5. Add source term to L/R states
136  ------------------------------------------------------- */
137 
138  for (nv = NFLX; nv--; ){
139  dwh = 0.5*g_dt*src[i][nv];
140  vp[nv] += dwh;
141  vm[nv] += dwh;
142  }
143 
144  /* -------------------------------------------------------
145  6. Repeat construction for passive scalars
146  ------------------------------------------------------- */
147 
148  #if NVAR != NFLX
149  if (nu[NFLX] >= 0.0) { /* -- scalars all move at the flow speed -- */
150  for (k = NFLX; k < NVAR; k++){
151  dwh = dwp[k] - dwm[k];
152  d2w = dwp[k] + dwm[k];
153  vp[k] -= 0.5*nu[k]*(dwh + d2w*(3.0 - 2.0*nu[k]));
154  }
155  }else{
156  for (k = NFLX; k < NVAR; k++){
157  dwh = dwp[k] - dwm[k];
158  d2w = dwp[k] + dwm[k];
159  vm[k] -= 0.5*nu[k]*(dwh - d2w*(3.0 + 2.0*nu[k]));
160  }
161  }
162  #endif
163  }
164 
165  #ifdef STAGGERED_MHD
166  for (i = beg-1; i <= end; i++) {
167  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
168  }
169  #endif
170 
171  CheckPrimStates (state->vm, state->vp, state->v, beg, end);
172 
173  for (i = beg; i <= end; i++) {
174  vp = state->vp[i];
175  vm = state->vm[i];
176  for (nv = NVAR; nv--; ) {
177  state->vh[i][nv] = 0.5*(vp[nv] + vm[nv]);
178  }
179  }
180 }
181 #undef REF_STATE
182 #endif /* INTERPOLATION == PARABOLIC || INTERPOLATION == WENO3 */
183 
184 
185 
186 
187 
188 #if RECONSTRUCTION == LINEAR
189 /* ------------------------------------------------------------
190  Set REF_STATE to 1,2,3 to use cell centered value (1),
191  interpolated states (2) or fastest wave (3, the standard
192  PPM prescription). Default is 1.
193  ------------------------------------------------------------ */
194 
195 #define REF_STATE 3
196 
197 /* *********************************************************************** */
198 void CharTracingStep(const State_1D *state, int beg, int end, Grid *grid)
199 /*
200  *
201  * PURPOSE:
202  *
203  * Compute 1D left and right interface states using piecewise
204  * linear reconstruction and the characteristic decomposition of the
205  * quasi-linear form of the equations.
206  *
207  * This is done by first extrapolating the cell center value to the
208  * interface using piecewise limited linear reconstruction
209  * on the characteristic variables.
210  *
211  * Left and right states are then evolved for the half time step
212  * using characteristic tracing if necessary.
213  *
214  * LAST MODIFIED
215  *
216  * May 22, 2012 by A. Mignone (mignone@ph.unito.it)
217  *
218  ************************************************************************* */
219 {
220  int i, j, k, nv;
221  double dx, dtdx, scrh;
222  double dv[NVAR], dw[NVAR];
223  double nu_max=1.0, nu_min=-1.0, nu[NVAR];
224  double *vc, *vp, *vm, **L, **R, *lambda;
225  double *dfR, *dfL;
226  #if GEOMETRY != CARTESIAN
227  double betaL[NVAR], betaR[NVAR];
228  #endif
229  static double **src;
230 
231 /* --------------------------------------------
232  allocate memory and set pointer shortcuts
233  -------------------------------------------- */
234 
235  if (src == NULL){
236  src = ARRAY_2D(NMAX_POINT, NVAR, double);
237  }
238 
239  #if GEOMETRY != CARTESIAN
240  dfL = grid[g_dir].dfL;
241  dfR = grid[g_dir].dfR;
242  #endif
243 
244 /* ---------------------------------------------
245  define some useful quantities, compute
246  source term and undivided differences
247  --------------------------------------------- */
248 
249  #if CHAR_LIMITING == NO
250  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
251  #endif
252 
253  PrimSource (state, beg, end, state->a2, state->h, src, grid);
254  for (i = beg; i <= end; i++){
255 
256  dx = grid[g_dir].dx[i];
257  dtdx = g_dt/dx;
258 
259  vc = state->v[i];
260  vp = state->vp[i];
261  vm = state->vm[i];
262  L = state->Lp[i];
263  R = state->Rp[i];
264  lambda = state->lambda[i];
265 
266  /* --------------------------------------------------------------
267  1. Compute eigenvectors and eigenvalues if not yet defined
268  -------------------------------------------------------------- */
269 
270  #if CHAR_LIMITING == NO
271  PrimEigenvectors (vc, state->a2[i], state->h[i], lambda, L, R);
272  #if NVAR != NFLX
273  for (k = NFLX; k < NVAR; k++) lambda[k] = vc[V1];
274  #endif
275  #endif
276  for (k = 0; k < NVAR; k++) nu[k] = dtdx*lambda[k];
277  nu_max = MAX(nu[1], 0.0); nu_min = MIN(nu[0], 0.0);
278 
279  /* ------------------------------------------------------------
280  1a. Geometry
281  ------------------------------------------------------------ */
282 
283 
284  #if GEOMETRY == CYLINDRICAL || GEOMETRY == POLAR
285  if (g_dir == IDIR) {
286  double *xR = grid[IDIR].xr;
287  for (k = NVAR; k--; ){
288 betaR[k] = betaL[k] = 0.0;
289 /*
290  scrh = nu[k]*dx;
291  betaR[k] = (nu[k] > 1.e-12 ? scrh/(6.0*xR[i] - 3.0*scrh):0.0);
292  betaL[k] = (nu[k] < -1.e-12 ? -scrh/(6.0*xR[i-1] - 3.0*scrh):0.0);
293 */
294  }
295  nu_max *= (1.0 - betaR[1]);
296  nu_min *= (1.0 + betaL[0]);
297  }else{
298  for (k = NVAR; k--; ) betaR[k] = betaL[k] = 0.0;
299  }
300  #endif
301 
302 
303  /* --------------------------------------------------------------
304  2. Obtain characteristic variable increment dw
305  -------------------------------------------------------------- */
306 
307  for (nv = NVAR; nv--; ) dv[nv] = vp[nv] - vm[nv];
308  PrimToChar(L, dv, dw);
309 
310  /* ----------------------------------------------------------------
311  3. Initialize vp and vm to the reference state.
312  Since this is somewhat arbitrary we set it using REF_STATE:
313 
314  REF_STATE==1: use cell-center value;
315  REF_STATE==2: interpolated value at base time level
316  REF_STATE==3: traditional PPM reference state (fastest
317  wave), minimize the size of the term
318  subject to characteristic limiting.
319 
320  Passive scalars use always REF_STATE == 2.
321  ---------------------------------------------------------------- */
322 
323  #if REF_STATE == 1
324  for (nv = NFLX; nv--; ) vp[nv] = vm[nv] = vc[nv];
325  #elif REF_STATE == 3
326  for (nv = NFLX; nv--; ) {
327  #if GEOMETRY == CARTESIAN
328  vp[nv] = vc[nv] + 0.5*dv[nv]*(1.0 - nu_max);
329  vm[nv] = vc[nv] - 0.5*dv[nv]*(1.0 + nu_min);
330  #else
331  vp[nv] = vc[nv] + dv[nv]*(dfR[i] - 0.5*nu_max);
332  vm[nv] = vc[nv] + dv[nv]*(dfL[i] - 0.5*nu_min);
333  #endif
334  }
335  #endif
336 
337  /* --------------------------------------------------------------------
338  4. Compute L/R states in primitive variables:
339 
340  Evolve interface states for dt/2 using characteristic tracing.
341  Depending on the choice of the reference state (i.e.
342  REF_STATE = 1,2,3) we compute, e.g. Vp(t+dt/2), as
343 
344  Vp(t+dt/2) = \sum [w + dw*dfR - nu*dw*(1-beta)]*R =
345 
346  = V + \sum [dfR - nu*(1-beta)]*dw*R (REF_STATE = 1)
347  = Vp + \sum [ - nu*(1-beta)]*dw*R (REF_STATE = 2)
348  = Vmax + \sum [nu_max*(1-beta_max) - nu*(1-beta)]*dw*R (REF_STATE = 3)
349 
350  where Vmax = Vp - nu_max*(1-beta_max) is the reference state
351  of the original PPM algorithm.
352  -------------------------------------------------------------------- */
353 
354  #ifdef GLM_MHD /* -- hot fix for backward test compatibility
355  must review this at some point !! -- */
356  nu_max = nu[1]; nu_min = nu[0];
357  #endif
358 
359  for (k = 0; k < NFLX; k++){
360  if (nu[k] >= 0.0) {
361  #if GEOMETRY != CARTESIAN
362  nu[k] *= (1.0 - betaR[k]);
363  #endif
364  #if REF_STATE == 1
365  dw[k] *= 0.5*(1.0 - nu[k]);
366  #elif REF_STATE == 2
367  dw[k] *= -0.5*nu[k];
368  #elif REF_STATE == 3
369  dw[k] *= 0.5*(nu_max - nu[k]);
370  #endif
371  for (nv = 0; nv < NFLX; nv++) vp[nv] += dw[k]*R[nv][k];
372  }else{
373  #if GEOMETRY != CARTESIAN
374  nu[k] *= (1.0 + betaL[k]);
375  #endif
376  #if REF_STATE == 1
377  dw[k] *= -0.5*(1.0 + nu[k]);
378  #elif REF_STATE == 2
379  dw[k] *= -0.5*nu[k];
380  #elif REF_STATE == 3
381  dw[k] *= 0.5*(nu_min - nu[k]);
382  #endif
383  for (nv = 0; nv < NFLX; nv++) vm[nv] += dw[k]*R[nv][k];
384  }
385  }
386 
387  /* -------------------------------------------------------
388  5. Add source term to L/R states
389  ------------------------------------------------------- */
390 
391  for (nv = NFLX; nv--; ) {
392  vp[nv] += 0.5*g_dt*src[i][nv];
393  vm[nv] += 0.5*g_dt*src[i][nv];
394  }
395 
396  /* -------------------------------------------------------
397  6. Repeat construction for passive scalars
398  ------------------------------------------------------- */
399 
400  #if NVAR != NFLX
401  for (nv = NFLX; nv < NVAR; nv++) {
402  #if GEOMETRY == CARTESIAN
403  vp[nv] = vc[nv] + 0.5*dv[nv];
404  vm[nv] = vc[nv] - 0.5*dv[nv];
405  #else
406  vp[nv] = vc[nv] + dv[nv]*dfR[i];
407  vm[nv] = vc[nv] + dv[nv]*dfL[i];
408  #endif
409  }
410  if (nu[NFLX] >= 0.0){ /* -- scalars all move at the flow speed -- */
411  scrh = -0.5*nu[NFLX];
412  #if GEOMETRY != CARTESIAN
413  scrh *= (1.0 - betaR[k]);
414  #endif
415  for (k = NFLX; k < NVAR; k++) vp[k] += dw[k]*scrh;
416  }else {
417  scrh = -0.5*nu[NFLX];
418  #if GEOMETRY != CARTESIAN
419  scrh *= (1.0 + betaL[k]);
420  #endif
421  for (k = NFLX; k < NVAR; k++) vm[k] += dw[k]*scrh;
422  }
423  #endif
424 
425  } /* -- end main loop on grid points -- */
426 
427 /* -------------------------------------------
428  Shock flattening (only 1D)
429  ------------------------------------------- */
430 
431  #if SHOCK_FLATTENING == ONED
432  Flatten (state, beg, end, grid);
433  #endif
434 
435 /* -------------------------------------------
436  Assign face-centered magnetic field
437  ------------------------------------------- */
438 
439  #ifdef STAGGERED_MHD
440  for (i = beg-1; i <= end; i++) {
441  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
442  }
443  #endif
444 
445 /* --------------------------------------------------------
446  evolve cell-center values by dt/2
447  -------------------------------------------------------- */
448 
449  CheckPrimStates (state->vm, state->vp, state->v, beg, end);
450 
451  for (i = beg; i <= end; i++){
452  vp = state->vp[i];
453  vm = state->vm[i];
454  for (nv = NVAR; nv--; ) {
455  state->vh[i][nv] = 0.5*(vp[nv] + vm[nv]);
456  }
457  }
458 
459 }
460 #undef REF_STATE
461 #endif /* INTERPOLATION == LINEAR */
MAX
#define MAX(a, b)
Definition: macros.h:101
STATE_1D::v
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
STATE_1D::vh
double ** vh
Primitive state at n+1/2 (only for one step method)
Definition: structs.h:162
Flatten
void Flatten(const State_1D *, int, int, Grid *)
Definition: flatten.c:4
STATE_1D::Lp
double *** Lp
Definition: structs.h:155
PrimEigenvectors
void PrimEigenvectors(double *q, double cs2, double h, double *lambda, double **LL, double **RR)
Definition: eigenv.c:91
GRID::xr
double * xr
Definition: structs.h:81
configure.scrh
tuple scrh
Definition: configure.py:200
CheckPrimStates
void CheckPrimStates(double **vM, double **vP, double **v0, int beg, int end)
Definition: check_states.c:4
g_dt
double g_dt
The current integration time step.
Definition: globals.h:118
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
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
MIN
#define MIN(a, b)
Definition: macros.h:104
CELL_CENTER
#define CELL_CENTER
Definition: pluto.h:205
CharTracingStep
void CharTracingStep(const State_1D *state, int beg, int end, Grid *grid)
Definition: char_tracing.c:46
NFLX
#define NFLX
Definition: mod_defs.h:32
IDIR
#define IDIR
Definition: pluto.h:193
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
GRID
Definition: structs.h:78
j
int j
Definition: analysis.c:2
STATE_1D::lambda
double ** lambda
Characteristic speed associated to Lp and Rp.
Definition: structs.h:156
k
int k
Definition: analysis.c:2
PrimSource
void PrimSource(const State_1D *, int, int, double *, double *, double **, Grid *)
Definition: prim_eqn.c:71
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
pluto.h
PLUTO main header file.
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
STATE_1D
Definition: structs.h:133
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::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
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
STATE_1D::Rp
double *** Rp
Left and right primitive eigenvectors.
Definition: structs.h:155
NVAR
#define NVAR
Definition: pluto.h:609