PLUTO
Macros | Functions
parabolic_rhs.c File Reference

Compute right hand side for diffusion terms. More...

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

Go to the source code of this file.

Macros

#define ADD_RESISTIVITY
 
#define ADD_VISCOSITY
 
#define ADD_TC
 

Functions

double ParabolicRHS (const Data *d, Data_Arr dU, double dt, Grid *grid)
 

Detailed Description

Compute right hand side for diffusion terms.

The ParabolicRHS function computes the right hand side, in conservative or divergence form, of the parabolic (diffusion) operators only:

\[ \pd{q}{t} = - \nabla\cdot\Big(D\nabla q\Big) \]

Here q is a generic cell-centered quantity (like momentum, magnetic field or total energy) and $\vec{F} = D\nabla q$ is the corresponding flux computed elsewhere. Contributions may simultaneously come from viscosity, magnetic resistivity and thermal conduction for which the flux functions are computed in the function ViscousFlux(), ResistiveFlux() or TC_Flux().

This function is intended for operator split algorithms (STS, RKC) in order to solve for the parabolic part of the equations only.

This function also computes the inverse diffusion time step by adding, for each diffusion equation, contributions coming from different directions:

\[ \Delta t_{p}^{-1} = \max\left[\frac{\eta_{i+\HALF} + \eta_{i-\HALF}}{2\Delta x^2} + \frac{\eta_{j+\HALF} + \eta_{j-\HALF}}{2\Delta y^2} + \frac{\eta_{k+\HALF} + \eta_{k-\HALF}}{2\Delta z^2}\right] \]

where $\eta_{x,y,z}$ are the diffusion coefficients available at cell interfaces in the three directions, respectively, and the maximum is taken over the local processor grid.

Date
Sept 1, 2014
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
P. Tzeferacos

Definition in file parabolic_rhs.c.

Macro Definition Documentation

#define ADD_RESISTIVITY
Value:
(RESISTIVITY == SUPER_TIME_STEPPING || \
RESISTIVITY == RK_CHEBYSHEV)
SUPER_TIME_STEPPING
#define SUPER_TIME_STEPPING
Definition: pluto.h:68
RESISTIVITY
#define RESISTIVITY
Definition: pluto.h:357
RK_CHEBYSHEV
#define RK_CHEBYSHEV
Definition: pluto.h:69

Definition at line 40 of file parabolic_rhs.c.

#define ADD_TC
Value:
(THERMAL_CONDUCTION == SUPER_TIME_STEPPING || \
THERMAL_CONDUCTION == RK_CHEBYSHEV)
SUPER_TIME_STEPPING
#define SUPER_TIME_STEPPING
Definition: pluto.h:68
THERMAL_CONDUCTION
#define THERMAL_CONDUCTION
Definition: pluto.h:365
RK_CHEBYSHEV
#define RK_CHEBYSHEV
Definition: pluto.h:69

Definition at line 46 of file parabolic_rhs.c.

#define ADD_VISCOSITY
Value:
(VISCOSITY == SUPER_TIME_STEPPING || \
VISCOSITY == RK_CHEBYSHEV)
VISCOSITY
#define VISCOSITY
Definition: pluto.h:385
SUPER_TIME_STEPPING
#define SUPER_TIME_STEPPING
Definition: pluto.h:68
RK_CHEBYSHEV
#define RK_CHEBYSHEV
Definition: pluto.h:69

Definition at line 43 of file parabolic_rhs.c.

Function Documentation

double ParabolicRHS ( const Data *  d,
Data_Arr  dU,
double  dt,
Grid *  grid 
)
Parameters
[in]V3D array containing primitive variables
[out]dU3D array containing the conservative right hand sides
[in]dtthe time step
[in]grida pointer to an array of Grid structure
Returns
On output it returns the maximum diffusion coefficients among all dissipative term over the local processor grid.

Definition at line 50 of file parabolic_rhs.c.

