24 double *cmax,
Grid *grid)
44 int nv,
i,
j,
k, np_tot;
46 double **
v, **flux, *press;
47 double fs,
fp, fm, flx[
NFLX], dx;
49 static double *Fp, *Fm, **F, *a2, **u;
50 static double **Uave, **Vave, **lp;
52 static double *psim, *Bm;
53 double (*REC)(
double *, double, int);
59 #if TIME_STEPPING != RK3 && TIME_STEPPING != SSP_RK4
60 print1 (
"! Finite Difference schemes work with RK3 only \n");
63 #if GEOMETRY != CARTESIAN
64 print1 (
"! Finite Difference schemes work in Cartesian coordinates only\n");
67 #if PHYSICS == RMHD || PHYSICS == RHD
68 print1 (
"! Finite difference schemes work only for HD od MHD modules\n");
72 print1 (
"! Finite difference schemes work only with cell-centered schemes\n");
75 #if PHYSICS == MHD && BACKGROUND_FIELD == YES
76 print1 (
"! Background field splitting not supported with FD schemes\n");
87 #if RECONSTRUCTION == WENO3_FD
90 #elif RECONSTRUCTION == LIMO3_FD
93 #elif RECONSTRUCTION == WENOZ_FD
96 #elif RECONSTRUCTION == MP5_FD
116 press = state->
press;
130 Flux (u, v, a2, NULL, F, press, 0, np_tot-1);
132 Flux (u, v, a2, F, press, 0, np_tot-1);
135 for (i = 0; i < np_tot; i++){
136 F[
i][
MXn] += press[
i];
138 #if EOS != ISOTHERMAL
141 fs = g_isoSoundSpeed*g_isoSoundSpeed;
143 fs = fabs(v[i][
VXn])/sqrt(fs);
153 for (i = beg; i <= end; i++){
155 fm = REC(Fm, dx, np_tot - i - 2);
156 #if SHOCK_FLATTENING == MULTID
164 #if COMPUTE_DIVB == YES
165 state->
bn[
i] = Bm[
i];
168 #if COMPUTE_DIVB == YES
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]);
203 for (i = beg; i <= end; i++){
210 for (k = Kmax; k--; ){
211 fs =
MAX(fs, fabs(state->
lmax[k]));
213 cmax[0] =
MAX(cmax[0], fs);
230 for (i = beg; i <= end; i++){
234 for (k = 0; k < Kmax; k++){
236 for (j = i-S; j <= i+S; j++) {
238 for (nv =
NFLX; nv--; ){
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]);
243 #if SHOCK_FLATTENING == MULTID
248 flx[
k] = REC(Fp, dx, i) + REC(Fm, dx,i);
251 for (nv = 0; nv <
NFLX; nv++){
253 for (k=0; k < Kmax; k++) fs += flx[k]*R[nv][k];
258 flux[
i][
BXn] = psim[
i];
275 if (flux[i][
RHO] >= 0.0){
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];
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);
304 double *x1, *x2, *x3;
305 static double **
v, **lambda, *a2;
313 for (nv =
NVAR; nv--; ) state->
lmax[nv] = 0.0;
324 for (nv =
NFLX; nv--; )
325 state->
lmax[nv] =
MAX(fabs(lambda[i][nv]), state->
lmax[nv]);
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];
double WENO3_Reconstruct(double *F, double dx, int j)
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
double * lmax
Define the maximum k-characteristic speed over the domain.
double ** flux
upwind flux computed with the Riemann solver
void FD_GetMaxEigenvalues(const Data *d, State_1D *state, Grid *grid)
void Eigenvalues(double **v, double *csound2, double **lambda, int beg, int end)
void print1(const char *fmt,...)
void Flux(double **u, double **w, double *a2, double **fx, double *p, int beg, int end)
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
double WENOZ_Reconstruct(double *F, double dx, int j)
double **** Vc
The main four-index data array used for cell-centered primitive variables.
double ** vR
Primitive variables to the right of the interface, .
#define FLAG_MINMOD
Reconstruct using MINMOD limiter.
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
double MP5_Reconstruct(double *F, double dx, int j)
double g_maxMach
The maximum Mach number computed during integration.
int g_dir
Specifies the current sweep or direction of integration.
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
#define ARRAY_1D(nx, type)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
void ConsEigenvectors(double *u, double *v, double a2, double **LL, double **RR, double *lambda)
double * bn
Face magentic field, bn = bx(i+1/2)
double LIMO3_Reconstruct(double *v, double dx, int i)
double ** vL
Primitive variables to the left of the interface, .
double * press
Upwind pressure term computed with the Riemann solver.
#define ARRAY_2D(nx, ny, type)
double *** Rp
Left and right primitive eigenvectors.
int np_tot
Total number of points in the local domain (boundaries included).
double glm_ch
The propagation speed of divergence error.
#define QUIT_PLUTO(e_code)
double LIN_Reconstruct(double *F, double dx, int j)
void GLM_ComputeDivB(const State_1D *state, Grid *grid)
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
void FD_Flux(const State_1D *state, int beg, int end, double *cmax, Grid *grid)