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