60 {
61  int i, j, k, nv;
62  double r, th, r_1, s_1;
63  double dtdl, dtdV, dt_rdr;
64  double *inv_dl, inv_dl2;
65  double *rp, *A, *du, **vh, *C, **tc_flx;
66  double dvp[NVAR], dvm[NVAR], dvl;
67  double ****V;
68  static double **res_flx;
69  static double **dcoeff, ***T;
70  static double **ViF, **ViS;
71 
72  static Data_Arr C_dtp;
73  static State_1D state;
74  double max_inv_dtp[3],inv_dtp;
75 
76  i = j = k = 0;
77 
78  if (C_dtp == NULL) {
79  C_dtp = ARRAY_4D(NX3_TOT, NX2_TOT, NX1_TOT, NVAR, double);
80  #if ADD_TC
81  MakeState (&state);
82  #ifdef CH_SPACEDIM
83  i = j = k = 1;
84  D_EXPAND(i = NMAX_POINT;, j = NMAX_POINT;, k = NMAX_POINT;)
85  T = ARRAY_3D(k, j, i, double);
86  #else
87  T = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
88  #endif
89  #endif
90 
91  ViF = ARRAY_2D(NMAX_POINT, NVAR, double);
92  ViS = ARRAY_2D(NMAX_POINT, NVAR, double);
93  res_flx = ARRAY_2D(NMAX_POINT, NVAR, double);
94  dcoeff = ARRAY_2D(NMAX_POINT, NVAR, double);
95  }
96 
97  V = d->Vc; /* -- set a pointer to the primitive vars array -- */
98 
99  #if ADD_RESISTIVITY && (defined STAGGERED_MHD)
100  GetCurrent (d, -1, grid);
101  #endif
102 
103 /* ------------------------------------------------
104  compute the temperature array if TC is needed
105  ------------------------------------------------ */
106 
107  #if ADD_TC
108  TOT_LOOP(k,j,i) T[k][j][i] = V[PRS][k][j][i]/V[RHO][k][j][i];
109  vh = state.vh;
110  tc_flx = state.tc_flux;
111  #endif
112 
113 max_inv_dtp[0] = max_inv_dtp[1] = max_inv_dtp[2] = 0.0;
114 
115 /* ----------------------------------------------------------------------
116  L O O P O N X1 D I R E C T I O N ( I D I R )
117  ---------------------------------------------------------------------- */
118 
119  g_dir = IDIR;
120  A = grid[g_dir].A;
121  rp = grid[g_dir].xr;
122  #if ADD_RESISTIVITY && !(defined STAGGERED_MHD)
123  GetCurrent(d, g_dir, grid);
124  #endif
125  KDOM_LOOP (k){
126  JDOM_LOOP (j){
127 
128  g_j = j; g_k = k;
129 
130  /* -------------------------------------------------------------------
131  Compute viscous, resistive and thermal cond. fluxes in the x1 dir
132  ------------------------------------------------------------------- */
133 
134  #if ADD_VISCOSITY
135  ViscousFlux (V, ViF, ViS, dcoeff, IBEG-1, IEND, grid);
136  #endif
137 
138  #if ADD_RESISTIVITY
139  ResistiveFlux (V, d->J, res_flx, dcoeff, IBEG - 1, IEND, grid);
140  #endif
141 
142  #if ADD_TC
143  ITOT_LOOP (i) for (nv = NVAR; nv--; ) vh[i][nv] = V[nv][k][j][i];
144 
145  for (nv = NVAR; nv--; ) dvm[nv] = vh[1][nv] - vh[0][nv];
146  for(i=1; i<NX1_TOT-1; i++){
147  for (nv = NVAR; nv--; ) {
148  dvp[nv] = vh[i+1][nv] - vh[i][nv];
149  dvl = VAN_LEER(dvp[nv], dvm[nv]);
150  state.vp[i][nv] = vh[i][nv] + 0.5*dvl;
151  state.vm[i][nv] = vh[i][nv] - 0.5*dvl;
152  dvm[nv] = dvp[nv];
153  }}
154  TC_Flux (T, &state, dcoeff, IBEG - 1, IEND, grid);
155  #endif
156 
157  /* ---------------------------
158  compute inverse time-step
159  --------------------------- */
160 
161  if (g_intStage == 1){
162  inv_dl = GetInverse_dl(grid);
163  IDOM_LOOP(i){
164  C = C_dtp[k][j][i];
165  inv_dl2 = 0.5*inv_dl[i]*inv_dl[i];
166  #if ADD_VISCOSITY
167  C[MX1] = (dcoeff[i-1][MX1] + dcoeff[i][MX1])*inv_dl2;
168  #endif
169  #if ADD_RESISTIVITY
170  EXPAND(C[BX1] = (dcoeff[i-1][BX1] + dcoeff[i][BX1])*inv_dl2; ,
171  C[BX2] = (dcoeff[i-1][BX2] + dcoeff[i][BX2])*inv_dl2; ,
172  C[BX3] = (dcoeff[i-1][BX3] + dcoeff[i][BX3])*inv_dl2;)
173  #endif
174  #if ADD_TC
175  C[ENG] = (dcoeff[i-1][ENG] + dcoeff[i][ENG])*inv_dl2;
176  inv_dtp = dcoeff[i][ENG]*inv_dl[i]*inv_dl[i];
177  max_inv_dtp[g_dir] = MAX(max_inv_dtp[g_dir], inv_dtp);
178  #endif
179  }
180  }
181 
182  /* ---------------------------
183  Main X1-sweep
184  --------------------------- */
185 
186  IDOM_LOOP (i){
187 
188  du = dU[k][j][i];
189  r = grid[IDIR].x[i];
190  dtdV = dt/grid[IDIR].dV[i];
191  dtdl = dt/grid[IDIR].dx[i];
192 
193  #if HAVE_ENERGY
194  du[ENG] = 0.0; /* energy contribution may come from any of the
195  dissipative terms. It's better to initialize it
196  to zero and then add contribution separately. */
197  #endif
198 
199  /* -------------------------------------------
200  VISCOSITY: build rhs in the x1 direction
201  ------------------------------------------- */
202 
203  #if ADD_VISCOSITY
204  r_1 = 1.0/grid[IDIR].x[i];
205  #if GEOMETRY == CARTESIAN
206  EXPAND(du[MX1] = (ViF[i][MX1] - ViF[i-1][MX1])*dtdl + dt*ViS[i][MX1];,
207  du[MX2] = (ViF[i][MX2] - ViF[i-1][MX2])*dtdl + dt*ViS[i][MX2];,
208  du[MX3] = (ViF[i][MX3] - ViF[i-1][MX3])*dtdl + dt*ViS[i][MX3];)
209  #if HAVE_ENERGY
210  du[ENG] += (ViF[i][ENG] - ViF[i - 1][ENG])*dtdl;
211  #endif
212  #elif GEOMETRY == CYLINDRICAL
213  EXPAND(du[MX1] = (A[i]*ViF[i][MX1] - A[i-1]*ViF[i-1][MX1])*dtdV
214  + dt*ViS[i][MX1]; ,
215  du[MX2] = (A[i]*ViF[i][MX2] - A[i-1]*ViF[i-1][MX2])*dtdV
216  + dt*ViS[i][MX2]; ,
217  du[MX3] = (A[i]*A[i]*ViF[i][MX3] - A[i-1]*A[i-1]*ViF[i-1][MX3])*r_1*dtdV;)
218  #if HAVE_ENERGY
219  du[ENG] += (A[i]*ViF[i][ENG] - A[i-1]*ViF[i-1][ENG])*dtdV;
220  #endif
221  #elif GEOMETRY == POLAR
222  EXPAND(du[MX1] = (A[i]*ViF[i][MX1] - A[i-1]*ViF[i-1][MX1])*dtdV
223  + dt*ViS[i][MX1]; ,
224  du[MX2] = (A[i]*A[i]*ViF[i][MX2] - A[i-1]*A[i-1]*ViF[i-1][MX2])*r_1*dtdV; ,
225  du[MX3] = (A[i]*ViF[i][MX3] - A[i-1]*ViF[i-1][MX3])*dtdV
226  + dt*ViS[i][MX3];)
227  #if HAVE_ENERGY
228  du[ENG] += (A[i]*ViF[i][ENG] - A[i-1]*ViF[i-1][ENG])*dtdV;
229  #endif
230  #elif GEOMETRY == SPHERICAL
231  EXPAND(du[MX1] = (A[i]*ViF[i][MX1] - A[i-1]*ViF[i-1][MX1])*dtdV
232  + dt*ViS[i][MX1];,
233  du[MX2] = (A[i]*ViF[i][MX2] - A[i-1]*ViF[i-1][MX2])*dtdV
234  + dt*ViS[i][MX2];,
235  du[MX3] = (rp[i]*A[i]*ViF[i][MX3] - rp[i-1]*A[i-1]*ViF[i-1][MX3])*r_1*dtdV;)
236  #if HAVE_ENERGY
237  du[ENG] += (A[i]*ViF[i][ENG] - A[i-1]*ViF[i-1][ENG])*dtdV;
238  #endif
239  #endif
240  #endif /* -- VISCOSITY -- */
241 
242  /* ----------------------------------------------
243  RESISTIVITY: build rhs in the x1 direction
244  ---------------------------------------------- */
245 
246  #if ADD_RESISTIVITY
247  #if GEOMETRY == CARTESIAN /* -- x coordinate -- */
248  EXPAND(
249  du[BX1] = 0.0; ,
250  du[BX2] = -(res_flx[i][BX2] - res_flx[i-1][BX2])*dtdl; ,
251  du[BX3] = -(res_flx[i][BX3] - res_flx[i-1][BX3])*dtdl; )
252  #if EOS != ISOTHERMAL
253  du[ENG] += -(res_flx[i][ENG] - res_flx[i-1][ENG])*dtdl;
254  #endif
255 
256  #elif GEOMETRY == CYLINDRICAL /* -- r coordinate -- */
257  EXPAND(
258  du[BX1] = 0.0; ,
259  du[BX2] = -(A[i]*res_flx[i][BX2] - A[i-1]*res_flx[i-1][BX2])*dtdV; ,
260  du[BX3] = -(res_flx[i][BX3] - res_flx[i-1][BX3])*dtdl;)
261  #if EOS != ISOTHERMAL
262  du[ENG] += -(A[i]*res_flx[i][ENG] - A[i-1]*res_flx[i-1][ENG])*dtdV;
263  #endif
264 
265  #elif GEOMETRY == POLAR /* -- r coordinate -- */
266  EXPAND(
267  du[BX1] = 0.0; ,
268  du[BX2] = -(res_flx[i][BX2] - res_flx[i-1][BX2])*dtdl; ,
269  du[BX3] = -(A[i]*res_flx[i][BX3] - A[i-1]*res_flx[i-1][BX3])*dtdV; )
270  #if EOS != ISOTHERMAL
271  du[ENG] += -(A[i]*res_flx[i][ENG] - A[i-1]*res_flx[i-1][ENG])*dtdV;
272  #endif
273 
274  #elif GEOMETRY == SPHERICAL /* -- r coordinate -- */
275  dt_rdr = dtdl/r;
276  EXPAND(
277  du[BX1] = 0.0; ,
278  du[BX2] = -(rp[i]*res_flx[i][BX2] - rp[i-1]*res_flx[i-1][BX2])*dt_rdr; ,
279  du[BX3] = -(rp[i]*res_flx[i][BX3] - rp[i-1]*res_flx[i-1][BX3])*dt_rdr; )
280  #if EOS != ISOTHERMAL
281  du[ENG] += -(A[i]*res_flx[i][ENG] - A[i-1]*res_flx[i-1][ENG])*dtdV;
282  #endif
283  #endif
284  #endif /* -- RESISTIVITY -- */
285 
286  /* ---------------------------------------------------
287  THERMAL_CONDUCTION: build rhs in the x1 direction
288  --------------------------------------------------- */
289 
290  #if ADD_TC
291  #if GEOMETRY == CARTESIAN
292  du[ENG] += (tc_flx[i][ENG] - tc_flx[i-1][ENG])*dtdl;
293  #else
294  du[ENG] += (A[i]*tc_flx[i][ENG] - A[i-1]*tc_flx[i-1][ENG])*dtdV;
295  #endif
296  #endif
297  }
298  }}
299 
300 /* ----------------------------------------------------------------------
301  L O O P O N X2 D I R E C T I O N ( J D I R )
302  ---------------------------------------------------------------------- */
303 
304 #if DIMENSIONS > 1
305  g_dir = JDIR;
306  A = grid[g_dir].A;
307  #if ADD_RESISTIVITY && !(defined STAGGERED_MHD)
308  GetCurrent(d, g_dir, grid);
309  #endif
310  KDOM_LOOP (k){
311  IDOM_LOOP (i){
312 
313  g_i = i; g_k = k;
314 
315  /* -------------------------------------------------------------------
316  Compute viscous, resistive and thermal cond. fluxes in the x2 dir
317  ------------------------------------------------------------------- */
318 
319  #if ADD_VISCOSITY
320  ViscousFlux (V, ViF, ViS, dcoeff, JBEG-1, JEND, grid);
321  #endif
322  #if ADD_RESISTIVITY
323  ResistiveFlux (V, d->J, res_flx, dcoeff, JBEG - 1, JEND, grid);
324  #endif
325 
326  #if ADD_TC
327  JTOT_LOOP (j) for (nv = NVAR; nv--; ) vh[j][nv] = V[nv][k][j][i];
328 
329  for (nv = NVAR; nv--; ) dvm[nv] = vh[1][nv] - vh[0][nv];
330  for(j = 1; j < NX2_TOT-1; j++){
331  for (nv = NVAR; nv--; ) {
332  dvp[nv] = vh[j+1][nv] - vh[j][nv];
333  dvl = VAN_LEER(dvp[nv], dvm[nv]);
334  state.vp[j][nv] = vh[j][nv] + 0.5*dvl;
335  state.vm[j][nv] = vh[j][nv] - 0.5*dvl;
336  dvm[nv] = dvp[nv];
337  }}
338  TC_Flux (T, &state, dcoeff, JBEG - 1, JEND, grid);
339  #endif
340 
341  /* ---------------------------
342  compute inverse time-step
343  --------------------------- */
344 
345  if (g_intStage == 1){
346  inv_dl = GetInverse_dl(grid);
347  JDOM_LOOP(j){
348  C = C_dtp[k][j][i];
349  inv_dl2 = 0.5*inv_dl[j]*inv_dl[j];
350  #if ADD_VISCOSITY
351  C[MX1] += (dcoeff[j-1][MX1] + dcoeff[j][MX1])*inv_dl2;
352  #endif
353  #if ADD_RESISTIVITY
354  EXPAND(C[BX1] += (dcoeff[j-1][BX1] + dcoeff[j][BX1])*inv_dl2; ,
355  C[BX2] += (dcoeff[j-1][BX2] + dcoeff[j][BX2])*inv_dl2; ,
356  C[BX3] += (dcoeff[j-1][BX3] + dcoeff[j][BX3])*inv_dl2;)
357  #endif
358  #if ADD_TC
359  C[ENG] += (dcoeff[j-1][ENG] + dcoeff[j][ENG])*inv_dl2;
360 
361  inv_dtp = dcoeff[j][ENG]*inv_dl[j]*inv_dl[j];
362  max_inv_dtp[g_dir] = MAX(max_inv_dtp[g_dir], inv_dtp);
363  #endif
364  }
365  }
366 
367  /* ---------------------------
368  Main X2-sweep
369  --------------------------- */
370 
371  r = grid[IDIR].x[i];
372  JDOM_LOOP (j){
373 
374  du = dU[k][j][i];
375  dtdV = dt/grid[JDIR].dV[j];
376  dtdl = dt/grid[JDIR].dx[j];
377 
378  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
379  dtdl /= r;
380  dtdV /= r;
381  #endif
382 
383  /* ------------------------------------------
384  VISCOSITY: build rhs in the x2 direction
385  ------------------------------------------ */
386 
387  #if ADD_VISCOSITY
388  s_1 = 1.0/sin(grid[JDIR].x[j]);
389  #if GEOMETRY != SPHERICAL
390  EXPAND(du[MX1] += (ViF[j][MX1] - ViF[j-1][MX1])*dtdl + dt*ViS[j][MX1];,
391  du[MX2] += (ViF[j][MX2] - ViF[j-1][MX2])*dtdl + dt*ViS[j][MX2];,
392  du[MX3] += (ViF[j][MX3] - ViF[j-1][MX3])*dtdl + dt*ViS[j][MX3];)
393  #if HAVE_ENERGY
394  du[ENG] += (ViF[j][ENG] - ViF[j-1][ENG])*dtdl;
395  #endif
396  #elif GEOMETRY == SPHERICAL
397  EXPAND(du[MX1] += (A[j]*ViF[j][MX1] - A[j-1]*ViF[j-1][MX1])*dtdV
398  + dt*ViS[j][MX1];,
399  du[MX2] += (A[j]*ViF[j][MX2] - A[j-1]*ViF[j-1][MX2])*dtdV
400  + dt*ViS[j][MX2];,
401  du[MX3] += (A[j]*A[j]*ViF[j][MX3] - A[j-1]*A[j-1]*ViF[j-1][MX3])*s_1*dtdV;)
402  #if HAVE_ENERGY
403  du[ENG] += (A[j]*ViF[j][ENG] - A[j-1]*ViF[j-1][ENG])*dtdV;
404  #endif
405  #endif
406  #endif/* -- VISCOSITY -- */
407 
408  /* ----------------------------------------------
409  RESISTIVITY: build rhs in the x2 direction
410  ---------------------------------------------- */
411 
412  #if ADD_RESISTIVITY
413  #if GEOMETRY != SPHERICAL /* -- y coordinate -- */
414  EXPAND(
415  du[BX1] -= (res_flx[j][BX1] - res_flx[j-1][BX1])*dtdl; ,
416  ,
417  du[BX3] -= (res_flx[j][BX3] - res_flx[j-1][BX3])*dtdl; )
418  #if HAVE_ENERGY
419  du[ENG] -= (res_flx[j][ENG] - res_flx[j-1][ENG])*dtdl;
420  #endif
421  #elif GEOMETRY == SPHERICAL /* -- theta coordinate -- */
422  EXPAND(
423  du[BX1] -= (A[j]*res_flx[j][BX1] - A[j-1]*res_flx[j-1][BX1])*dtdV; ,
424  ,
425  du[BX3] -= (res_flx[j][BX3] - res_flx[j-1][BX3])*dtdl;)
426  #if HAVE_ENERGY
427  du[ENG] -= (A[j]*res_flx[j][ENG] - A[j-1]*res_flx[j-1][ENG])*dtdV;
428  #endif
429  #endif
430  #endif /* -- RESISTIVE MHD -- */
431 
432  /* ---------------------------------------------------
433  THERMAL_CONDUCTION: build rhs in the x2 direction
434  --------------------------------------------------- */
435 
436  #if ADD_TC
437  #if GEOMETRY != SPHERICAL
438  du[ENG] += (tc_flx[j][ENG] - tc_flx[j-1][ENG])*dtdl;
439  #else
440  du[ENG] += (A[j]*tc_flx[j][ENG] - A[j-1]*tc_flx[j-1][ENG])*dtdV;
441  #endif
442  #endif
443  }
444 
445  }}
446 #endif
447 
448 /* ----------------------------------------------------------------------
449  L O O P O N X3 D I R E C T I O N ( K D I R )
450  ---------------------------------------------------------------------- */
451 
452 #if DIMENSIONS == 3
453  g_dir = KDIR;
454  A = grid[g_dir].A;
455  #if ADD_RESISTIVITY && !(defined STAGGERED_MHD)
456  GetCurrent(d, g_dir, grid);
457  #endif
458  JDOM_LOOP (j){
459  IDOM_LOOP (i){
460 
461  g_i = i; g_j = j;
462 
463  /* -------------------------------------------------------------------
464  Compute viscous, resistive and thermal cond. fluxes in the x3 dir
465  ------------------------------------------------------------------- */
466 
467  #if ADD_VISCOSITY
468  ViscousFlux (V, ViF, ViS, dcoeff, KBEG-1, KEND, grid);
469  #endif
470  #if ADD_RESISTIVITY
471  ResistiveFlux (V, d->J, res_flx, dcoeff, KBEG - 1, KEND, grid);
472  #endif
473 
474  #if ADD_TC
475  KTOT_LOOP (k) for (nv = NVAR; nv--; ) vh[k][nv] = V[nv][k][j][i];
476 
477  for (nv = NVAR; nv--; ) dvm[nv] = vh[1][nv] - vh[0][nv];
478  for(k=1; k<NX3_TOT-1; k++){
479  for (nv = NVAR; nv--; ) {
480  dvp[nv] = vh[k+1][nv]- vh[k][nv];
481  dvl = VAN_LEER(dvp[nv], dvm[nv]);
482  state.vp[k][nv] = vh[k][nv] + 0.5*dvl;
483  state.vm[k][nv] = vh[k][nv] - 0.5*dvl;
484  dvm[nv] = dvp[nv];
485  }}
486  TC_Flux (T, &state, dcoeff, KBEG - 1, KEND, grid);
487  #endif
488 
489  /* ---------------------------
490  compute inverse time-step
491  --------------------------- */
492 
493  if (g_intStage == 1){
494  inv_dl = GetInverse_dl(grid);
495  KDOM_LOOP(k){
496  C = C_dtp[k][j][i];
497  inv_dl2 = 0.5*inv_dl[k]*inv_dl[k];
498  #if ADD_VISCOSITY
499  C[MX1] += (dcoeff[k-1][MX1] + dcoeff[k][MX1])*inv_dl2;
500  #endif
501  #if ADD_RESISTIVITY
502  EXPAND(C[BX1] += (dcoeff[k-1][BX1] + dcoeff[k][BX1])*inv_dl2; ,
503  C[BX2] += (dcoeff[k-1][BX2] + dcoeff[k][BX2])*inv_dl2; ,
504  C[BX3] += (dcoeff[k-1][BX3] + dcoeff[k][BX3])*inv_dl2;)
505  #endif
506  #if ADD_TC
507  C[ENG] += (dcoeff[k-1][ENG] + dcoeff[k][ENG])*inv_dl2;
508 
509  inv_dtp = dcoeff[k][ENG]*inv_dl[k]*inv_dl[k];
510  max_inv_dtp[g_dir] = MAX(max_inv_dtp[g_dir], inv_dtp);
511 
512  #endif
513  }
514  }
515 
516  /* ---------------------------
517  Main X3-sweep
518  --------------------------- */
519 
520  r = grid[IDIR].x[i];
521  th = grid[JDIR].x[j];
522  KDOM_LOOP (k){
523 
524  du = dU[k][j][i];
525  dtdV = dt/grid[KDIR].dV[k];
526  dtdl = dt/grid[KDIR].dx[k];
527  #if GEOMETRY == SPHERICAL
528  dtdl /= r*sin(th);
529  #endif
530 
531  /* ------------------------------------------
532  VISCOSITY: build rhs in the x3 direction
533  ------------------------------------------ */
534 
535  #if ADD_VISCOSITY
536  du[MX1] += (ViF[k][MX1] - ViF[k-1][MX1])*dtdl + dt*ViS[k][MX1];
537  du[MX2] += (ViF[k][MX2] - ViF[k-1][MX2])*dtdl + dt*ViS[k][MX2];
538  du[MX3] += (ViF[k][MX3] - ViF[k-1][MX3])*dtdl + dt*ViS[k][MX3];
539  #if HAVE_ENERGY
540  du[ENG] += (ViF[k][ENG] - ViF[k-1][ENG])*dtdl;
541  #endif
542  #endif/* -- VISCOSITY -- */
543 
544  /* ----------------------------------------------
545  RESISTIVITY: build rhs in the x3 direction
546  ---------------------------------------------- */
547 
548  #if ADD_RESISTIVITY
549  du[BX1] -= (res_flx[k][BX1] - res_flx[k-1][BX1])*dtdl;
550  du[BX2] -= (res_flx[k][BX2] - res_flx[k-1][BX2])*dtdl;
551  #if HAVE_ENERGY
552  du[ENG] -= (res_flx[k][ENG] - res_flx[k-1][ENG])*dtdl;
553  #endif
554  #endif /* -- RESISTIVITY -- */
555 
556  /* ---------------------------------------------------
557  THERMAL_CONDUCTION: build rhs in the x3 direction
558  --------------------------------------------------- */
559 
560  #if ADD_TC
561  du[ENG] += (tc_flx[k][ENG] - tc_flx[k-1][ENG])*dtdl;
562  #endif
563  }
564  }}
565 #endif
566 
567 /* OLD VERSION
568 D_EXPAND(th = max_inv_dtp[IDIR]; ,
569  th = MAX(th, max_inv_dtp[JDIR]); ,
570  th = MAX(th, max_inv_dtp[KDIR]);)
571 return th;
572 */
573 
574 /* --------------------------------------------------------------
575  take the maximum of inverse dt over domain and zero
576  right hand side in the internal boundary zones.
577  -------------------------------------------------------------- */
578 
579  th = 0.0;
580  DOM_LOOP(k,j,i){
581  #if ADD_VISCOSITY
582  th = MAX(th, C_dtp[k][j][i][MX1]);
583  #endif
584  #if ADD_RESISTIVITY
585  EXPAND(th = MAX(th, C_dtp[k][j][i][BX1]); ,
586  th = MAX(th, C_dtp[k][j][i][BX2]); ,
587  th = MAX(th, C_dtp[k][j][i][BX3]);)
588  #endif
589  #if ADD_TC
590  th = MAX(th, C_dtp[k][j][i][ENG]);
591  #endif
592  #if INTERNAL_BOUNDARY == YES
593  if (d->flag[k][j][i] & FLAG_INTERNAL_BOUNDARY) {
594  for (nv = NVAR; nv--; ) dU[k][j][i][nv] = 0.0;
595  }
596  #endif
597  }
598  return th;
599 }
FLAG_INTERNAL_BOUNDARY
#define FLAG_INTERNAL_BOUNDARY
Zone belongs to an internal boundary region and should be excluded from being updated in time...
Definition: pluto.h:184
MX3
#define MX3
Definition: mod_defs.h:22
EOS
#define EOS
Definition: pluto.h:341
MAX
#define MAX(a, b)
Definition: macros.h:101
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
KDOM_LOOP
#define KDOM_LOOP(k)
Definition: macros.h:36
ViscousFlux
void ViscousFlux(Data_Arr, double **, double **, double **, int, int, Grid *)
Definition: viscous_flux.c:14
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
DATA::J
double **** J
Electric current defined as curl(B).
Definition: structs.h:54
GetCurrent
void GetCurrent(const Data *d, int dir, Grid *grid)
Definition: res_functions.c:97
TC_Flux
void TC_Flux(double ***, const State_1D *, double **, int, int, Grid *)
Definition: tc_flux.c:57
GRID::xr
double * xr
Definition: structs.h:81
RHO
#define RHO
Definition: mod_defs.h:19
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
KTOT_LOOP
#define KTOT_LOOP(k)
Definition: macros.h:40
GetInverse_dl
double * GetInverse_dl(const Grid *)
Definition: set_geometry.c:205
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
VAN_LEER
#define VAN_LEER(a, b)
Definition: macros.h:113
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
GRID::dx
double * dx
Definition: structs.h:83
KDIR
#define KDIR
Definition: pluto.h:195
DATA::flag
unsigned char *** flag
Pointer to a 3D array setting useful integration flags that are retrieved during integration.
Definition: structs.h:55
MakeState
void MakeState(State_1D *)
Definition: tools.c:51
JDOM_LOOP
#define JDOM_LOOP(j)
Definition: macros.h:35
STATE_1D::vp
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
Definition: structs.h:142
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
TOT_LOOP
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
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
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
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
ISOTHERMAL
#define ISOTHERMAL
Definition: pluto.h:49
STATE_1D::vm
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
BX3
#define BX3
Definition: mod_defs.h:27
ITOT_LOOP
#define ITOT_LOOP(i)
Definition: macros.h:38
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
JTOT_LOOP
#define JTOT_LOOP(j)
Definition: macros.h:39
STATE_1D
Definition: structs.h:133
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
ADD_TC
#define ADD_TC
Definition: parabolic_rhs.c:46
KBEG
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
BX2
#define BX2
Definition: mod_defs.h:26
KEND
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:45
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
ResistiveFlux
void ResistiveFlux(Data_Arr V, Data_Arr curlB, double **res_flux, double **dcoeff, int beg, int end, Grid *grid)
Definition: res_flux.c:52
NVAR
#define NVAR
Definition: pluto.h:609
IDOM_LOOP
#define IDOM_LOOP(i)
Definition: macros.h:34
JEND
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:41
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
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87
HAVE_ENERGY
#define HAVE_ENERGY
Definition: pluto.h:395

Here is the call graph for this function:

Here is the caller graph for this function: