PLUTO
Functions
fd_flux.c File Reference

Compute finite difference fluxes. More...

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

Go to the source code of this file.

Functions

void FD_Flux (const State_1D *state, int beg, int end, double *cmax, Grid *grid)
 
void FD_GetMaxEigenvalues (const Data *d, State_1D *state, Grid *grid)
 

Detailed Description

Compute finite difference fluxes.

The function FD_Flux() computes a high-order approximation to the flux divergence using WENO or Monotonicity-preserving reconstruction in characteristic variables.

It is called as a regular Riemann solver function except that the left and right states do not need to be defined since high-order reconstruction is done on the characteristic projection of the cell-centered fluxes.

Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
P. Tzeferacos (petro.nosp@m.s.tz.nosp@m.efera.nosp@m.cos@.nosp@m.ph.un.nosp@m.ito..nosp@m.it)
Date
April 02, 2015

Definition in file fd_flux.c.

Function Documentation

void FD_Flux ( const State_1D *  state,
int  beg,
int  end,
double *  cmax,
Grid *  grid 
)

Compute interface flux by suitable high-order finite difference non-oscillatory interpolants.

Parameters
[in]statepointer to State_1D structure
[in]beginitial index of computation
[in]endfinal index of computation
[out]cmaxarray of maximum characteristic speeds
[in]gridpointer to an array of Grid structures

Reference

  • "High-order conservative finite difference GLM-MHD schemes for cell-centered MHD"
    Mignone et al., JCP (2010), 229, 5896
Returns
This function has no return value.

Definition at line 23 of file fd_flux.c.

43 {
44  int nv, i, j, k, np_tot;
45  int Kmax, S;
46  double **v, **flux, *press;
47  double fs, fp, fm, flx[NFLX], dx;
48 
49  static double *Fp, *Fm, **F, *a2, **u;
50  static double **Uave, **Vave, **lp;
51  double **L, **R;
52  static double *psim, *Bm;
53  double (*REC)(double *, double, int);
54 
55 /* ------------------------------------------------------------------
56  Check compatibility with other modules
57  ------------------------------------------------------------------ */
58 
59  #if TIME_STEPPING != RK3 && TIME_STEPPING != SSP_RK4
60  print1 ("! Finite Difference schemes work with RK3 only \n");
61  QUIT_PLUTO(1);
62  #endif
63  #if GEOMETRY != CARTESIAN
64  print1 ("! Finite Difference schemes work in Cartesian coordinates only\n");
65  QUIT_PLUTO(1);
66  #endif
67  #if PHYSICS == RMHD || PHYSICS == RHD
68  print1 ("! Finite difference schemes work only for HD od MHD modules\n");
69  QUIT_PLUTO(1);
70  #endif
71  #ifdef STAGGERED_MHD
72  print1 ("! Finite difference schemes work only with cell-centered schemes\n");
73  QUIT_PLUTO(1);
74  #endif
75  #if PHYSICS == MHD && BACKGROUND_FIELD == YES
76  print1 ("! Background field splitting not supported with FD schemes\n");
77  QUIT_PLUTO(1);
78  #endif
79 
80 /* --------------------------------------------------------------
81  - Define pointer to reconstruction function
82  - Determine the stencil for interpolation
83  For a given S, the stencil is: i-S <= i <= i+S+1
84  -------------------------------------------------------------- */
85 
86  dx = grid[g_dir].dx[beg];
87  #if RECONSTRUCTION == WENO3_FD
88  REC = WENO3_Reconstruct;
89  S = 1;
90  #elif RECONSTRUCTION == LIMO3_FD
91  REC = LIMO3_Reconstruct;
92  S = 1;
93  #elif RECONSTRUCTION == WENOZ_FD
94  REC = WENOZ_Reconstruct;
95  S = 2;
96  #elif RECONSTRUCTION == MP5_FD
97  REC = MP5_Reconstruct;
98  S = 2;
99  #endif
100 
101  if (F == NULL){
102  Fp = ARRAY_1D(NMAX_POINT, double);
103  Fm = ARRAY_1D(NMAX_POINT, double);
104  F = ARRAY_2D(NMAX_POINT, NVAR, double);
105  lp = ARRAY_2D(NMAX_POINT, NVAR, double);
106  psim = ARRAY_1D(NMAX_POINT, double);
107  Bm = ARRAY_1D(NMAX_POINT, double);
108  a2 = ARRAY_1D(NMAX_POINT, double);
109  Uave = ARRAY_2D(NMAX_POINT, NVAR, double);
110  Vave = ARRAY_2D(NMAX_POINT, NVAR, double);
111  u = ARRAY_2D(NMAX_POINT, NVAR, double);
112  }
113 
114  np_tot = grid[g_dir].np_tot;
115  flux = state->flux;
116  press = state->press;
117  v = state->v;
118 
119 /* ------------------------------------------------------
120  Compute cell-centered fluxes. Pressure is added to
121  the normal component of momentum.
122  To save computational time, we also solve the linear
123  GLM-MHD subsystem here and use the solution in the
124  nonlinear fluxes later.
125  ------------------------------------------------------ */
126 
127  PrimToCons (v, u, 0, np_tot - 1);
128  SoundSpeed2 (v, a2, NULL, 0, np_tot-1, CELL_CENTER, grid);
129  #if PHYSICS == MHD
130  Flux (u, v, a2, NULL, F, press, 0, np_tot-1);
131  #else
132  Flux (u, v, a2, F, press, 0, np_tot-1);
133  #endif
134 
135  for (i = 0; i < np_tot; i++){
136  F[i][MXn] += press[i];
137  press[i] = 0.0;
138  #if EOS != ISOTHERMAL
139  fs = g_gamma*v[i][PRS]/v[i][RHO];
140  #else
141  fs = g_isoSoundSpeed*g_isoSoundSpeed;
142  #endif
143  fs = fabs(v[i][VXn])/sqrt(fs);
144  g_maxMach = MAX(g_maxMach, fs);
145  #ifdef GLM_MHD
146  j = np_tot - 1 - i;
147  Fp[i] = 0.5*(u[i][PSI_GLM] + glm_ch*u[i][BXn]);
148  Fm[j] = 0.5*(u[i][PSI_GLM] - glm_ch*u[i][BXn]);
149  #endif
150  }
151 
152  #ifdef GLM_MHD
153  for (i = beg; i <= end; i++){
154  fp = REC(Fp, dx, i);
155  fm = REC(Fm, dx, np_tot - i - 2);
156  #if SHOCK_FLATTENING == MULTID
157  if (state->flag[i] & FLAG_MINMOD){
158  fp = WENO3_Reconstruct(Fp, dx, i);
159  fm = WENO3_Reconstruct(Fm, dx, np_tot - i - 2);
160  }
161  #endif
162  Bm[i] = (fp - fm)/glm_ch;
163  psim[i] = fp + fm;
164  #if COMPUTE_DIVB == YES
165  state->bn[i] = Bm[i];
166  #endif
167  }
168  #if COMPUTE_DIVB == YES
169  if (g_intStage == 1) GLM_ComputeDivB(state, grid);
170  #endif
171 
172  Kmax = KPSI_GLMM;
173  #else
174  Kmax = NFLX;
175  #endif
176 
177 /* ----------------------------------------------------
178  Compute average state
179  ---------------------------------------------------- */
180 
181  for (i = beg; i <= end; i++){
182  for (nv = NFLX; nv--; ){
183  Uave[i][nv] = 0.5*(u[i][nv] + u[i + 1][nv]);
184  Vave[i][nv] = 0.5*(v[i][nv] + v[i + 1][nv]);
185  }
186 
187  #ifdef GLM_MHD
188  Uave[i][BXn] = Vave[i][BXn] = Bm[i];
189  Uave[i][PSI_GLM] = Vave[i][PSI_GLM] = psim[i];
190  #endif
191  }
192 
193 /* ---------------------------------------
194  compute sound speed
195  --------------------------------------- */
196 
197  SoundSpeed2 (Vave, a2, NULL, beg, end, FACE_CENTER, grid);
198 
199 /* ----------------------------------------------------
200  Compute eigenvectors and eigenvalues
201  ---------------------------------------------------- */
202 
203  for (i = beg; i <= end; i++){
204  L = state->Lp[i];
205  R = state->Rp[i];
206 
207  ConsEigenvectors(Uave[i], Vave[i], a2[i], L, R, lp[i]);
208 
209  fs = 0.0;
210  for (k = Kmax; k--; ){
211  fs = MAX(fs, fabs(state->lmax[k]));
212  }
213  cmax[0] = MAX(cmax[0], fs);
214  cmax[i] = fs;
215  }
216 
217 /* ---------------------------------------------------------
218  Main spatial loop.
219  At each interface construct, for each characteristic
220  field k, the positive and negative part of the flux:
221 
222  F_{+,j} = 1/2 L.[F(j) + a*u(j)]
223  F_{-,j} = 1/2 L.[F(j') - a*u(j')]
224 
225  a is the global maximum eigenvalue for the k-th
226  characteristic.
227  Reconstruct F_{+,j} and F_{-,j}.
228  --------------------------------------------------------- */
229 
230  for (i = beg; i <= end; i++){
231  L = state->Lp[i];
232  R = state->Rp[i];
233 
234  for (k = 0; k < Kmax; k++){
235  fs = state->lmax[k];
236  for (j = i-S; j <= i+S; j++) {
237  Fp[j] = Fm[j] = 0.0;
238  for (nv = NFLX; nv--; ){ /* -- LF flux splitting -- */
239  Fp[j] += 0.5*L[k][nv]*(F[j][nv] + fs*u[j][nv]);
240  Fm[j] += 0.5*L[k][nv]*(F[2*i-j+1][nv] - fs*u[2*i-j+1][nv]);
241  }
242  }
243  #if SHOCK_FLATTENING == MULTID
244  if ( (state->flag[i] & FLAG_MINMOD) || (state->flag[i+1] & FLAG_MINMOD)){
245  fs = LIN_Reconstruct(Fp, dx, i) + LIN_Reconstruct(Fm, dx, i);
246  }else
247  #endif
248  flx[k] = REC(Fp, dx, i) + REC(Fm, dx,i);
249  }
250 
251  for (nv = 0; nv < NFLX; nv++){
252  fs = 0.0;
253  for (k=0; k < Kmax; k++) fs += flx[k]*R[nv][k];
254  flux[i][nv] = fs;
255  }
256 
257  #ifdef GLM_MHD
258  flux[i][BXn] = psim[i];
259  flux[i][PSI_GLM] = glm_ch*glm_ch*Bm[i];
260  #endif
261 
262  /* -------------------------------------------
263  Passive scalars are treated separately in
264  a way similar to the finite volume case:
265 
266  F(T) = F(\rho)*TL if F(rho) >= 0
267  F(T) = F(\rho)*TR if F(rho) < 0
268 
269  (inside adv_flux).
270  This seems to preserve better the
271  original tracer bounds (e.g., 0<T<1).
272  ------------------------------------------ */
273 
274 #if NSCL > 0
275  if (flux[i][RHO] >= 0.0){
276  NSCL_LOOP(k){
277  for (j = i-S; j <= i+S; j++) Fp[j] = v[j][k];
278  state->vL[i][k] = REC(Fp, dx, i);
279  state->vR[i][k] = v[i+1][k];
280  }
281  }else{
282  NSCL_LOOP(k){
283  for (j = i-S; j <= i+S; j++) Fm[j] = v[2*i-j+1][k];
284  state->vL[i][k] = v[i][k];
285  state->vR[i][k] = REC(Fm, dx, i);
286  }
287  }
288 #endif
289  }
290 }
MAX
#define MAX(a, b)
Definition: macros.h:101
WENO3_Reconstruct
double WENO3_Reconstruct(double *F, double dx, int j)
Definition: fd_reconstruct.c:192
STATE_1D::v
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
g_gamma
double g_gamma
Definition: globals.h:112
STATE_1D::lmax
double * lmax
Define the maximum k-characteristic speed over the domain.
Definition: structs.h:157
STATE_1D::flux
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
STATE_1D::Lp
double *** Lp
Definition: structs.h:155
Flux
void Flux(double **u, double **w, double *a2, double **fx, double *p, int beg, int end)
Definition: fluxes.c:23
RHO
#define RHO
Definition: mod_defs.h:19
g_intStage
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
Definition: globals.h:98
WENOZ_Reconstruct
double WENOZ_Reconstruct(double *F, double dx, int j)
Definition: fd_reconstruct.c:140
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
func_param::v
double v[NVAR]
Definition: eos.h:106
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
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
MP5_Reconstruct
double MP5_Reconstruct(double *F, double dx, int j)
Definition: fd_reconstruct.c:4
g_maxMach
double g_maxMach
The maximum Mach number computed during integration.
Definition: globals.h:119
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
FACE_CENTER
#define FACE_CENTER
Definition: pluto.h:206
j
int j
Definition: analysis.c:2
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
k
int k
Definition: analysis.c:2
MXn
int MXn
Definition: globals.h:74
VXn
int VXn
Definition: globals.h:73
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
ConsEigenvectors
void ConsEigenvectors(double *u, double *v, double a2, double **LL, double **RR, double *lambda)
Definition: eigenv.c:279
STATE_1D::bn
double * bn
Face magentic field, bn = bx(i+1/2)
Definition: structs.h:165
LIMO3_Reconstruct
double LIMO3_Reconstruct(double *v, double dx, int i)
Definition: fd_reconstruct.c:103
fp
FILE * fp
Definition: analysis.c:7
STATE_1D::vL
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
STATE_1D::press
double * press
Upwind pressure term computed with the Riemann solver.
Definition: structs.h:164
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
GRID::np_tot
int np_tot
Total number of points in the local domain (boundaries included).
Definition: structs.h:100
glm_ch
double glm_ch
The propagation speed of divergence error.
Definition: glm.c:21
NVAR
#define NVAR
Definition: pluto.h:609
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
LIN_Reconstruct
double LIN_Reconstruct(double *F, double dx, int j)
Definition: fd_reconstruct.c:242
GLM_ComputeDivB
void GLM_ComputeDivB(const State_1D *state, Grid *grid)
Definition: glm.c:411
KPSI_GLMM
Definition: mod_defs.h:75

