PLUTO
Functions
rkc.c File Reference

Runge-Kutta-Chebyshev driver for integration of diffusion terms. More...

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

Go to the source code of this file.

Functions

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

Detailed Description

Runge-Kutta-Chebyshev driver for integration of diffusion terms.

Take one step in the solution of the diffusion equation $ dQ/dt = F$ using Runge-Kutta-Chebyshev (RKC) method.

\[ \begin{array}{lcl} Y_0 & = & Q^n \\ \noalign{\medskip} Y_1 & = & Y_0 + \tilde{\mu}_1\tau F_0 \\ \noalign{\medskip} Y_j & = & (1 - \mu_j - \nu_j) Y_0 + \mu_j Y_{j-1} + \nu_j Y_{j-2} + \tilde{\mu}_j \tau F_{j-1} + \tilde{\gamma}_j\tau F_0 \quad \rm{for} \quad j = 2,...,s \\ \noalign{\medskip} Q^{n+1} & = & Y_s \end{array} \]

where

\[\tau = \Delta t,\, \tilde{\mu}_1 = b_1 w_1,\, \mu_j = 2 b_j \frac{w_0}{b_{j-1}}, \, \nu_j = -\frac{b_j}{b_{j-2}},\, \tilde{\mu}_j = 2 b_j \frac{w_1}{b_{j-1}},\, \tilde{\gamma}_j = -\alpha_{j-1} \tilde{\mu}_j,\, \]

\[ \alpha_j = 1 - b_j T_j(w_0),\, w_0 = 1 + \epsilon/s^2,\, w_1 = T'_s(w_0)/T''_s(w_0),\, b_j = T''_j(w_0)/(T'_j(w_0))^2,\, b0=b1=b2,\, \epsilon = 2/13 \]

Here the T's are the Chebyshev polynomial of the first kind:

\[ T_j(x) = \cos(j {\rm arccos(x)}) ,\, T_0(x) = 1,\, T_1(x)=x,\, T_j(x)=2*x*T_jm1(x) - T_jm2(x), \quad \rm{for}\quad j = 2,...,s \]

while sprad=Spectral_Radius-> 4,8,12 kappa/dx/dx (check)

Only the parabolic terms of the equations are advanced by a single step g_dt equal to the current integration time step and the number of RKC steps (Nrkc) is computed... The explicit parabolic time step is computed in the same manner as in STS.

References

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

Definition in file rkc.c.

Function Documentation

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

Solve diffusion equation using Runge-Kutta-Chebyshev (RKC) method

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 58 of file rkc.c.

67 {
68  int i, j, k, nv, s, s_RKC = 0;
69  int nv_indx, var_list[NVAR], nvar_rkc;
70  double kappa_dl2, eta_dl2, nu_dl2;
71  double absh, sprad, Y;
72  double epsilon = 2./13.;
73  double w_0, w_1;
74  double scrh, scr1, scr2;
75  double mu_j, nu_j, mu_tilde,a_jm1;
76  double Tj, dTj, d2Tj, Tjm1, Tjm2, dTjm1, dTjm2, d2Tjm1, d2Tjm2;
77  double b_j, b_jm1, b_jm2, arg;
78  double dt_par;
79  static Data_Arr Y_jm1, Y_jm2, UU_0, F_0, F_jm1;
80  static double **v;
81  static unsigned char *flag;
82 
83  if (Y_jm1 == NULL) {
84  Y_jm1 = ARRAY_4D(NX3_TOT, NX2_TOT, NX1_TOT, NVAR, double);
85  Y_jm2 = ARRAY_4D(NX3_TOT, NX2_TOT, NX1_TOT, NVAR, double);
86  UU_0 = ARRAY_4D(NX3_TOT, NX2_TOT, NX1_TOT, NVAR, double);
87  F_0 = ARRAY_4D(NX3_TOT, NX2_TOT, NX1_TOT, NVAR, double);
88  F_jm1 = ARRAY_4D(NX3_TOT, NX2_TOT, NX1_TOT, NVAR, double);
89  flag = ARRAY_1D(NMAX_POINT, unsigned char);
90  v = ARRAY_2D(NMAX_POINT, NVAR, double);
91  }
92 
93 /* -------------------------------------------------------
94  set the number of variables to be evolved with RKC
95  ------------------------------------------------------- */
96 
97  i = 0;
98  for (nv = 0; nv < NVAR; nv++) var_list[nv] = 0;
99  #if VISCOSITY == RK_CHEBYSHEV
100  EXPAND(var_list[i++] = MX1; ,
101  var_list[i++] = MX2; ,
102  var_list[i++] = MX3;)
103  #endif
104 
105  #if RESISTIVITY == RK_CHEBYSHEV
106  EXPAND(var_list[i++] = BX1; ,
107  var_list[i++] = BX2; ,
108  var_list[i++] = BX3;)
109  #endif
110  #if (THERMAL_CONDUCTION == RK_CHEBYSHEV) \
111  || (VISCOSITY == RK_CHEBYSHEV && EOS == IDEAL) \
112  || (RESISTIVITY == RK_CHEBYSHEV && EOS == IDEAL)
113 
114  var_list[i++] = ENG;
115  #endif
116  nvar_rkc = i;
117 
118 /* -------------------------------------------------
119  get conservative vector UU
120  ------------------------------------------------- */
121 
122  KDOM_LOOP(k){
123  JDOM_LOOP(j){
124  IDOM_LOOP(i){
125  for (nv = NVAR; nv--; ) {
126  v[i][nv] = d->Vc[nv][k][j][i];
127  }}
128  PrimToCons (v, UU_0[k][j], IBEG, IEND);
129  }}
130 
131 /* -------------------------------------------------
132  Compute the parabolic time step by calling
133  the RHS function once outside the main loop.
134  ------------------------------------------------- */
135 
136  g_intStage = 1;
137  Boundary(d, ALL_DIR, grid);
138  #if SHOCK_FLATTENING == MULTID
139  FindShock (d, grid);
140  #endif
141 
142  Dts->inv_dtp = ParabolicRHS(d, F_0, 1.0, grid);
143  Dts->inv_dtp /= (double) DIMENSIONS;
144  #ifdef PARALLEL
145  MPI_Allreduce (&Dts->inv_dtp, &scrh, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
146  Dts->inv_dtp = scrh;
147  #endif
148  dt_par = Dts->cfl_par/(2.0*Dts->inv_dtp); /* -- explicit parabolic time step -- */
149 
150 /* -----------------------------------------------------
151  Compute spectral radius (placeholder for now...)
152  -----------------------------------------------------*/
153 
154  D_EXPAND( sprad = 4.*Dts->inv_dtp;,
155  sprad = 8.*Dts->inv_dtp;,
156  sprad = 12.*Dts->inv_dtp;)
157 
158 D_EXPAND( sprad = 2./dt_par; ,
159  sprad = 4./dt_par; ,
160  sprad = 6./dt_par;)
161 
162 /* -----------------------------------------------------
163  Compute number of RKC steps from sprad
164  -----------------------------------------------------*/
165 
166  absh = g_dt;
167  s_RKC = 1 + (int)(sqrt(1.54*absh*sprad + 1));
168  Dts->Nrkc = s_RKC;
169 /* limit dt */
170  if(absh > 0.653*s_RKC*s_RKC/sprad){
171  print1 ("! RKC: outside parameter range\n");
172  QUIT_PLUTO(1);
173  }
174 
175 /* ------------------------------------------------------------------
176  RKC main stage loop
177  ------------------------------------------------------------------ */
178 
179 /* -- define coefficients -- */
180 
181  w_0 = 1. + epsilon/(1.*s_RKC)/(1.*s_RKC);
182  scr1 = w_0*w_0 - 1.;
183  scr2 = sqrt(scr1);
184  arg = s_RKC*log(w_0 + scr2);
185  w_1 = sinh(arg)*scr1/(cosh(arg)*s_RKC*scr2 - w_0*sinh(arg));
186  b_jm1 = 1./(2.*w_0)/(2.*w_0);
187  b_jm2 = b_jm1;
188  mu_tilde = w_1*b_jm1;
189 
190 /* -------------------------------------------------------------
191  First RKC stage here:
192  - take one step in conservative variables,
193  - convert to primitive
194  ------------------------------------------------------------- */
195 
196  DOM_LOOP (k,j,i){
197  for (nv = 0; nv < NVAR; nv++) {
198  Y_jm1[k][j][i][nv] = Y_jm2[k][j][i][nv] = UU_0[k][j][i][nv];
199  }
200  for (nv_indx = 0; nv_indx < nvar_rkc; nv_indx++){
201  nv = var_list[nv_indx];
202  Y_jm1[k][j][i][nv] = Y_jm2[k][j][i][nv] + mu_tilde*absh*F_0[k][j][i][nv];
203  }
204  }
205 
206  KDOM_LOOP(k){
207  JDOM_LOOP(j){
208  ConsToPrim (Y_jm1[k][j], v, IBEG, IEND, flag);
209  IDOM_LOOP(i){
210  for (nv = NVAR; nv--; ){
211  d->Vc[nv][k][j][i] = v[i][nv];
212  }}
213  }}
214 
215 /* -----------------------------------------------------------------
216  Loop over remaining RKC stages s = 2,..,s_RKC steps
217  ----------------------------------------------------------------- */
218 
219 /* -- Initialize Chebyshev polynomials (recursive) -- */
220 
221  Tjm1 = w_0; Tjm2 = 1.0; dTjm1 = 1.0;
222  dTjm2 = 0.0; d2Tjm1 = 0.0; d2Tjm2 = 0.0;
223 
224  for (s = 2; s <= s_RKC; s++){
225 
226  /* ---- compute coefficients ---- */
227 
228  Tj = 2.*w_0*Tjm1 - Tjm2;
229  dTj = 2.*w_0*dTjm1 - dTjm2 + 2.*Tjm1;
230  d2Tj = 2.*w_0*d2Tjm1 - d2Tjm2 + 4.*dTjm1;
231  b_j = d2Tj/dTj/dTj;
232  a_jm1 = 1. - Tjm1*b_jm1;
233  mu_j = 2.*w_0*b_j/b_jm1;
234  nu_j = - b_j/b_jm2;
235  mu_tilde = mu_j*w_1/w_0;
236 
237  /* -- call again boundary and take a new step -- */
238 
239  g_intStage = s;
240  Boundary(d, ALL_DIR, grid);
241  ParabolicRHS (d, F_jm1, 1.0, grid);
242  DOM_LOOP (k,j,i){
243  for (nv_indx = 0; nv_indx < nvar_rkc; nv_indx++){
244  nv = var_list[nv_indx];
245  Y = mu_j*Y_jm1[k][j][i][nv] + nu_j*Y_jm2[k][j][i][nv];
246  Y_jm2[k][j][i][nv] = Y_jm1[k][j][i][nv];
247  Y_jm1[k][j][i][nv] = Y + (1. - mu_j - nu_j)*UU_0[k][j][i][nv]
248  + absh*mu_tilde*(F_jm1[k][j][i][nv]
249  - a_jm1*F_0[k][j][i][nv]);
250  }
251  } /* END DOM_LOOP */
252 
253  /* -- Put parameters of s -> s-1 (outside of dom loop) -- */
254 
255  b_jm2 = b_jm1; b_jm1 = b_j; Tjm2 = Tjm1;
256  Tjm1 = Tj; dTjm2 = dTjm1; dTjm1 = dTj;
257  d2Tjm2 = d2Tjm1; d2Tjm1 = d2Tj;
258 
259  /* ----------------------------------------------------------
260  convert conservative variables to primitive for the next
261  s steps (or timestep if s == s_RKC)
262  ---------------------------------------------------------- */
263 
264  KDOM_LOOP(k){
265  JDOM_LOOP(j){
266  ConsToPrim (Y_jm1[k][j], v, IBEG, IEND, flag);
267  IDOM_LOOP(i){
268  for (nv = NVAR; nv--; ){
269  d->Vc[nv][k][j][i] = v[i][nv];
270  }}
271  }}
272  }/* s loop */
273 }/* func */
IDEAL
#define IDEAL
Definition: pluto.h:45
MX3
#define MX3
Definition: mod_defs.h:22
Boundary
void Boundary(const Data *d, int idim, Grid *grid)
Definition: boundary.c:36
EOS
#define EOS
Definition: pluto.h:341
KDOM_LOOP
#define KDOM_LOOP(k)
Definition: macros.h:36
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
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
configure.scrh
tuple scrh
Definition: configure.py:200
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
VISCOSITY
#define VISCOSITY
Definition: pluto.h:385
Y
#define Y
Definition: radiat.modified_by_Bhargav.c:4
FindShock
void FindShock(const Data *, Grid *)
TIME_STEP::inv_dtp
double inv_dtp
Inverse of diffusion (parabolic) time step .
Definition: structs.h:217
JDOM_LOOP
#define JDOM_LOOP(j)
Definition: macros.h:35
THERMAL_CONDUCTION
#define THERMAL_CONDUCTION
Definition: pluto.h:365
TIME_STEP::cfl_par
double cfl_par
Courant number for diffusion (STS only).
Definition: structs.h:221
ConsToPrim
int ConsToPrim(double **ucons, double **uprim, int ibeg, int iend, unsigned char *flag)
Definition: mappers.c:89
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
if
if(divB==NULL)
Definition: analysis.c:10
ARRAY_4D
#define ARRAY_4D(nx, ny, nz, nv, type)
Definition: prototypes.h:173
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
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
D_EXPAND
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;)
Definition: analysis.c:27
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
RESISTIVITY
#define RESISTIVITY
Definition: pluto.h:357
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
BX2
#define BX2
Definition: mod_defs.h:26
RK_CHEBYSHEV
#define RK_CHEBYSHEV
Definition: pluto.h:69
NVAR
#define NVAR
Definition: pluto.h:609
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
IDOM_LOOP
#define IDOM_LOOP(i)
Definition: macros.h:34
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2
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: