PLUTO
Macros | Functions
char_tracing.c File Reference
#include "pluto.h"
Include dependency graph for char_tracing.c:

Go to the source code of this file.

Macros

#define REF_STATE   3
 
#define REF_STATE   3
 

Functions

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

Macro Definition Documentation

#define REF_STATE   3

Definition at line 195 of file char_tracing.c.

#define REF_STATE   3

Definition at line 195 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 18 of file char_tracing.c.

23 {
24  int i, nv, k;
25  double dx, dtdx;
26  double dwh, d2w, dwp[NVAR], dwm[NVAR], dvp[NVAR], dvm[NVAR];
27  double nup=1.0, num=-1.0, nu[NVAR];
28  double Spp, Smm;
29  double *vc, *vp, *vm, **L, **R, *lambda;
30  static double **src;
31 
32  if (src == NULL){
33  src = ARRAY_2D(NMAX_POINT, NVAR, double);
34  }
35 
36  #if CHAR_LIMITING == NO
37  SoundSpeed2 (state->v, state->a2, state->h, beg, end, CELL_CENTER, grid);
38  #endif
39 
40  PrimSource (state, beg, end, state->a2, state->h, src, grid);
41  for (i = beg; i <= end; i++){
42 
43  dx = grid[g_dir].dx[beg];
44  dtdx = g_dt/dx;
45 
46  vc = state->v[i];
47  vp = state->vp[i];
48  vm = state->vm[i];
49  L = state->Lp[i];
50  R = state->Rp[i];
51  lambda = state->lambda[i];
52 
53  /* --------------------------------------------------------------
54  1. Compute eigenvectors and eigenvalues if not yet defined
55  -------------------------------------------------------------- */
56 
57  #if CHAR_LIMITING == NO
58  PrimEigenvectors (vc, state->a2[i], state->h[i], lambda, L, R);
59  #if NVAR != NFLX
60  for (k = NFLX; k < NVAR; k++) lambda[k] = vc[V1];
61  #endif
62  #endif
63  for (k = 0; k < NVAR; k++) nu[k] = dtdx*lambda[k];
64 
65  /* --------------------------------------------------------------
66  2. Obtain characteristic variable increments dwp and dwm.
67  -------------------------------------------------------------- */
68 
69  for (nv = NVAR; nv--; ) {
70  dvp[nv] = vp[nv] - vc[nv];
71  dvm[nv] = vm[nv] - vc[nv];
72  }
73  PrimToChar(L, dvp, dwp);
74  PrimToChar(L, dvm, dwm);
75 
76  /* ----------------------------------------------------------------
77  3. Initialize vp and vm to the reference state.
78  Since this is somewhat arbitrary we set it using REF_STATE:
79 
80  REF_STATE==1: use cell-center value;
81  REF_STATE==2: interpolated value at base time level
82  REF_STATE==3: traditional PPM reference state (fastest
83  wave), minimize the size of the term
84  subject to characteristic limiting.
85 
86  Passive scalars use always REF_STATE == 2.
87  ---------------------------------------------------------------- */
88 
89  #if REF_STATE == 1
90  for (nv = NFLX; nv--; ) vp[nv] = vm[nv] = vc[nv];
91  #elif REF_STATE == 3
92  nup = MAX(nu[1], 0.0); num = MIN(nu[0], 0.0);
93  for (nv = NFLX; nv--; ){
94  dwh = vp[nv] - vm[nv];
95  d2w = vp[nv] + vm[nv] - 2.0*vc[nv];
96  vp[nv] -= 0.5*nup*(dwh + d2w*(3.0 - 2.0*nup));
97  vm[nv] -= 0.5*num*(dwh - d2w*(3.0 + 2.0*num));
98  }
99  #endif
100 
101  /* ------------------------------------------------------------------
102  4. Compute L/R states in primitive variables:
103 
104  - evolve characteristic variable increments by dt/2;
105  - discard contributions from waves not reaching the interface;
106  - project characteristic differences dwp and dwm onto right
107  eigenvectors
108  ------------------------------------------------------------------ */
109 
110  for (k = 0; k < NFLX; k++){
111  dwh = dwp[k] - dwm[k];
112  d2w = dwp[k] + dwm[k];
113  if (nu[k] >= 0.0) {
114  #if REF_STATE == 1
115  Spp = dwp[k] - 0.5*nu[k]*(dwh + d2w*(3.0 - 2.0*nu[k]));
116  #elif REF_STATE == 2
117  Spp = - 0.5*nu[k]*(dwh + d2w*(3.0 - 2.0*nu[k]));
118  #elif REF_STATE == 3
119  Spp = -0.5*(nu[k]-nup)*(dwh + 3.0*d2w) + (nu[k]*nu[k] - nup*nup)*d2w;
120  #endif
121  for (nv = NFLX; nv--; ) vp[nv] += Spp*R[nv][k];
122  } else {
123  #if REF_STATE == 1
124  Smm = dwm[k] - 0.5*nu[k]*(dwh - d2w*(3.0 + 2.0*nu[k]));
125  #elif REF_STATE == 2
126  Smm = - 0.5*nu[k]*(dwh - d2w*(3.0 + 2.0*nu[k]));
127  #elif REF_STATE == 3
128  Smm = -0.5*(nu[k]-num)*(dwh - 3.0*d2w) + (nu[k]*nu[k] - num*num)*d2w;
129  #endif
130  for (nv = NFLX; nv--; ) vm[nv] += Smm*R[nv][k];
131  }
132  }
133 
134  /* -------------------------------------------------------
135  5. Add source term to L/R states
136  ------------------------------------------------------- */
137 
138  for (nv = NFLX; nv--; ){
139  dwh = 0.5*g_dt*src[i][nv];
140  vp[nv] += dwh;
141  vm[nv] += dwh;
142  }
143 
144  /* -------------------------------------------------------
145  6. Repeat construction for passive scalars
146  ------------------------------------------------------- */
147 
148  #if NVAR != NFLX
149  if (nu[NFLX] >= 0.0) { /* -- scalars all move at the flow speed -- */
150  for (k = NFLX; k < NVAR; k++){
151  dwh = dwp[k] - dwm[k];
152  d2w = dwp[k] + dwm[k];
153  vp[k] -= 0.5*nu[k]*(dwh + d2w*(3.0 - 2.0*nu[k]));
154  }
155  }else{
156  for (k = NFLX; k < NVAR; k++){
157  dwh = dwp[k] - dwm[k];
158  d2w = dwp[k] + dwm[k];
159  vm[k] -= 0.5*nu[k]*(dwh - d2w*(3.0 + 2.0*nu[k]));
160  }
161  }
162  #endif
163  }
164 
165  #ifdef STAGGERED_MHD
166  for (i = beg-1; i <= end; i++) {
167  state->vR[i][BXn] = state->vL[i][BXn] = state->bn[i];
168  }
169  #endif
170 
171  CheckPrimStates (state->vm, state->vp, state->v, beg, end);
172 
173  for (i = beg; i <= end; i++) {
174  vp = state->vp[i];
175  vm = state->vm[i];
176  for (nv = NVAR; nv--; ) {
177  state->vh[i][nv] = 0.5*(vp[nv] + vm[nv]);
178  }
179  }
180 }
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
STATE_1D::lambda
double ** lambda
Characteristic speed associated to Lp and Rp.
Definition: structs.h:156
k
int k
Definition: analysis.c:2
PrimSource
void PrimSource(const State_1D *, int, int, double *, double *, double **, Grid *)
Definition: prim_eqn.c:71
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:

Here is the caller graph for this function: