PLUTO
Macros | Functions
hancock.c File Reference

MUSCL-Hancock predictor step. More...

#include "pluto.h"
Include dependency graph for hancock.c:

Go to the source code of this file.

Macros

#define CHAR_UPWIND   NO
 

Functions

void HancockStep (const State_1D *state, int beg, int end, Grid *grid)
 

Detailed Description

MUSCL-Hancock predictor step.

Use time-extrapolation to compute interface states and cell-centered value at the half-time step, $ V_{i,\pm}^{n+\HALF} = V_{i,\pm}^n + \partial_t V \Delta t^n/2 $

This is done using the one-dimensional primitive form ot the equations for the HD, RHD and MHD modules since we have at disposal $ \partial_t V = -A\partial V_x + S $.

Conversely, for relativistic MHD, we adopt the one-dimensional conservative form of the equations (conservative Hancock predictor) and compute the primitive values using $ \partial_t V \approx (V^{n+\HALF}_i-V^n_i)/(\Delta t/2) $. This requires taking the following steps:

  • Convert to conservative: vp -> up, vm -> um
  • Evolve u(n+1/2) = u(n) - dt/(2dx)*[F(vp) - F(vm)]
  • Convert to primitive u(n+1/2) -> v(n+1/2)
  • Add time incerement v(n+1/2)-v(n) to left/right states
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
June 24, 2015

Definition in file hancock.c.

Macro Definition Documentation

#define CHAR_UPWIND   NO

Function Documentation

void HancockStep ( const State_1D *  state,
int  beg,
int  end,
Grid *  grid 
)

Definition at line 136 of file hancock.c.

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 }
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
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
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
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
ConsToPrim
int ConsToPrim(double **ucons, double **uprim, int ibeg, int iend, unsigned char *flag)
Definition: mappers.c:89
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
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
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
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::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
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

Here is the call graph for this function:

Here is the caller graph for this function: