PLUTO
Macros | Functions
char_tracing.c File Reference

Compute time-centered interface states using characteristic tracing. More...

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

Go to the source code of this file.

Macros

#define CHTR_REF_STATE   2
 

Functions

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

Detailed Description

Compute time-centered interface states using characteristic tracing.

Advance 1-D left and right interface states, previously computed with any of the States functions, to the half time level (n+1/2) by extrapolating characteristic variables. This is done using an upwind selection rule that discards waves not reaching the interface in dt/2.

References:

  • "The PLUTO Code for Adaptive Mesh Refinement in Astrophysical Fluid Dynamics"
    Mignone et al, ApJS (2012) 198, 7M
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Jan 28, 2014

Definition in file char_tracing.c.

Macro Definition Documentation

#define CHTR_REF_STATE   2

Flag to control the choice of the reference state in the upwind selection rule. Set CHTR_REF_STATE to 1,2,3 to use

  • cell centered value (1),
  • interpolated states (2),
  • fastest wave (3, the usual PPM rescription).

Definition at line 37 of file char_tracing.c.

Function Documentation

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

Compute interface states using characteristic tracing step.

Parameters
[in]statepointer to a State_1D structure
[in]beginitial index of computation
[in]endfinal index of computation
[in]gridpointer to an array of Grid structures
Returns
This function has no return value.

Tasks are numbered below.

1) Compute eigenvectors and eigenvalues if not yet defined.

2) Obtain characteristic variable increments dwp and dwm.

3) Initialize vp and vm to the reference state. Since this is somewhat arbitrary we use the value of CHTR_REF_STATE to select one of the following cases:

  • CHTR_REF_STATE==1: use cell-center value;
  • CHTR_REF_STATE==2: interpolated value at base time level
  • CHTR_REF_STATE==3: traditional PPM reference state (fastest wave), minimize the size of the term subject to characteristic limiting.

Passive scalars use always CHTR_REF_STATE == 2.

4) Compute left and right states in primitive variables. This step also depends on the value of CHTR_REF_STATE and include:

  • evolve characteristic variable increments by dt/2;
  • discard contributions from waves not reaching the interface;
  • project characteristic differences dwp and dwm onto right eigenvectors

5) Add source term to L/R states

6) Repeat construction for passive scalars

Definition at line 46 of file char_tracing.c.

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 }
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
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
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
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
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
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
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
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

Here is the call graph for this function: