PLUTO
plm_states.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief Piecewise linear reconstruction.
5 
6  Compute interface states using piecewise linear reconstruction
7  inside each zone.
8  Reconstruction is performed in primitive variable when
9  <tt> CHAR_LIMITING == NO </tt>) or characteristic variables
10  when <tt> CHAR_LIMITING == YES </tt>.
11 
12  The convention used throughout is that \b vL and \b vR are left and
13  right states with respect to the \e interface, while \b vm and \b vp
14  refer to the cell \e center:
15  \verbatim
16  vL(i)-> <-vR(i)
17  |----------*----------|----------*----------|
18  <-vm(i) (i) vp(i)-> (i+1)
19  \endverbatim
20 
21  The default setting (LIMITER == DEFAULT) applies a different limiter
22  to each variable:
23  - MC for density
24  - VanLeer for velcoity and magnetic field
25  - MinMod for pressure.
26 
27  Otherwise the same limiter can be imposed to all variables from
28  definitions.h.
29 
30  The approach followed here is taken from Mignone (JCP, 2014) "High-order
31  conservative reconstruction schemes for finite volume methods in cylindrical
32  and spherical geometries" where interface states are constructed as
33  \f[
34  V_i^{\pm} = V_i \;\pm\; \overline{\Delta V}_i \; d^\pm_i\,,
35  \qquad
36  \overline{\Delta V}_i = {\rm Lim}\left(\Delta V_i^+,\,\Delta V_i^-,\,
37  c^+_i, c^-_i\right)\,,
38  \qquad
39  \Delta V^\pm_i = \pm\left(V_{i\pm1}-V_i\right)w^\pm_i
40  \f]
41  where
42  \f[
43  w^\pm_i = \left|\frac{\Delta\xi_i}{\bar{\xi}_{i\pm1} - \bar{\xi}_i}\right| \,,
44  \qquad
45  d^\pm_i = \left|\frac{\xi_{i\pm\HALF} - \bar{\xi}_i}{\Delta\xi_i}\right|\,,
46  \qquad
47  c^\pm_i = \left|\frac{\bar{\xi}_{i\pm1} - \bar{\xi}_i}
48  {\xi_{i\pm\HALF}- \bar{\xi}_i}\right|
49  \f]
50  are interpolation coefficients returned by the PLM_Coeffs structure
51  (see the ::PLM_CoefficientsGet() function).
52  The slope limiter \f$\overline{\Delta V}_i=\f$ <tt> dv\_lim[i][nv] </tt>
53  is a function of the forward and backward (undivided) derivatives
54  and geometrical coefficients in case of non-uniform or non-Cartesian grids.
55  Different limiters are implemented through small macros defined in
56  plm_coeffs.h.
57  Note that \f$ \bar{\xi} \f$ defines the centroid of volume which differs
58  from the cell center in cylindrical and spherical geometires.
59 
60  A stencil of 3 zones is required for all limiters except for
61  the FOURTH_ORDER_LIM which requires 5 zones.
62 
63  \author A. Mignone (mignone@ph.unito.it)
64  \date June 16, 2015
65 
66  \b References
67  - "High-order conservative reconstruction schemes for finite
68  volume methods in cylindrical and spherical coordinates"
69  A. Mignone, JCP (2014), 270, 784.
70 */
71 /* ///////////////////////////////////////////////////////////////////// */
72 #include "pluto.h"
73 
74 static void MonotonicityTest(double **, double **, double **, int, int);
75 
76 #if CHAR_LIMITING == NO
77 static void FourthOrderLinear(const State_1D *, int, int, Grid *);
78 
79 /* ********************************************************************* */
80 void States (const State_1D *state, int beg, int end, Grid *grid)
81 /*!
82  * Compute states using piecewise linear interpolation.
83  *
84  * \param [in] state pointer to a State_1D structure
85  * \param [in] beg starting point where vp and vm must be computed
86  * \param [in] end final point where vp and vm must be computed
87  * \param [in] grid pointer to array of Grid structures
88  *
89  * \return This function has no return value.
90  *
91  ************************************************************************ */
92 {
93  int nv, i;
94  double **v, **vp, **vm;
95  double dv_lim[NVAR], dvp[NVAR], dvm[NVAR];
96  double cp, cm, wp, wm, dp, dm;
97  PLM_Coeffs plm_coeffs;
98  static double **dv;
99 
100  #if LIMITER == FOURTH_ORDER_LIM
101  FourthOrderLinear(state, beg, end, grid);
102  return;
103  #endif
104 
105 /* -----------------------------------------------------------
106  0. Memory allocation and pointer shortcuts, geometrical
107  coefficients and conversion to 4vel (if required)
108  ----------------------------------------------------------- */
109 
110  if (dv == NULL) {
111  dv = ARRAY_2D(NMAX_POINT, NVAR, double);
112  }
113 
114  v = state->v;
115  vp = state->vp;
116  vm = state->vm;
117 
118 #if UNIFORM_CARTESIAN_GRID == NO
119  PLM_CoefficientsGet (&plm_coeffs, g_dir);
120 #endif
121 
122 #if RECONSTRUCT_4VEL
123  ConvertTo4vel (state->v, beg-1, end+1);
124 #endif
125 
126 /* -------------------------------------------
127  compute undivided differences
128  ------------------------------------------- */
129 
130  for (i = beg-1; i <= end; i++){
131  NVAR_LOOP(nv) dv[i][nv] = v[i+1][nv] - v[i][nv];
132  }
133 
134 /* -------------------------------------------
135  2. Main spatial loop
136  ------------------------------------------- */
137 
138  for (i = beg; i <= end; i++){
139 
140  /* ---------------------------------------------------------
141  2a. compute forward (dvp) and backward (dvm) derivatives
142  --------------------------------------------------------- */
143 
144  #if UNIFORM_CARTESIAN_GRID == YES
145  cp = cm = 2.0;
146  wp = wm = 1.0;
147  dp = dm = 0.5;
148  NVAR_LOOP(nv) {
149  dvp[nv] = dv[i][nv];
150  dvm[nv] = dv[i-1][nv];
151  }
152  #else
153  cp = plm_coeffs.cp[i]; cm = plm_coeffs.cm[i];
154  wp = plm_coeffs.wp[i]; wm = plm_coeffs.wm[i];
155  dp = plm_coeffs.dp[i]; dm = plm_coeffs.dm[i];
156  NVAR_LOOP(nv) {
157  dvp[nv] = dv[i][nv]*wp;
158  dvm[nv] = dv[i-1][nv]*wm;
159  }
160  #endif
161 
162  /* ---------------------------------------------------------------
163  2b. if shock-flattening is enabled, revert to minmod limiter
164  --------------------------------------------------------------- */
165 
166 #if SHOCK_FLATTENING == MULTID
167  if (state->flag[i] & FLAG_FLAT) {
168  NVAR_LOOP(nv)vp[i][nv] = vm[i][nv] = v[i][nv];
169  continue;
170  }else if (state->flag[i] & FLAG_MINMOD) {
171  NVAR_LOOP(nv){
172  SET_MM_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
173  vp[i][nv] = v[i][nv] + dv_lim[nv]*dp;
174  vm[i][nv] = v[i][nv] - dv_lim[nv]*dm;
175  }
176  #if (PHYSICS == RHD || PHYSICS == RMHD)
177  VelocityLimiter (v[i], vp[i], vm[i]);
178  #endif
179  continue;
180  }
181 #endif
182 
183  /* ---------------------------------------------------------
184  2c. DEFAULT limiting: combination of different limiters
185  (this has some hystorical reasons)
186  --------------------------------------------------------- */
187 
188  #if LIMITER == DEFAULT
189  SET_MC_LIMITER(dv_lim[RHO], dvp[RHO], dvm[RHO], cp, cm);
190  #if PHYSICS != ADVECTION
191  EXPAND(SET_VL_LIMITER(dv_lim[VX1], dvp[VX1], dvm[VX1], cp, cm); ,
192  SET_VL_LIMITER(dv_lim[VX2], dvp[VX2], dvm[VX2], cp, cm); ,
193  SET_VL_LIMITER(dv_lim[VX3], dvp[VX3], dvm[VX3], cp, cm);)
194  #endif
195 
196  #if PHYSICS == MHD || PHYSICS == RMHD
197  EXPAND(SET_VL_LIMITER(dv_lim[BX1], dvp[BX1], dvm[BX1], cp, cm); ,
198  SET_VL_LIMITER(dv_lim[BX2], dvp[BX2], dvm[BX2], cp, cm); ,
199  SET_VL_LIMITER(dv_lim[BX3], dvp[BX3], dvm[BX3], cp, cm);)
200  #ifdef GLM_MHD
201  SET_MC_LIMITER(dv_lim[PSI_GLM], dvp[PSI_GLM], dvm[PSI_GLM], cp, cm);
202  #endif
203  #endif
204 
205  #if HAVE_ENERGY
206  SET_MM_LIMITER(dv_lim[PRS], dvp[PRS], dvm[PRS], cp, cm);
207  #endif
208 
209  #if NFLX != NVAR /* -- scalars: MC lim -- */
210  for (nv = NFLX; nv < NVAR; nv++){
211  SET_MC_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
212  }
213  #endif
214  #endif /* LIMITER == DEFAULT */
215 
216  /* ----------------------------------------
217  2d. construct (+) and (-) states
218  ---------------------------------------- */
219 
220  for (nv = 0; nv < NVAR; nv++){
221  #if LIMITER != DEFAULT /* -- same limiter for all variables -- */
222  SET_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
223  #endif
224 
225  vp[i][nv] = v[i][nv] + dv_lim[nv]*dp;
226  vm[i][nv] = v[i][nv] - dv_lim[nv]*dm;
227  }
228 
229  /* --------------------------------------
230  2e. Relativistic Limiter
231  -------------------------------------- */
232 
233  #if (PHYSICS == RHD || PHYSICS == RMHD)
234  VelocityLimiter (v[i], vp[i], vm[i]);
235  #endif
236  } /* -- end loop on zones -- */
237 
238 /* ----------------------------------------------------
239  3a. Check monotonicity
240  ---------------------------------------------------- */
241 
242 #if CHECK_MONOTONICITY == YES
243  MonotonicityTest(state->v, state->vp, state->vm, beg, end);
244 #endif
245 
246 /* -------------------------------------------
247  3b. Shock flattening
248  ------------------------------------------- */
249 
250 #if SHOCK_FLATTENING == ONED
251  Flatten (state, beg, end, grid);
252 #endif
253 
254 /* -------------------------------------------
255  4. Assign face-centered magnetic field
256  ------------------------------------------- */
257 
258 #ifdef STAGGERED_MHD
259  for (i = beg - 1; i <= end; i++) {
260  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
261  }
262 #endif
263 
264 /* --------------------------------------------------------
265  5. Evolve L/R states and center value by dt/2
266  -------------------------------------------------------- */
267 
268 #if TIME_STEPPING == CHARACTERISTIC_TRACING
269  CharTracingStep(state, beg, end, grid);
270 #elif TIME_STEPPING == HANCOCK
271  HancockStep(state, beg, end, grid);
272  /* The conservative hancock scheme is used only for RMHD with linear
273  interpolation in primitive variables.
274  In this case, the conversion to 3-vel is already carried out in
275  the predictor step and we need to skip it here. */
276  #if PRIMITIVE_HANCOCK == NO
277  PrimToCons (state->vp, state->up, beg, end);
278  PrimToCons (state->vm, state->um, beg, end);
279  return;
280  #endif
281 #endif
282 
283 /* --------------------------------------------------------
284  6. Convert back to 3-velocity
285  -------------------------------------------------------- */
286 
287 #if RECONSTRUCT_4VEL
288  ConvertTo3vel (state->v, beg-1, end+1);
289  ConvertTo3vel (state->vp, beg, end);
290  ConvertTo3vel (state->vm, beg, end);
291 #endif
292 
293 /* ----------------------------------------------
294  7. Obtain L/R states in conservative variables
295  ---------------------------------------------- */
296 
297  PrimToCons (state->vp, state->up, beg, end);
298  PrimToCons (state->vm, state->um, beg, end);
299 }
300 /* ********************************************************************** */
301 void FourthOrderLinear(const State_1D *state, int beg, int end, Grid *grid)
302 /*
303  * PURPOSE
304  *
305  * Compute interface states using Colella's fourth-order slope limiter
306  *
307  * Ref: Miller, G.H and P. COlella,
308  * "A high-order Eulerian Godunov Method for
309  * Elastic-Plastic Flow in Solids", JCP 167,131-176 (2001)
310  *
311  * +
312  *
313  * Saltzman, J, " An Unsplit 3D Upwind Method for
314  * Hyperbolic Conservation Laws",
315  * JCP 115, 153-168 (1994)
316  *
317  *********************************************************************** */
318 {
319  int i, nv;
320  static double **s;
321  static double **dv, **dvf, **dvc, **dvlim;
322  double scrh, dvp, dvm, dvl;
323  double **v, **vp, **vm;
324 
325  v = state->v;
326  vp = state->vp;
327  vm = state->vm;
328 
329  if (s == NULL){
330  s = ARRAY_2D(NMAX_POINT, NVAR, double);
331  dv = ARRAY_2D(NMAX_POINT, NVAR, double);
332  dvf = ARRAY_2D(NMAX_POINT, NVAR, double);
333  dvc = ARRAY_2D(NMAX_POINT, NVAR, double);
334  dvlim = ARRAY_2D(NMAX_POINT, NVAR, double);
335  }
336 
337  #if TIME_STEPPING == HANCOCK && PHYSICS != RMHD
338  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
339  #endif
340 /*
341  #if GEOMETRY != CARTESIAN
342  print1 ("! FourthOrderLinear: only Cartesian geometry supported\n");
343  QUIT_PLUTO(1);
344  #endif
345 */
346 /* -----------------------------------------------------------
347  compute undivided differences
348  ----------------------------------------------------------- */
349 
350  for (i = beg-2; i <= end+1; i++){
351  for (nv = 0; nv < NVAR; nv++) dv[i][nv] = v[i+1][nv] - v[i][nv];
352  }
353 
354  for (i = beg - 1; i <= end + 1; i++){
355  for (nv = 0; nv < NVAR; nv++){
356  dvp = dv[i][nv]; dvm = dv[i-1][nv];
357  dvc[i][nv] = 0.5*(dvp + dvm);
358  s[i][nv] = (dvp > 0.0 ? 0.5:-0.5)
359  + (dvm > 0.0 ? 0.5:-0.5);
360  dvlim[i][nv] = 2.0*MIN(fabs(dvp), fabs(dvm));
361  dvf[i][nv] = MIN(fabs(dvc[i][nv]), dvlim[i][nv])*s[i][nv];
362  }}
363 
364  for (i = beg; i <= end; i++){
365  for (nv = 0; nv < NVAR; nv++){
366  if (dv[i][nv]*dv[i-1][nv] > 0.0) {
367  scrh = 4.0/3.0*dvc[i][nv] - (dvf[i+1][nv] + dvf[i-1][nv])/6.0;
368  dvlim[i][nv] = MIN(fabs(scrh), dvlim[i][nv])*s[i][nv];
369  }else{
370  dvlim[i][nv] = 0.0;
371  }
372  }
373 
374  for (nv = NVAR; nv--; ) {
375  vp[i][nv] = v[i][nv] + 0.5*dvlim[i][nv];
376  vm[i][nv] = v[i][nv] - 0.5*dvlim[i][nv];
377  }
378  #if (PHYSICS == RHD || PHYSICS == RMHD)
379  VelocityLimiter(v[i], vp[i], vm[i]);
380  #endif
381  }
382 
383 /* -------------------------------------------
384  Shock flattening
385  ------------------------------------------- */
386 
387  #if SHOCK_FLATTENING == MULTID && CHAR_LIMITING == NO
388  Flatten (state, beg, end, grid);
389  #endif
390 
391 /* -------------------------------------------
392  Shock flattening
393  ------------------------------------------- */
394 
395  #if SHOCK_FLATTENING == ONED
396  Flatten (state, beg, end, grid);
397  #endif
398 
399 /* -------------------------------------------
400  Assign face-centered magnetic field
401  ------------------------------------------- */
402 
403  #ifdef STAGGERED_MHD
404  for (i = beg - 1; i <= end; i++) {
405  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
406  }
407  #endif
408 
409 /* --------------------------------------------------------
410  evolve center values by dt/2
411  -------------------------------------------------------- */
412 
413 
414  #if TIME_STEPPING == HANCOCK
415  HancockStep(state, beg, end, grid);
416  #endif
417 
418 /* -------------------------------------------
419  compute states in conservative variables
420  ------------------------------------------- */
421 
422  PrimToCons (state->vp, state->up, beg, end);
423  PrimToCons (state->vm, state->um, beg, end);
424 }
425 
426 #endif
427 
428 #if CHAR_LIMITING == YES
429 /* *********************************************************************** */
430 void States (const State_1D *state, int beg, int end, Grid *grid)
431 /*!
432  * Compute 1D left and right interface states using piecewise
433  * linear reconstruction and the characteristic decomposition of the
434  * quasi-linear form of the equations.
435  *
436  * This is done by first extrapolating the cell center value to the
437  * interface using piecewise limited linear reconstruction
438  * on the characteristic variables.
439  *
440  * Left and right states are then evolved for the half time step
441  * using characteristic tracing if necessary.
442  *
443  ************************************************************************* */
444 {
445  int i, j, k, nv;
446  double dvp[NVAR], dvm[NVAR], dv_lim[NVAR], dvc[NVAR], d2v;
447  double dw_lim[NVAR], dwp[NVAR], dwm[NVAR];
448  double dp, dm, dc;
449  double **vp, **vm, **v;
450  double **L, **R, *lambda;
451  double cp, cm, wp, wm, cpk[NVAR], cmk[NVAR];
452  double kstp[NVAR];
453  PLM_Coeffs plm_coeffs;
454  static double **dv;
455 
456 /* ---------------------------------------------
457  0. Allocate memory and set pointer shortcuts
458  --------------------------------------------- */
459 
460  if (dv == NULL) {
461  dv = ARRAY_2D(NMAX_POINT, NVAR, double);
462  }
463 
464  v = state->v;
465  vp = state->vp;
466  vm = state->vm;
467 
468  #if UNIFORM_CARTESIAN_GRID == NO
469  PLM_CoefficientsGet(&plm_coeffs, g_dir);
470  #endif
471 
472 /* ---------------------------------------------
473  Compute sound speed and then convert
474  to 4-vel if necessary.
475  --------------------------------------------- */
476 
477  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
478 
479 #if RECONSTRUCT_4VEL
480  ConvertTo4vel (state->v, beg-1, end+1); /* From now up to the end of the
481  * function, v contains the 4-vel */
482 #endif
483 
484  for (i = beg-1; i <= end; i++){
485  NVAR_LOOP(nv) dv[i][nv] = v[i+1][nv] - v[i][nv];
486  }
487 
488 /* --------------------------------------------------------------
489  Set the amount of steepening for each characteristic family.
490  Default is 2, but nonlinear fields may be safely set to 1
491  for strongly nonlinear problems.
492  -------------------------------------------------------------- */
493 
494  for (k = NVAR; k--; ) kstp[k] = 2.0;
495  #if PHYSICS == HD || PHYSICS == RHD
496  kstp[0] = kstp[1] = 1.0;
497  #elif PHYSICS == MHD
498  kstp[KFASTP] = kstp[KFASTM] = 1.0;
499  #if COMPONENTS > 1
500  kstp[KSLOWP] = kstp[KSLOWM] = 1.0;
501  #endif
502  #endif
503 
504 /* --------------------------------------------------------------
505  2. Start main spatial loop
506  -------------------------------------------------------------- */
507 
508  for (i = beg; i <= end; i++){
509 
510  L = state->Lp[i];
511  R = state->Rp[i];
512  lambda = state->lambda[i];
513 
514  PrimEigenvectors(v[i], state->a2[i], state->h[i], lambda, L, R);
515 
516  /* ---------------------------------------------------------------
517  2a. Project forward, backward (and centered) undivided
518  differences of primitive variables along characteristics,
519  dw(k) = L(k).dv
520  --------------------------------------------------------------- */
521 
522  #if UNIFORM_CARTESIAN_GRID == YES
523  cp = cm = 2.0;
524  wp = wm = 1.0;
525  dp = dm = 0.5;
526  NVAR_LOOP(nv) {
527  dvp[nv] = dv[i][nv];
528  dvm[nv] = dv[i-1][nv];
529 
530  k = nv;
531  cpk[k] = cmk[k] = kstp[k];
532  }
533  #else
534  cp = plm_coeffs.cp[i]; cm = plm_coeffs.cm[i];
535  wp = plm_coeffs.wp[i]; wm = plm_coeffs.wm[i];
536  dp = plm_coeffs.dp[i]; dm = plm_coeffs.dm[i];
537  NVAR_LOOP(nv) {
538  dvp[nv] = dv[i][nv]*wp;
539  dvm[nv] = dv[i-1][nv]*wm;
540 
541  /* -- Map 1 < kstp < 2 ==> 1 < ck < c -- */
542 
543  k = nv;
544  cpk[k] = (2.0 - cp) + (cp - 1.0)*kstp[k]; /* used only for general */
545  cmk[k] = (2.0 - cm) + (cm - 1.0)*kstp[k]; /* minmod limiter */
546  }
547  #endif
548 
549  PrimToChar(L, dvm, dwm);
550  PrimToChar(L, dvp, dwp);
551 
552  /* ----------------------------------------------------------
553  2b. Apply slope limiter to characteristic differences for
554  nv < NFLX, dw = Limiter(dwp, dwm).
555  ------------------------------------------------------- */
556 
557  #if SHOCK_FLATTENING == MULTID
558  if (state->flag[i] & FLAG_FLAT) {
559  for (k = NFLX; k--; ) dw_lim[k] = 0.0;
560  } else if (state->flag[i] & FLAG_MINMOD) {
561  for (k = NFLX; k--; ){
562  SET_MM_LIMITER(dw_lim[k], dwp[k], dwm[k], cp, cm);
563  }
564  } else
565  #endif
566  for (k = NFLX; k--; ){
567  #if LIMITER == DEFAULT
568  SET_GM_LIMITER(dw_lim[k], dwp[k], dwm[k], cpk[k], cmk[k]);
569  #else
570  SET_LIMITER(dw_lim[k], dwp[k], dwm[k], cp, cm);
571  #endif
572  }
573 
574  /* ------------------------------------------------------------------
575  2c. Project limited slopes in characteristic variables on right
576  eigenvectors to obtain primitive slopes: dv = \sum dw.R
577 
578  Also, enforce monotonicity in primitive variables as well.
579  ------------------------------------------------------------------ */
580 
581  for (nv = NFLX; nv--; ){
582  dc = 0.0;
583  for (k = 0; k < NFLX; k++) dc += dw_lim[k]*R[nv][k];
584 
585  if (dvp[nv]*dvm[nv] > 0.0){
586  d2v = ABS_MIN(cp*dvp[nv], cm*dvm[nv]);
587  dv_lim[nv] = MINMOD(d2v, dc);
588  }else dv_lim[nv] = 0.0;
589  }
590 
591  /* -----------------------------------------------------------------
592  2d. Repeat construction for passive scalars (tracers).
593  For a passive scalar, the primitive variable is the same as
594  the characteristic one. We use the MC limiter always.
595  ----------------------------------------------------------------- */
596 
597  #if NFLX != NVAR
598  for (nv = NFLX; nv < NVAR; nv++ ){
599  #if LIMITER == DEFAULT
600  SET_MC_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
601  #else
602  SET_LIMITER(dv_lim[nv], dvp[nv], dvm[nv], cp, cm);
603  #endif
604  }
605  #endif
606 
607  /* --------------------------------------------------------------------
608  2e. Build L/R states at time level t^n
609  -------------------------------------------------------------------- */
610 
611  for (nv = NVAR; nv--; ) {
612  vp[i][nv] = v[i][nv] + dv_lim[nv]*dp;
613  vm[i][nv] = v[i][nv] - dv_lim[nv]*dm;
614  }
615 
616  #if (PHYSICS == RHD || PHYSICS == RMHD)
617  VelocityLimiter (v[i], vp[i], vm[i]);
618  #endif
619 
620  } /* -- end main loop on grid points -- */
621 
622 /* ----------------------------------------------------
623  3a. Check monotonicity
624  ---------------------------------------------------- */
625 
626  #if CHECK_MONOTONICITY == YES
627  MonotonicityTest(state->v, state->vp, state->vm, beg, end);
628  #endif
629 
630 /* -------------------------------------------
631  3b. Shock flattening (only 1D)
632  ------------------------------------------- */
633 
634  #if SHOCK_FLATTENING == ONED
635  Flatten (state, beg, end, grid);
636  #endif
637 
638 /* -------------------------------------------
639  4. Assign face-centered magnetic field
640  ------------------------------------------- */
641 
642  #ifdef STAGGERED_MHD
643  for (i = beg-1; i <= end; i++) {
644  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
645  }
646  #endif
647 
648 /* --------------------------------------------------------
649  5. Evolve L/R states and center value by dt/2
650  -------------------------------------------------------- */
651 
652  #if TIME_STEPPING == CHARACTERISTIC_TRACING
653  CharTracingStep(state, beg, end, grid);
654  #elif TIME_STEPPING == HANCOCK
655  HancockStep(state, beg, end, grid);
656  #endif
657 
658 /* --------------------------------------------------------
659  6. Convert back to 3-velocity
660  -------------------------------------------------------- */
661 
662 #if RECONSTRUCT_4VEL
663  ConvertTo3vel (state->v, beg-1, end+1);
664  ConvertTo3vel (state->vp, beg, end);
665  ConvertTo3vel (state->vm, beg, end);
666 #endif
667 
668 /* ----------------------------------------------
669  7. Obtain L/R states in conservative variables
670  ---------------------------------------------- */
671 
672  PrimToCons (state->vp, state->up, beg, end);
673  PrimToCons (state->vm, state->um, beg, end);
674 }
675 #endif /* CHAR_LIMITING == YES */
676 
677 
678 /* ********************************************************************* */
679 void MonotonicityTest(double **v, double **vp, double **vm, int beg, int end)
680 /*
681  *
682  *********************************************************************** */
683 {
684  int i, nv;
685  double vp_max, vp_min, vm_max, vm_min;
686 
687  for (i = beg; i <= end; i++){ VAR_LOOP(nv) {
688  vp_max = MAX(v[i][nv], v[i+1][nv]) + 1.e-9;
689  vp_min = MIN(v[i][nv], v[i+1][nv]) - 1.e-9;
690 
691  vm_max = MAX(v[i][nv], v[i-1][nv]) + 1.e-9;
692  vm_min = MIN(v[i][nv], v[i-1][nv]) - 1.e-9;
693 
694  if (vp[i][nv] > vp_max || vp[i][nv] < vp_min){
695  print ("! States: monotonicity violated at + interface, i = %d\n",i);
696  print ("vp = %12.6e; vmax = %12.6e, vmin = %12.6e\n",
697  vp[i][nv], vp_max, vp_min);
698  QUIT_PLUTO(1);
699  }
700 
701  if (vm[i][nv] > vm_max || vm[i][nv] < vm_min){
702  print ("! States: monotonicity violated at - interface, i = %d\n",i);
703  print ("vm = %12.6e; vmax = %12.6e, vmin = %12.6e\n",
704  vm[i][nv], vm_max, vm_min);
705  QUIT_PLUTO(1);
706  }
707  }}
708 }
709 
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
VelocityLimiter
void VelocityLimiter(double *, double *, double *)
Definition: vel_limiter.c:16
ConvertTo4vel
void ConvertTo4vel(double **, int, int)
Definition: four_vel.c:4
VX2
#define VX2
Definition: mod_defs.h:29
SET_VL_LIMITER
#define SET_VL_LIMITER(dv, dvp, dvm, cp, cm)
Definition: plm_coeffs.h:114
Flatten
void Flatten(const State_1D *, int, int, Grid *)
Definition: flatten.c:4
FLAG_FLAT
#define FLAG_FLAT
Reconstruct using FLAT limiter.
Definition: pluto.h:180
ABS_MIN
#define ABS_MIN(a, b)
Definition: macros.h:107
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
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
CharTracingStep
void CharTracingStep(const State_1D *, int, int, Grid *)
Definition: char_tracing.c:198
PLM_COEFFS::cm
double * cm
Definition: plm_coeffs.h:40
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
VX1
#define VX1
Definition: mod_defs.h:28
STATE_1D::h
double * h
Enthalpy.
Definition: structs.h:159
SET_GM_LIMITER
#define SET_GM_LIMITER(dv, dvp, dvm, cp, cm)
Definition: plm_coeffs.h:96
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
MIN
#define MIN(a, b)
Definition: macros.h:104
CELL_CENTER
#define CELL_CENTER
Definition: pluto.h:205
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
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
States
void States(const State_1D *state, int beg, int end, Grid *grid)
Definition: plm_states.c:80
NVAR_LOOP
#define NVAR_LOOP(n)
Definition: pluto.h:618
GRID
Definition: structs.h:78
SET_MC_LIMITER
#define SET_MC_LIMITER(dv, dvp, dvm, cp, cm)
Definition: plm_coeffs.h:119
PHYSICS
#define PHYSICS
Definition: definitions_01.h:1
j
int j
Definition: analysis.c:2
STATE_1D::lambda
double ** lambda
Characteristic speed associated to Lp and Rp.
Definition: structs.h:156
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
PLM_COEFFS::cp
double * cp
Definition: plm_coeffs.h:39
KFASTP
Definition: mod_defs.h:58
k
int k
Definition: analysis.c:2
ConvertTo3vel
void ConvertTo3vel(double **, int, int)
Definition: four_vel.c:25
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
KFASTM
Definition: mod_defs.h:58
MINMOD
#define MINMOD(a, b)
Definition: macros.h:112
PrimToChar
void PrimToChar(double **L, double *v, double *w)
Definition: eigenv.c:593
SET_MM_LIMITER
#define SET_MM_LIMITER(dv, dvp, dvm, cp, cm)
Definition: plm_coeffs.h:76
STATE_1D::vm
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
MHD
#define MHD
Definition: pluto.h:111
PLM_COEFFS
Definition: plm_coeffs.h:38
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
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
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
MonotonicityTest
static void MonotonicityTest(double **, double **, double **, int, int)
Definition: plm_states.c:679
BX1
#define BX1
Definition: mod_defs.h:25
STATE_1D::um
double ** um
same as vm, in conservative vars
Definition: structs.h:146
VX3
#define VX3
Definition: mod_defs.h:30
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::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
STATE_1D::Rp
double *** Rp
Left and right primitive eigenvectors.
Definition: structs.h:155
BX2
#define BX2
Definition: mod_defs.h:26
NVAR
#define NVAR
Definition: pluto.h:609
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
RMHD
#define RMHD
Definition: pluto.h:112
FourthOrderLinear
static void FourthOrderLinear(const State_1D *, int, int, Grid *)
Definition: plm_states.c:301
GLM_MHD
#define GLM_MHD
Definition: glm.h:25
PLM_COEFFS::wp
double * wp
Definition: plm_coeffs.h:41
HancockStep
void HancockStep(const State_1D *, int, int, Grid *)
Definition: hancock.c:136