103 double dt2, inv_dtp, *inv_dl;
105 static unsigned char *flagm, *flagp;
110 static double **dtdV, **dcoeff, ***
T;
111 double *dtdV2, **rhs;
118 #if !(GEOMETRY == CARTESIAN || GEOMETRY == CYLINDRICAL)
119 print1 (
"! AdvanceStep(): ");
120 print1 (
"CTU only works in Cartesian or cylindrical coordinates\n");
165 #if (PARABOLIC_FLUX & EXPLICIT)
168 #if THERMAL_CONDUCTION == EXPLICIT
189 #if (SHOCK_FLATTENING == MULTID) || (ENTROPY_SWITCH)
200 #if THERMAL_CONDUCTION == EXPLICIT
201 TOT_LOOP(k,j,i) T[k][j][i] = d->Vc[PRS][k][j][i]/d->Vc[
RHO][k][j][i];
208 #if (RESISTIVITY == EXPLICIT) && (defined STAGGERED_MHD)
218 VAR_LOOP(nv) state.v[i][nv] = d->Vc[nv][k][j][i];
224 nv = NX1_TOT*
sizeof(double);
227 memcpy(Bs0[
BX2s][k][j], d->
Vs[
BX2s][k][j], nv); ,
228 memcpy(Bs0[
BX3s][k][j], d->
Vs[
BX3s][k][j], nv);)
242 #if (RESISTIVITY == EXPLICIT) && !(defined STAGGERED_MHD)
263 for (*in = 0; (*in) < indx.
ntot; (*in)++) {
264 VAR_LOOP(nv) state.v[*in][nv] = d->Vc[nv][k][j][i];
265 state.flag[*in] = d->flag[k][j][i];
273 #if !(PARABOLIC_FLUX & EXPLICIT)
276 Riemann (&state, indx.
beg - 1, indx.
end, Dts->
cmax, grid);
282 #if CTU_MHD_SOURCE == YES
317 for (nv =
NVAR; nv--; ) {
318 state.
rhs[indx.
beg - 1][nv] = 0.0;
319 state.
rhs[indx.
end + 1][nv] = 0.0;
322 for ((*in) = indx.
beg-1; (*in) <= indx.
end+1; (*in)++) {
323 for (nv =
NVAR; nv--; ){
324 dU[
k][
j][
i][nv] = state.
rhs[*in][nv];
325 UH[
k][
j][
i][nv] = d->
Uc[
k][
j][
i][nv] + state.
rhs[*in][nv];
326 state.
up[*in][nv] -= state.
rhs[*in][nv];
327 state.
um[*in][nv] -= state.
rhs[*in][nv];
330 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++) {
331 for (nv =
NVAR; nv--; ){
332 dU[
k][
j][
i][nv] += state.
rhs[*in][nv];
333 UH[
k][
j][
i][nv] += state.
rhs[*in][nv];
334 state.
up[*in][nv] -= state.
rhs[*in][nv];
335 state.
um[*in][nv] -= state.
rhs[*in][nv];
361 for ((*in) = 0; (*in) < indx.
ntot; (*in)++) {
362 for (nv =
NVAR; nv--; ) {
363 state.
vp[*in][nv] = state.
vm[*in][nv] = state.
vh[*in][nv] = state.
v[*in][nv];
366 for (*in = 0; (*in) < indx.
ntot-1; (*in)++) {
367 state.
vR[*in][
BXn] = state.
vL[*in][
BXn] = state.
bn[*in];
373 Riemann (&state, indx.
beg-1, indx.
end, Dts->
cmax, grid);
382 #if (VISCOSITY == EXPLICIT)
383 for ((*in) = 0; (*in) < indx.
ntot; (*in)++)
for (nv =
NVAR; nv--; )
384 state.par_src[*in][nv] = 0.0;
401 #if CTU_MHD_SOURCE == YES
406 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++) {
407 for (nv =
NVAR; nv--; ){
408 state.
up[*in][nv] -= state.
rhs[*in][nv];
409 state.
um[*in][nv] -= state.
rhs[*in][nv];
419 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++) {
420 for (nv =
NVAR; nv--; ){
421 dU[
k][
j][
i][nv] = state.
rhs[*in][nv];
422 UH[
k][
j][
i][nv] = d->
Uc[
k][
j][
i][nv] + state.
rhs[*in][nv];
425 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++) {
426 for (nv =
NVAR; nv--; ){
427 dU[
k][
j][
i][nv] += state.
rhs[*in][nv];
428 UH[
k][
j][
i][nv] += state.
rhs[*in][nv];
456 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++) {
457 VAR_LOOP(nv) d->Vc[nv][k][j][i] = state.v[*in][nv];
459 #if THERMAL_CONDUCTION == EXPLICIT
460 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++){
465 #if (RESISTIVITY == EXPLICIT) && (defined STAGGERED_MHD)
478 #if (RESISTIVITY == EXPLICIT) && !(defined STAGGERED_MHD)
492 for ((*in) = indx.
beg-1; (*in) <= indx.
end+1; (*in)++){
493 for (nv =
NVAR; nv--; ){
494 state.
up[*in][nv] += dU[
k][
j][
i][nv];
495 state.
um[*in][nv] += dU[
k][
j][
i][nv];
497 flagp[*in] = flagm[*in] = d->
flag[
k][
j][
i];
506 for ((*in) = indx.
beg-1; (*in) <= indx.
end+1; (*in)++) {
507 NVAR_LOOP(nv) state.vh[*in][nv] = d->Vc[nv][k][j][i];
508 state.flag[*in] = d->flag[k][j][i];
512 for ((*in) = indx.
beg - 2; (*in) <= indx.
end + 1; (*in)++){
514 state.
bn[*in] = state.
uL[*in][
BXn];
543 Riemann (&state, indx.
beg - 1, indx.
end, Dts->
cmax, grid);
547 #if (PARABOLIC_FLUX & EXPLICIT)
566 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++) {
568 #if VISCOSITY == EXPLICIT
569 inv_dtp =
MAX(inv_dtp, dcoeff[*in][
MX1]);
571 #if RESISTIVITY == EXPLICIT
572 EXPAND(inv_dtp =
MAX(inv_dtp, dcoeff[*in][
BX1]); ,
573 inv_dtp =
MAX(inv_dtp, dcoeff[*in][
BX2]); ,
574 inv_dtp =
MAX(inv_dtp, dcoeff[*in][
BX3]);)
577 inv_dtp =
MAX(inv_dtp, dcoeff[*in][ENG]);
579 inv_dtp *= inv_dl[*in]*inv_dl[*in];
585 #if UPDATE_VECTOR_POTENTIAL == YES
594 for ((*in) = indx.
beg; (*in) <= indx.
end; (*in)++) {
595 NVAR_LOOP(nv) d->Uc[k][j][i][nv] += state.rhs[*in][nv];
void Boundary(const Data *d, int idim, Grid *grid)
void CT_StoreEMF(const State_1D *state, int beg, int end, Grid *grid)
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
void ResetState(const Data *, State_1D *, Grid *)
double ** vh
Primitive state at n+1/2 (only for one step method)
double **** J
Electric current defined as curl(B).
void GetCurrent(const Data *d, int dir, Grid *grid)
void print1(const char *fmt,...)
double ** rhs
Conservative right hand side.
void States(const State_1D *, int, int, Grid *)
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
double **** Vc
The main four-index data array used for cell-centered primitive variables.
static void CTU_CT_Source(double **, double **, double **, double *, int, int, Grid *)
double g_dt
The current integration time step.
double * GetInverse_dl(const Grid *)
void SB_CorrectFluxes(Data_Arr U, double t, double dt, Grid *grid)
#define DOM
Computational domain (interior)
void CT_Update(const Data *d, Data_Arr Bs, double dt, Grid *grid)
void FARGO_ShiftSolution(Data_Arr U, Data_Arr Us, Grid *grid)
double ** vR
Primitive variables to the right of the interface, .
#define ARRAY_3D(nx, ny, nz, type)
double * cmax
Maximum signal velocity for hyperbolic eqns.
void FARGO_AddVelocity(const Data *, Grid *)
unsigned char *** flag
Pointer to a 3D array setting useful integration flags that are retrieved during integration.
void ConsToPrim3D(Data_Arr U, Data_Arr V, unsigned char ***flag, RBox *box)
void MakeState(State_1D *)
double inv_dtp
Inverse of diffusion (parabolic) time step .
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
int CheckNaN(double **, int, int, int)
void FARGO_SubtractVelocity(const Data *, Grid *)
#define THERMAL_CONDUCTION
#define TOT_LOOP(k, j, i)
int g_dir
Specifies the current sweep or direction of integration.
int ConsToPrim(double **ucons, double **uprim, int ibeg, int iend, unsigned char *flag)
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
#define ARRAY_4D(nx, ny, nz, nv, type)
void SetIndexes(Index *indx, Grid *grid)
void RightHandSide(const State_1D *state, Time_Step *Dts, int beg, int end, double dt, Grid *grid)
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
double ** uR
same as vR, in conservative vars
#define ARRAY_1D(nx, type)
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
double g_time
The current integration time.
double * bn
Face magentic field, bn = bx(i+1/2)
double **** Uc
The main four-index data array used for cell-centered conservative variables.
void SB_SaveFluxes(State_1D *state, Grid *grid)
double ** um
same as vm, in conservative vars
void ParabolicFlux(Data_Arr V, Data_Arr J, double ***T, const State_1D *state, double **dcoeff, int beg, int end, Grid *grid)
double ** vL
Primitive variables to the left of the interface, .
double ** up
same as vp, in conservative vars
#define ARRAY_2D(nx, ny, type)
void FlagShock(const Data *, Grid *)
void VectorPotentialUpdate(const Data *d, const void *vp, const State_1D *state, const Grid *grid)
#define QUIT_PLUTO(e_code)
#define TRANSVERSE_LOOP(indx, ip, i, j, k)
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
double ** uL
same as vL, in conservative vars
void CT_AverageMagneticField(double ****bf, double ****UU, Grid *grid)