Here is the call graph for this function:

Here is the caller graph for this function:

void FD_GetMaxEigenvalues ( const Data *  d,
State_1D *  state,
Grid *  grid 
)

Definition at line 293 of file fd_flux.c.

302 {
303  int i, j, k, nv;
304  double *x1, *x2, *x3;
305  static double **v, **lambda, *a2;
306 
307  if (v == NULL){
308  v = ARRAY_2D(NMAX_POINT, NFLX, double);
309  lambda = ARRAY_2D(NMAX_POINT, NFLX, double);
310  a2 = ARRAY_1D(NMAX_POINT, double);
311  }
312 
313  for (nv = NVAR; nv--; ) state->lmax[nv] = 0.0;
314 
315  x1 = grid[IDIR].x;
316  x2 = grid[JDIR].x;
317  x3 = grid[KDIR].x;
318 
319  KDOM_LOOP(k) JDOM_LOOP(j){
320  IDOM_LOOP(i) for (nv = NVAR; nv--; ) v[i][nv] = d->Vc[nv][k][j][i];
321  SoundSpeed2 (v, a2, NULL, IBEG, IEND, CELL_CENTER, grid);
322  Eigenvalues (v, a2, lambda, IBEG, IEND);
323  IDOM_LOOP(i){
324  for (nv = NFLX; nv--; )
325  state->lmax[nv] = MAX(fabs(lambda[i][nv]), state->lmax[nv]);
326  }
327  }
328 
329  #ifdef PARALLEL
330  MPI_Allreduce (state->lmax, lambda[0], NFLX, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
331  for (nv = 0; nv < NFLX; nv++) state->lmax[nv] = lambda[0][nv];
332  #endif
333 }
MAX
#define MAX(a, b)
Definition: macros.h:101
KDOM_LOOP
#define KDOM_LOOP(k)
Definition: macros.h:36
STATE_1D::lmax
double * lmax
Define the maximum k-characteristic speed over the domain.
Definition: structs.h:157
Eigenvalues
void Eigenvalues(double **v, double *csound2, double **lambda, int beg, int end)
Definition: eigenv.c:67
func_param::v
double v[NVAR]
Definition: eos.h:106
KDIR
#define KDIR
Definition: pluto.h:195
SoundSpeed2
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
Definition: eos.c:16
JDOM_LOOP
#define JDOM_LOOP(j)
Definition: macros.h:35
CELL_CENTER
#define CELL_CENTER
Definition: pluto.h:205
NFLX
#define NFLX
Definition: mod_defs.h:32
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
IEND
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:37
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
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
i
int i
Definition: analysis.c:2
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
IDOM_LOOP
#define IDOM_LOOP(i)
Definition: macros.h:34
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35

Here is the call graph for this function:

Here is the caller graph for this function: