PLUTO
Macros | Functions
sts.c File Reference

Super Time Stepping driver for integration of diffusion terms. More...

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

Go to the source code of this file.

Macros

#define STS_NU   0.01
 Sets the nu parameter used in STS. More...
 
#define STS_MAX_STEPS   1024
 

Functions

static void STS_ComputeSubSteps (double, double tau[], int)
 
static double STS_FindRoot (double, double, double)
 
static double STS_CorrectTimeStep (int, double)
 
void STS (const Data *d, Time_Step *Dts, Grid *grid)
 

Detailed Description

Super Time Stepping driver for integration of diffusion terms.

Take one step in the solution of the diffusion equation $ dU/dt = R(\partial^2U) $ where R is a nonlinear right hand side involving second derivatives. The super step is taken to be equal to the current time step g_dt and the number of substeps Nsts is given by solving the following nonlinear equation:

\[ \frac{\Delta t}{\Delta t_{\rm par}} = \frac{N}{2\sqrt{\nu}}\frac{(1+\sqrt{\nu})^{2N} - (1-\sqrt{\nu})^{2N}} {(1+\sqrt{\nu})^{2N} + (1-\sqrt{\nu})^{2N}} \]

where $\nu$ is set by the macro STS_NU (default equal to 0.01). The previous relation is given by Eq. [2.10] of Alexiades et al. with the explicit parabolic time step $\Delta t_{\rm par}$ being computed from

\[ \frac{2}{N_d} \max\left[ \frac{\eta_x}{\Delta x^2} + \frac{\eta_y}{\Delta y^2} + \frac{\eta_z}{\Delta z^2} \right] \Delta t_{\rm par} = C_p < \frac{1}{N_d} \]

where $C_p$ is the parabolic Courant number, $ N_d $ is the number of spatial dimensions and the maximum of the square bracket is computed during the call to ParabolicRHS.

This function is called in an operator-split way before/after advection has been carried out.

References

  • Alexiades, V., Amiez, A., & Gremaud E.-A. 1996, Com. Num. Meth. Eng., 12, 31
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
T. Matsakos
Date
Aug 27, 2015

Definition in file sts.c.

Macro Definition Documentation

#define STS_MAX_STEPS   1024

Definition at line 49 of file sts.c.

#define STS_NU   0.01

Sets the nu parameter used in STS.

Definition at line 46 of file sts.c.

Function Documentation

void STS ( const Data *  d,
Time_Step *  Dts,
Grid *  grid 
)

Solve diffusion equation(s) using Super-Time-Stepping.

Parameters
[in,out]dpointer to Data structure
[in,out]Dtspointer to Time_Step structure
[in]gridpointer to an array of Grid structures

Definition at line 55 of file sts.c.

64 {
65  int i, j, k, nv, n, m;
66  double N, ts[STS_MAX_STEPS];
67  double dt_par, tau, tsave, inv_dtp;
68  static Data_Arr rhs;
69  RBox *box = GetRBox(DOM, CENTER);
70 
71  if (rhs == NULL) rhs = ARRAY_4D(NX3_TOT, NX2_TOT, NX1_TOT, NVAR, double);
72 
73 /* -------------------------------------------------------------
74  Compute the conservative vector in order to start the cycle.
75  This step will be useless if the data structure
76  contains Vc as well as Uc (for future improvements).
77  ------------------------------------------------------------- */
78 
79  PrimToCons3D(d->Vc, d->Uc, box);
80  tsave = g_time;
81 
82 /* ------------------------------------------------------------
83  Main STS Loop starts here
84  ------------------------------------------------------------ */
85 
86  m = 0;
87  n = STS_MAX_STEPS;
88  while (m < n){
89 
90  g_intStage = m + 1;
91  Boundary(d, ALL_DIR, grid);
92  inv_dtp = ParabolicRHS(d, rhs, 1.0, grid);
93 
94  /* --------------------------------------------------------------
95  At the first step (m=0) compute (explicit) parabolic time
96  step. Restriction on explicit time step should be
97 
98  [2*eta_x/dx^2 + 2*eta_y/dy^2 + 2*eta_z/dz^2]*dt < 1
99 
100  which in PLUTO is normally written as
101 
102  [2/dtp_x + 2/dtp_y + 2/dtp_z]*dt/Ndim = 1/Ndim = cfl_par
103 
104  where Ndim is the number of spatial dimensions.
105  -------------------------------------------------------------- */
106 
107  if (m == 0){
108  Dts->inv_dtp = inv_dtp/(double)DIMENSIONS;
109  #ifdef PARALLEL
110  MPI_Allreduce (&Dts->inv_dtp, &tau, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
111  Dts->inv_dtp = tau;
112  #endif
113  Dts->inv_dtp = MAX(Dts->inv_dtp, 1.e-18);
114  dt_par = Dts->cfl_par/(2.0*Dts->inv_dtp); /* explicit parabolic
115  time step */
116  /* -------------------------------------------
117  Compute the number of steps needed to fit
118  the supertimestep with the advection one.
119  ------------------------------------------- */
120 
121  N = STS_FindRoot(1.0, dt_par, g_dt);
122  N = floor(N+1.0);
123  n = (int)N;
124 
125  Dts->Nsts = n;
126  if (n > STS_MAX_STEPS){
127  print1 ("! STS: the number of substeps (%d) is > %d\n",n, STS_MAX_STEPS);
128  QUIT_PLUTO(1);
129  }
130 
131  /* ------------------------------------------
132  Adjust explicit timestep to make Nsts an
133  integer number and compute the substeps.
134  ------------------------------------------ */
135 
136  if (Dts->Nsts > 1){
137  dt_par = STS_CorrectTimeStep(Dts->Nsts, g_dt);
138  STS_ComputeSubSteps(dt_par, ts, Dts->Nsts);
139  }
140 
141  if (Dts->Nsts == 1) ts[0] = g_dt;
142 
143  } /* end if m == 0 */
144 
145  tau = ts[n-m-1];
146 
147  /* -----------------------------------------------
148  Update cell-centered conservative variables
149  ----------------------------------------------- */
150 
151  DOM_LOOP (k,j,i){
152  #if VISCOSITY == SUPER_TIME_STEPPING
153  EXPAND(d->Uc[k][j][i][MX1] += tau*rhs[k][j][i][MX1]; ,
154  d->Uc[k][j][i][MX2] += tau*rhs[k][j][i][MX2]; ,
155  d->Uc[k][j][i][MX3] += tau*rhs[k][j][i][MX3];)
156  #endif
157  #if (RESISTIVITY == SUPER_TIME_STEPPING)
158  EXPAND(d->Uc[k][j][i][BX1] += tau*rhs[k][j][i][BX1]; ,
159  d->Uc[k][j][i][BX2] += tau*rhs[k][j][i][BX2]; ,
160  d->Uc[k][j][i][BX3] += tau*rhs[k][j][i][BX3];)
161  #endif
162  #if HAVE_ENERGY
163  #if (THERMAL_CONDUCTION == SUPER_TIME_STEPPING) || \
164  (RESISTIVITY == SUPER_TIME_STEPPING) || \
165  (VISCOSITY == SUPER_TIME_STEPPING)
166  d->Uc[k][j][i][ENG] += tau*rhs[k][j][i][ENG];
167  #endif
168  #endif
169  }
170  #if (defined STAGGERED_MHD) && (RESISTIVITY == SUPER_TIME_STEPPING)
171 
172  /* -----------------------------------------------
173  Update staggered magnetic field variables
174  ----------------------------------------------- */
175 
176 /* This piece of code should be used for simple Cartesian 2D
177  configurations with constant resistivity.
178  It shows that the stability limit of the resistive part of the
179  induction equation is larger if the right hand side is written
180  as Laplacian(B) rather than curl(J)
181 
182 double ***Bx, ***By, ***Bz, Jz, dx, dy, dt;
183 static double ***Ez, ***Rx, ***Ry;
184 if (Ez == NULL) {
185  Ez = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
186  Rx = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
187  Ry = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
188 }
189 Bx = d->Vs[BX1s];
190 By = d->Vs[BX2s];
191 dx = grid[IDIR].dx[IBEG];
192 dy = grid[JDIR].dx[JBEG];
193 
194 dt = tau;
195 k=0;
196 for (j = 1; j < NX2_TOT-1; j++) for (i = 1; i < NX1_TOT-1; i++){
197  Jz = (By[k][j][i+1] - By[k][j][i])/dx - (Bx[k][j+1][i] - Bx[k][j][i])/dy;
198  Ez[k][j][i] = Jz*g_inputParam[ETAZ];
199  Rx[k][j][i] = (Bx[k][j][i+1] - 2.0*Bx[k][j][i] + Bx[k][j][i-1])/(dx*dx)
200  + (Bx[k][j+1][i] - 2.0*Bx[k][j][i] + Bx[k][j-1][i])/(dy*dy);
201 
202  Ry[k][j][i] = (By[k][j][i+1] - 2.0*By[k][j][i] + By[k][j][i-1])/(dx*dx)
203  + (By[k][j+1][i] - 2.0*By[k][j][i] + By[k][j-1][i])/(dy*dy);
204 }
205 
206 JDOM_LOOP(j) for (i = IBEG-1; i <= IEND; i++){
207 // Bx[k][j][i] += - dt/dy*(Ez[k][j][i] - Ez[k][j-1][i]);
208  Bx[k][j][i] += dt*g_inputParam[ETAZ]*Rx[k][j][i];
209 }
210 
211 for (j = JBEG-1; j <= JEND; j++) IDOM_LOOP(i){
212 // By[k][j][i] += dt/dx*(Ez[k][j][i] - Ez[k][j][i-1]);
213  By[k][j][i] += dt*g_inputParam[ETAZ]*Ry[k][j][i];
214 }
215 */
216  CT_Update (d, d->Vs, tau, grid);
217  CT_AverageMagneticField (d->Vs, d->Uc, grid);
218  #endif
219 
220  /* --------------------------------------------------
221  Unflag zone tagged with the ENTROPY_SWITCH since
222  only total energy can be evolved using STS
223  -------------------------------------------------- */
224 
225  #if ENTROPY_SWITCH
226  TOT_LOOP(k,j,i) d->flag[k][j][i] &= ~FLAG_ENTROPY;
227  #endif
228 
229  /* ----------------------------------------------
230  convert conservative variables to primitive
231  for next iteration. Increment loop index.
232  ---------------------------------------------- */
233 
234  ConsToPrim3D(d->Uc, d->Vc, d->flag, box);
235  g_time += ts[n-m-1];
236  m++;
237  }
238  g_time = tsave; /* restore initial time step */
239 }
MX3
#define MX3
Definition: mod_defs.h:22
Boundary
void Boundary(const Data *d, int idim, Grid *grid)
Definition: boundary.c:36
MAX
#define MAX(a, b)
Definition: macros.h:101
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
MX1
#define MX1
Definition: mod_defs.h:20
Data_Arr
double **** Data_Arr
Definition: pluto.h:492
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
n
static int n
Definition: analysis.c:3
GetRBox
RBox * GetRBox(int, int)
Definition: rbox.c:232
CENTER
#define CENTER
Definition: pluto.h:200
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
PrimToCons3D
void PrimToCons3D(Data_Arr V, Data_Arr U, RBox *box)
Definition: mappers3D.c:86
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
g_intStage
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
Definition: globals.h:98
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
g_dt
double g_dt
The current integration time step.
Definition: globals.h:118
ParabolicRHS
double ParabolicRHS(const Data *d, Data_Arr dU, double dt, Grid *grid)
Definition: parabolic_rhs.c:50
DOM
#define DOM
Computational domain (interior)
Definition: pluto.h:152
CT_Update
void CT_Update(const Data *d, Data_Arr Bs, double dt, Grid *grid)
Definition: ct.c:29
STS_ComputeSubSteps
static void STS_ComputeSubSteps(double, double tau[], int)
Definition: sts.c:242
DATA::flag
unsigned char *** flag
Pointer to a 3D array setting useful integration flags that are retrieved during integration.
Definition: structs.h:55
ConsToPrim3D
void ConsToPrim3D(Data_Arr U, Data_Arr V, unsigned char ***flag, RBox *box)
Definition: mappers3D.c:16
TIME_STEP::inv_dtp
double inv_dtp
Inverse of diffusion (parabolic) time step .
Definition: structs.h:217
STS_CorrectTimeStep
static double STS_CorrectTimeStep(int, double)
Definition: sts.c:289
TOT_LOOP
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
FLAG_ENTROPY
#define FLAG_ENTROPY
Update pressure using entropy equation.
Definition: pluto.h:182
TIME_STEP::cfl_par
double cfl_par
Courant number for diffusion (STS only).
Definition: structs.h:221
NX3_TOT
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
j
int j
Definition: analysis.c:2
MX2
#define MX2
Definition: mod_defs.h:21
ARRAY_4D
#define ARRAY_4D(nx, ny, nz, nv, type)
Definition: prototypes.h:173
k
int k
Definition: analysis.c:2
BX3
#define BX3
Definition: mod_defs.h:27
STS_FindRoot
static double STS_FindRoot(double, double, double)
Definition: sts.c:258
STS_MAX_STEPS
#define STS_MAX_STEPS
Definition: sts.c:49
i
int i
Definition: analysis.c:2
g_time
double g_time
The current integration time.
Definition: globals.h:117
DATA::Uc
double **** Uc
The main four-index data array used for cell-centered conservative variables.
Definition: structs.h:37
BX1
#define BX1
Definition: mod_defs.h:25
TIME_STEP::Nsts
int Nsts
Maximum number of substeps used in STS.
Definition: structs.h:223
BX2
#define BX2
Definition: mod_defs.h:26
NVAR
#define NVAR
Definition: pluto.h:609
RBOX
Definition: structs.h:346
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2
CT_AverageMagneticField
void CT_AverageMagneticField(double ****bf, double ****UU, Grid *grid)
Definition: ct_field_average.c:13
ALL_DIR
#define ALL_DIR
Definition: pluto.h:196

Here is the call graph for this function:

Here is the caller graph for this function:

void STS_ComputeSubSteps ( double  dtex,
double  tau[],
int  ssorder 
)
static

Definition at line 242 of file sts.c.

246 {
247  int i;
248  double S=0.0;
249 
250  for (i = 0; i < ssorder; i++) {
251  tau[i] = dtex / ((-1.0 + STS_NU)*cos(((2.0*i+1.0)*CONST_PI)/(2.0*ssorder))
252  + 1.0 + STS_NU);
253  S += tau[i];
254  }
255 }
STS_NU
#define STS_NU
Sets the nu parameter used in STS.
Definition: sts.c:46
i
int i
Definition: analysis.c:2
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269

Here is the caller graph for this function:

double STS_CorrectTimeStep ( int  n0,
double  dta 
)
static

Definition at line 289 of file sts.c.

293 {
294  double a,b,c;
295  double dtr;
296 
297  a = (1.0-sqrt(STS_NU))/(1.0+sqrt(STS_NU));
298  b = pow(a,2.0*n0);
299  c = (1.0-b)/(1.0+b);
300 
301  dtr = dta*2.0*sqrt(STS_NU)/(n0*c);
302  return(dtr);
303 }
a
static double a
Definition: init.c:135
STS_NU
#define STS_NU
Sets the nu parameter used in STS.
Definition: sts.c:46
menu.c
tuple c
Definition: menu.py:375

Here is the caller graph for this function:

double STS_FindRoot ( double  x0,
double  dtr,
double  dta 
)
static

Definition at line 258 of file sts.c.

262 {
263  double a,b,c;
264  double Ns, Ns1;
265  int n;
266 
267  n = 0;
268 
269  Ns = x0+1.0;
270  Ns1 = x0;
271 
272  while(fabs(Ns-Ns1) >= 1.0e-5){
273  Ns = Ns1;
274  a = (1.0-sqrt(STS_NU))/(1.0+sqrt(STS_NU));
275  b = pow(a,2.0*Ns);
276  c = (1.0-b)/(1.0+b);
277  Ns1 = Ns + (dta - dtr*Ns/(2.0*sqrt(STS_NU))*c)
278  /(dtr/(2.0*sqrt(STS_NU))*(c-2.0*Ns*b*log(a)*(1.0+c)/(1.0+b)));
279  n += 1;
280  if (n == 128){
281  print1 ("! STS_FindRoot: max number of iterations exceeded");
282  QUIT_PLUTO(1);
283  }
284  }
285  return(Ns);
286 }
a
static double a
Definition: init.c:135
n
static int n
Definition: analysis.c:3
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
STS_NU
#define STS_NU
Sets the nu parameter used in STS.
Definition: sts.c:46
menu.c
tuple c
Definition: menu.py:375
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function: