PLUTO
hancock.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief MUSCL-Hancock predictor step.
5 
6  Use time-extrapolation to compute interface states and cell-centered
7  value at the half-time step,
8  \f$ V_{i,\pm}^{n+\HALF} = V_{i,\pm}^n + \partial_t V \Delta t^n/2 \f$
9 
10  This is done using the one-dimensional primitive form ot the
11  equations for the HD, RHD and MHD modules since we have at disposal
12  \f$ \partial_t V = -A\partial V_x + S \f$.
13 
14  Conversely, for relativistic MHD, we adopt the one-dimensional
15  conservative form of the equations (conservative Hancock predictor) and
16  compute the primitive values using
17  \f$ \partial_t V \approx (V^{n+\HALF}_i-V^n_i)/(\Delta t/2) \f$.
18  This requires taking the following steps:
19 
20  - Convert to conservative: vp -> up, vm -> um
21  - Evolve u(n+1/2) = u(n) - dt/(2dx)*[F(vp) - F(vm)]
22  - Convert to primitive u(n+1/2) -> v(n+1/2)
23  - Add time incerement v(n+1/2)-v(n) to left/right states
24 
25  \author A. Mignone (mignone@ph.unito.it)
26  \date June 24, 2015
27 */
28 /* ///////////////////////////////////////////////////////////////////// */
29 #include "pluto.h"
30 
31 #if PHYSICS != RMHD
32 /* ********************************************************************* */
33 void HancockStep (const State_1D *state, int beg, int end, Grid *grid)
34 /*!
35  * Use Taylor expansion to compute time-centered states in primitive
36  * variables.
37  *
38  * \param [in,out] state pointer to a State_1D structure
39  * \param [in] beg initial index of computation
40  * \param [in] end final index of computation
41  * \param [in grid pointer to array of Grid structures
42  *
43  * \b References
44  * - "Riemann Solvers and Numerical Methods for Fluid Dynamics" \n
45  * E.F. Toro
46  *
47  *********************************************************************** */
48 {
49  int nv, i;
50  double scrh, dt_2, ch2;
51  double Adv[NVAR], dv[NVAR];
52  double *vp, *vm, *vc;
53  static double **src, *d_dl;
54 
55 /* -----------------------------------------
56  Check scheme compatibility
57  ----------------------------------------- */
58 
59  #if RECONSTRUCTION != LINEAR
60  print1 ("! MUSCL-Hancock scheme works with Linear reconstruction only\n");
61  QUIT_PLUTO(1);
62  #endif
63 
64 
65  if (src == NULL){
66  src = ARRAY_2D(NMAX_POINT, NVAR, double);
67  d_dl = ARRAY_1D(NMAX_POINT, double);
68  }
69 
70 /* -------------------------------------------------------
71  compute geometric factor 1/(h*dx)
72  ------------------------------------------------------- */
73 
74  #if GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL
75 
76  for (i = beg; i <= end; i++) {
77  d_dl[i] = 1.0/grid[g_dir].dx[i];
78  }
79 
80  #elif GEOMETRY == SPHERICAL || GEOMETRY == POLAR
81 
82  if (g_dir == IDIR){
83  for (i = beg; i <= end; i++) {
84  d_dl[i] = 1.0/grid[g_dir].dx[i];
85  }
86  }else if (g_dir == JDIR){
87  for (i = beg; i <= end; i++) {
88  d_dl[i] = 1.0/(grid[IDIR].x[g_i]*grid[JDIR].dx[i]);
89  }
90  }
91 
92  #endif
93 
94  #if CHAR_LIMITING == NO
95  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
96  #endif
97 
98  PrimSource (state, beg, end, state->a2, state->h, src, grid);
99  dt_2 = 0.5*g_dt;
100  for (i = beg; i <= end; i++) {
101  vc = state->v[i];
102  vm = state->vm[i];
103  vp = state->vp[i];
104 
105  for (nv = NVAR; nv--; ) dv[nv] = vp[nv] - vm[nv];
106 
107  PrimRHS (vc, dv, state->a2[i], state->h[i], Adv);
108  for (nv = NVAR; nv--; ) {
109  scrh = dt_2*(d_dl[i]*Adv[nv] - src[i][nv]);
110  vp[nv] -= scrh;
111  vm[nv] -= scrh;
112  }
113  }
114 
115  CheckPrimStates (state->vm, state->vp, state->v, beg, end);
116 
117 /* -----------------------------------
118  evolve center value by dt/2
119  ----------------------------------- */
120 
121  for (i = beg; i <= end; i++) {
122  vc = state->vh[i];
123  vp = state->vp[i];
124  vm = state->vm[i];
125  NVAR_LOOP(nv) vc[nv] = 0.5*(vp[nv] + vm[nv]);
126  }
127 #if RECONSTRUCT_4VEL
128  ConvertTo3vel (state->vh, beg, end);
129 #endif
130 
131 }
132 
133 #elif PHYSICS == RMHD
134 
135 /* ********************************************************************* */
136 void HancockStep (const State_1D *state, int beg, int end, Grid *grid)
137 /*
138  *
139  *********************************************************************** */
140 #define CHAR_UPWIND NO
141  {
142  int nv, i, ifail;
143  double dBx, dpsi, dtdx, dt, dtdV;
144  double *A, *dx, *dV, r;
145  double *vp, *vm, *vc, dvp[NVAR], dvm[NVAR];
146  double P0, P1, P2;
147  double **rhs;
148  static double **fp, **fm, **uh, **u;
149  static double *pp, *pm, *hp, *hm;
150  static double *lambda_max, *lambda_min;
151 
152 #if GEOMETRY != CARTESIAN && GEOMETRY != CYLINDRICAL
153  print1 ("! Hancock does not work in this geometry \n");
154  QUIT_PLUTO(1);
155 #endif
156 #if BODY_FORCE != NO
157  print ("! Conservative Hancock scheme does not support gravity\n");
158  QUIT_PLUTO(1);
159 #endif
160 #if RECONSTRUCTION != LINEAR
161  print1 ("! The MUSCL-Hancock scheme works with Linear reconstruction only\n");
162  QUIT_PLUTO(1);
163 #endif
164 #if (PARABOLIC_FLUX & EXPLICIT)
165  print1 ("! Conservative MUSCL-Hancock incompatible with Explicit Diffusion\n");
166  QUIT_PLUTO(1);
167 #endif
168 
169  if (fp == NULL){
170  fp = ARRAY_2D(NMAX_POINT, NVAR, double);
171  fm = ARRAY_2D(NMAX_POINT, NVAR, double);
172  pp = ARRAY_1D(NMAX_POINT, double);
173  pm = ARRAY_1D(NMAX_POINT, double);
174  u = ARRAY_2D(NMAX_POINT, NVAR, double);
175  uh = ARRAY_2D(NMAX_POINT, NVAR, double);
176 
177  hp = ARRAY_1D(NMAX_POINT, double);
178  hm = ARRAY_1D(NMAX_POINT, double);
179  lambda_max = ARRAY_1D(NMAX_POINT, double);
180  lambda_min = ARRAY_1D(NMAX_POINT, double);
181  }
182  rhs = state->rhs;
183 
184 #if RECONSTRUCT_4VEL
185  ConvertTo3vel (state->v, beg-1, end+1);
186  ConvertTo3vel (state->vp, beg, end);
187  ConvertTo3vel (state->vm, beg, end);
188 #endif
189 
190 /* --------------------------------
191  make fluxes
192  -------------------------------- */
193 
194  PrimToCons (state->v, u, beg, end);
195  PrimToCons (state->vp, state->up, beg, end);
196  PrimToCons (state->vm, state->um, beg, end);
197 
198  Enthalpy (state->vp, hp, beg, end);
199  Enthalpy (state->vm, hm, beg, end);
200 
201  Flux(state->up, state->vp, hp, fp, pp, beg, end);
202  Flux(state->um, state->vm, hm, fm, pm, beg, end);
203 
204  dx = grid[g_dir].dx; A = grid[g_dir].A; dV = grid[g_dir].dV;
205 
206 /* ---------------------------------------------------
207  Compute passive scalar fluxes
208  --------------------------------------------------- */
209 
210  #if NFLX != NVAR
211  for (i = beg; i <= end; i++){
212  vp = state->vp[i]; vm = state->vm[i];
213  #if NSCL > 0
214  for (nv = NFLX; nv < (NFLX + NSCL); nv++){
215  fp[i][nv] = fp[i][RHO]*vp[nv];
216  fm[i][nv] = fm[i][RHO]*vm[nv];
217  }
218  #endif
219  }
220  #endif
221 
222  #if CHAR_UPWIND == YES
223  MAX_CH_SPEED (state->v, lambda_min, lambda_max, beg, end);
224  #endif
225 
226 /* ------------------------------------------------------
227  Initialize right-hand-side with flux difference
228  ------------------------------------------------------ */
229 
230  dt = 0.5*g_dt;
231  for (i = beg; i <= end; i++) {
232 
233  dtdx = dt/dx[i];
234  vp = state->vp[i];
235  vm = state->vm[i];
236  #if GEOMETRY == CARTESIAN
237  for (nv = NVAR; nv--; ) rhs[i][nv] = -dtdx*(fp[i][nv] - fm[i][nv]);
238  #elif GEOMETRY == CYLINDRICAL
239  dtdV = dt/dV[i];
240  for (nv = NVAR; nv--; ) {
241  rhs[i][nv] = -dtdV*(A[i]*fp[i][nv] - A[i-1]*fm[i][nv]);
242  }
243  #ifdef GLM_MHD
244  rhs[i][BXn] = -dtdx*(vp[PSI_GLM] - vm[PSI_GLM]);
245  #endif
246  #endif
247  rhs[i][MXn] -= dtdx*(pp[i] - pm[i]);
248 
249  /* ---- add Powell source term ---- */
250 
251  #if (PHYSICS == MHD || PHYSICS == RMHD) && (DIVB_CONTROL == EIGHT_WAVES)
252  dBx = (vp[BXn]*A[i] - vm[BXn]*A[i-1])/dV[i];
253  EXPAND(rhs[i][BX1] -= dt*state->v[i][VX1]*dBx; ,
254  rhs[i][BX2] -= dt*state->v[i][VX2]*dBx; ,
255  rhs[i][BX3] -= dt*state->v[i][VX3]*dBx;)
256  #endif
257 
258  /* ---- Extended GLM source term ---- */
259 
260  #ifdef GLM_MHD
261  #if GLM_EXTENDED == YES
262  dBx = (vp[BXn]*A[i] - vm[BXn]*A[i-1])/dV[i];
263  EXPAND(rhs[i][MX1] -= dt*state->v[i][BX1]*dBx; ,
264  rhs[i][MX2] -= dt*state->v[i][BX2]*dBx; ,
265  rhs[i][MX3] -= dt*state->v[i][BX3]*dBx;)
266  #if HAVE_ENERGY
267  dpsi = (vp[PSI_GLM] - vm[PSI_GLM])/dx[i];
268  rhs[i][ENG] -= dt*state->v[i][BXn]*dpsi;
269  #endif
270  #endif
271  #endif
272  }
273 
274 /* ------------------------------------------------
275  Source terms: geometry
276  ------------------------------------------------ */
277 
278  #if GEOMETRY == CYLINDRICAL && COMPONENTS == 3
279  if (g_dir == IDIR) {
280  double vB, vel2, g_2, r_1, *uc;
281 
282  for (i = beg; i <= end; i++){
283  r_1 = grid[IDIR].r_1[i];
284  vc = state->v[i];
285  uc = u[i];
286 
287  vB = EXPAND(vc[VX1]*vc[BX1], + vc[VX2]*vc[BX2], + vc[VX3]*vc[BX3]);
288  vel2 = EXPAND(vc[VX1]*vc[VX1], + vc[VX2]*vc[VX2], + vc[VX3]*vc[VX3]);
289  g_2 = 1.0 - vel2;
290 
291  rhs[i][iMR] += dt*((uc[MX3]*vc[VX3] - (vc[BX3]*g_2 + vB*vc[VX3])*vc[BX3]))*r_1;
292 
293  rhs[i][iMPHI] = -dtdV*(A[i]*A[i]*fp[i][iMPHI] - A[i-1]*A[i-1]*fm[i][iMPHI]);
294  rhs[i][iMPHI] *= r_1;
295 
296  rhs[i][iBPHI] -= dt*(vc[VX3]*uc[BX1] - uc[BX3]*vc[VX1])*r_1;
297  }
298  }
299  #endif
300 
301 /* ------------------------------------------------
302  Update solution vector to get u(n+1/2)
303  ------------------------------------------------ */
304 
305  for (i = beg; i <= end; i++) {
306  for (nv = NVAR; nv--; ) uh[i][nv] = u[i][nv] + rhs[i][nv];
307  }
308 
309 /* -------------------------------------------------
310  check if U(x_i, n+1/2) has physical meaning.
311  Revert to 1st order otherwise.
312  ------------------------------------------------- */
313 
314  if (ConsToPrim (uh, state->vh, beg, end, state->flag) != 0){
315  for (i = beg; i <= end; i++){
316  if (state->flag[i] != 0){
317  for (nv = NVAR; nv--; ){
318  uh[i][nv] = u[i][nv];
319  state->vh[i][nv] = state->v[i][nv];
320  state->vp[i][nv] = state->vm[i][nv] = state->v[i][nv];
321  }
322  #ifdef STAGGERED_MHD
323  state->vp[i][BXn] = state->bn[i];
324  state->vm[i][BXn] = state->bn[i-1];
325  #endif
326  }
327  }
328  }
329 
330 /* ----------------------------------------------------
331  evolve interface values accordingly
332  ---------------------------------------------------- */
333 
334  for (i = beg; i <= end; i++) {
335  vp = state->vp[i];
336  vm = state->vm[i];
337  vc = state->v[i];
338 
339  for (nv = NVAR; nv--; ){
340 /*
341  dvp[nv] = dvm[nv] = state->vh[i][nv] - state->v[i][nv];
342 
343  P0 = 1.5*vc[nv] - 0.25*(vp[nv] + vm[nv]);
344  P1 = vp[nv] - vm[nv];
345  P2 = 3.0*(vp[nv] - 2.0*vc[nv] + vm[nv]);
346 
347  dvp[nv] += P0 + 0.5*P1 + 0.25*P2 - vc[nv];
348  dvm[nv] += P0 - 0.5*P1 + 0.25*P2 - vc[nv];
349 */
350  P1 = vp[nv] - vm[nv];
351  dvp[nv] = dvm[nv] = state->vh[i][nv] - state->v[i][nv];
352  dvp[nv] += 0.5*P1;
353  dvm[nv] -= 0.5*P1;
354  }
355 
356  #if CHAR_UPWIND == YES
357 
358  if (lambda_min[i] > 0.0)
359  for (nv = NVAR; nv--; ) dvm[nv] = 0.0;
360  else if (lambda_max[i] < 0.0)
361  for (nv = NVAR; nv--; ) dvp[nv] = 0.0;
362 
363  #endif
364 
365  for (nv = NVAR; nv--; ){
366  vp[nv] = vc[nv] + dvp[nv];
367  vm[nv] = vc[nv] + dvm[nv];
368  }
369  }
370  #ifdef STAGGERED_MHD
371  for (i = beg-1; i <= end; i++){
372  state->vp[i][BXn] = state->vm[i+1][BXn] = state->bn[i];
373  }
374  #endif
375 
376  CheckPrimStates (state->vm, state->vp, state->v, beg, end);
377 }
378 #undef CHAR_UPWIND
379 #endif
380 
MX3
#define MX3
Definition: mod_defs.h:22
iMPHI
#define iMPHI
Definition: mod_defs.h:71
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
MX1
#define MX1
Definition: mod_defs.h:20
STATE_1D::vh
double ** vh
Primitive state at n+1/2 (only for one step method)
Definition: structs.h:162
VX2
#define VX2
Definition: mod_defs.h:29
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
Flux
void Flux(double **u, double **w, double *a2, double **fx, double *p, int beg, int end)
Definition: fluxes.c:23
STATE_1D::rhs
double ** rhs
Conservative right hand side.
Definition: structs.h:163
RHO
#define RHO
Definition: mod_defs.h:19
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
Enthalpy
void Enthalpy(double **v, real *h, int beg, int end)
Definition: eos.c:48
g_dt
double g_dt
The current integration time step.
Definition: globals.h:118
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
GRID::dx
double * dx
Definition: structs.h:83
VX1
#define VX1
Definition: mod_defs.h:28
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
iMR
#define iMR
Definition: mod_defs.h:69
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
CELL_CENTER
#define CELL_CENTER
Definition: pluto.h:205
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
ConsToPrim
int ConsToPrim(double **ucons, double **uprim, int ibeg, int iend, unsigned char *flag)
Definition: mappers.c:89
GRID
Definition: structs.h:78
NSCL
#define NSCL
Definition: pluto.h:572
MX2
#define MX2
Definition: mod_defs.h:21
Sph_disk.vel2
list vel2
Definition: Sph_disk.py:39
iBPHI
#define iBPHI
Definition: mod_defs.h:152
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
ConvertTo3vel
void ConvertTo3vel(double **, int, int)
Definition: four_vel.c:25
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MXn
int MXn
Definition: globals.h:74
PrimSource
void PrimSource(const State_1D *, int, int, double *, double *, double **, Grid *)
Definition: prim_eqn.c:71
GRID::x
double * x
Definition: structs.h:80
HancockStep
void HancockStep(const State_1D *state, int beg, int end, Grid *grid)
Definition: hancock.c:136
GRID::r_1
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
STATE_1D::vm
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
BX3
#define BX3
Definition: mod_defs.h:27
pluto.h
PLUTO main header file.
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
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
fp
FILE * fp
Definition: analysis.c:7
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
STATE_1D::a2
double * a2
Sound speed squared.
Definition: structs.h:158
PrimRHS
void PrimRHS(double *, double *, double, double, double *)
Definition: prim_eqn.c:31
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
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87