PLUTO
Functions
cooling_source.c File Reference

Integrate cooling and reaction source terms. More...

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

Go to the source code of this file.

Functions

void CoolingSource (const Data *d, double dt, Time_Step *Dts, Grid *GXYZ)
 
void Numerical_Jacobian (double *v, double **J)
 

Detailed Description

Integrate cooling and reaction source terms.

Solve the system of ordinary differential equations describing optically thin cooling and ionization network (if any). We use an adaptive step size integration method and follow the strategy outlined in section 5.1 of Tesileanu et al. (2008) for handling stiff cells.

On output, a time-step estimate for the next time level is computed using the relative or absolute variation obtained during the integration of the ODE (from t(n) –> t(n+1))

\[ \Delta t_c = \min_{ijk}\left[\frac{\Delta t^n M_R}{\epsilon}\right] \,\quad\rm{where}\quad \epsilon = \max\left(\left|\frac{p^{n+1}}{p^n} - 1\right|,\, |X^{n+1}-X^n|\right) \]

where $ M_R $ is the maximum cooling rate (defined by the global variable g_maxCoolingRate) and X are the chemical species.

References

  • "Simulating radiative astrophysical flows with the PLUTO code: a non-equilibrium, multi-species cooling function"
    Tesileanu, Mignone and Massaglia, A&A (2008) 488, 429
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
O. Tesileanu B. Vaidya
Date
April 22, 2014

Definition in file cooling_source.c.

Function Documentation

void CoolingSource ( const Data *  d,
double  dt,
Time_Step *  Dts,
Grid *  GXYZ 
)

Integrate cooling and reaction source terms.

Parameters
[in,out]dpointer to Data structure
[in]dtthe time step to be taken
[out]Dtspointer to the Time_Step structure
[in]GXYZpointer to an array of Grid structures

Definition at line 38 of file cooling_source.c.

48 {
49  int nv, k, j, i, stiff, status;
50  double err, scrh, min_tol = 2.e-5;
51  double mu0, T0, T1, mu1, prs;
52  double v0[NVAR], v1[NVAR], k1[NVAR];
53  double maxrate;
54  intList var_list;
55 
56 /* --------------------------------------------------------
57  Set number and indices of the time-dependent variables
58  -------------------------------------------------------- */
59 
60  var_list.nvar = NIONS+1;
61  var_list.indx[0] = PRS;
62  for (nv = 0; nv < NIONS; nv++) var_list.indx[nv+1] = NFLX+nv;
63 
64 /* -----------------------------------------------------
65  Zero k1
66  ----------------------------------------------------- */
67 
68  NVAR_LOOP(nv) k1[nv] = 0.0;
69 
70 /* -----------------------------------------------------------
71  Begin Integration
72  ----------------------------------------------------------- */
73 
74  DOM_LOOP(k,j,i){ /* -- span the computational domain -- */
75 
76  /* --------------------------------------------------
77  Skip integration if cell has been tagged with
78  FLAG_INTERNAL_BOUNDARY or FLAG_SPLIT_CELL
79  (only for AMR)
80  -------------------------------------------------- */
81 
82  #if INTERNAL_BOUNDARY == YES
83  if (d->flag[k][j][i] & FLAG_INTERNAL_BOUNDARY) continue;
84  #endif
85  if (d->flag[k][j][i] & FLAG_SPLIT_CELL) continue;
86 
87  /* ----------------------------------------------
88  Compute temperature and internal energy from
89  density, pressure and concentrations.
90  ---------------------------------------------- */
91 
92  NVAR_LOOP(nv) v0[nv] = v1[nv] = d->Vc[nv][k][j][i];
93  prs = v0[PRS];
94  mu0 = MeanMolecularWeight(v0);
95  T0 = v0[PRS]/v0[RHO]*KELVIN*mu0;
96  #if EOS == IDEAL
97  v0[RHOE] = v1[RHOE] = prs/(g_gamma-1.0);
98  #else
99  v1[RHOE] = v0[RHOE] = InternalEnergy(v0, T0);
100  #endif
101 
102  if (T0 <= 0.0){
103  print ("! CoolingSource: negative initial temperature\n");
104  print (" %12.6e %12.6e\n",v0[RHOE], v0[RHO]);
105  print (" at: %f %f\n",GXYZ[IDIR].x[i], GXYZ[JDIR].x[j]);
106  QUIT_PLUTO(1);
107  }
108 
109  /* -------------------------------------------
110  Get estimated time step based on
111  the max rate of the reaction network.
112  ------------------------------------------- */
113 
114  Radiat(v0, k1);
115 
116  maxrate = GetMaxRate (v0, k1, T0);
117  stiff = (dt*maxrate > 1.0 ? 1:0);
118 
119  /* ---------------------------------------------------
120  If the system is not stiff, then try to advance
121  with an explicit 2-nd order midpoint rule
122  --------------------------------------------------- */
123 
124  if (!stiff){ /* -- no stiffness: try explicit -- */
125 
126  scrh = SolveODE_RKF12 (v0, k1, v1, dt, min_tol, &var_list);
127 /* scrh = SolveODE_RKF23 (v0, k1, v1, dt, min_tol); */
128 
129 /* -- error is too big ? --> use some other integrator -- */
130 
131  if (scrh < 0.0) SolveODE_CK45 (v0, k1, v1, dt, min_tol, &var_list);
132 
133  } else { /* -- if stiff = 1 use more sophisticated integrators -- */
134 
135 /* SolveODE_ROS34 (v0, k1, v1, dtsub, min_tol); */
136 
137  int nsub, k;
138  double dtsub, dtnew, t;
139 
140  nsub = ceil(dt*maxrate);
141  dtsub = dt/(double)nsub;
142 
143  /* ----------------------------------
144  use sub-time stepping
145  ---------------------------------- */
146 
147  t = 0.0;
148  for (k = 1; 1; k++){
149  dtnew = SolveODE_CK45 (v0, k1, v1, dtsub, min_tol, &var_list);
150  /* dtnew = SolveODE_ROS34(v0, k1, v1, dtsub, min_tol); */
151  t += dtsub;
152  /* printf ("%d / %d %12.6e %12.6e, t/dt = %f\n",
153  k,nsub,dtsub,dtnew, t/dt);*/
154 
155  if (fabs(t/dt - 1.0) < 1.e-9) break;
156 
157  v0[RHOE] = v1[RHOE];
158  NIONS_LOOP(nv) v0[nv] = v1[nv];
159 
160  Radiat(v0, k1);
161  dtsub = MIN (dtnew, dt - t);
162  }
163 
164  if (k > 100) print ("! CoolingSource: Number of substeps exceeded 100 (%d)\n",k);
165  if (fabs(t/dt - 1.0) > 1.e-12) {
166  print ("! CoolingSource: dt mismatch\n");
167  QUIT_PLUTO(1);
168  }
169  } /* -- end if (stiff) -- */
170 
171  /* -- Constrain ions to lie between [0,1] -- */
172 
173  NIONS_LOOP(nv){
174  v1[nv] = MAX(v1[nv], 0.0);
175  v1[nv] = MIN(v1[nv], 1.0);
176  }
177  #if COOLING == H2_COOL
178  v1[X_H2] = MIN(v1[X_H2], 0.5);
179  #endif
180 
181  /* -- pressure must be positive -- */
182 
183  mu1 = MeanMolecularWeight(v1);
184  #if EOS == IDEAL
185  prs = v1[RHOE]*(g_gamma - 1.0);
186  T1 = prs/v1[RHO]*KELVIN*mu1;
187  #elif EOS == PVTE_LAW
188  status = GetEV_Temperature(v1[RHOE], v1, &T1);
189  prs = v1[RHO]*T1/(KELVIN*mu1);
190  #endif
191 
192  if (prs < 0.0) prs = g_smallPressure;
193 
194  /* -- Check final temperature -- */
195 
196 
197  if (T1 < g_minCoolingTemp && T0 > g_minCoolingTemp)
198  prs = g_minCoolingTemp*v1[RHO]/(KELVIN*mu1);
199 
200  /* ------------------------------------------
201  Suggest next time step based on
202  fractional variaton.
203  ------------------------------------------ */
204 
205  err = fabs(prs/d->Vc[PRS][k][j][i] - 1.0);
206 
207  #if COOLING == MINEq
208  for (nv = NFLX; nv < NFLX + NIONS - Fe_IONS; nv++)
209  #else
210  NIONS_LOOP(nv)
211  #endif
212  err = MAX(err, fabs(d->Vc[nv][k][j][i] - v1[nv]));
213 
214  scrh = dt*g_maxCoolingRate/err;
215 
216  Dts->dt_cool = MIN(Dts->dt_cool, scrh);
217 
218  /* ---- Update solution array ---- */
219 
220  d->Vc[PRS][k][j][i] = prs;
221  NIONS_LOOP(nv) d->Vc[nv][k][j][i] = v1[nv];
222 
223  } /* -- end loop on points -- */
224 }
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
MAX
#define MAX(a, b)
Definition: macros.h:101
NIONS_LOOP
#define NIONS_LOOP(n)
Definition: pluto.h:614
g_gamma
double g_gamma
Definition: globals.h:112
GetMaxRate
double GetMaxRate(double *, double *, double)
Definition: maxrate.c:4
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
InternalEnergy
double InternalEnergy(double *v, double T)
Definition: internal_energy.c:255
n
static int n
Definition: analysis.c:3
Radiat
void Radiat(double *, double *)
Definition: radiat.c:94
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
CoolingSource
void CoolingSource(const Data *d, double dt, Time_Step *Dts, Grid *GXYZ)
Definition: cooling_source.c:38
KELVIN
#define KELVIN
Definition: pluto.h:401
DATA::flag
unsigned char *** flag
Pointer to a 3D array setting useful integration flags that are retrieved during integration.
Definition: structs.h:55
INT_LIST::nvar
int nvar
Number of variables.
Definition: structs.h:330
FLAG_SPLIT_CELL
#define FLAG_SPLIT_CELL
Zone is covered by a finer level (AMR only).
Definition: pluto.h:183
SolveODE_CK45
double SolveODE_CK45(double *v0, double *k1, double *v5th, double dt, double tol, intList *vars)
Definition: cooling_ode_solver.c:7
MIN
#define MIN(a, b)
Definition: macros.h:104
SolveODE_RKF12
double SolveODE_RKF12(double *v0, double *k1, double *v2nd, double dt, double tol, intList *vars)
Definition: cooling_ode_solver.c:226
TIME_STEP::dt_cool
double dt_cool
Cooling time step.
Definition: structs.h:219
NFLX
#define NFLX
Definition: mod_defs.h:32
IDIR
#define IDIR
Definition: pluto.h:193
GetEV_Temperature
int GetEV_Temperature(double rhoe, double *v, double *T)
Definition: internal_energy.c:290
NVAR_LOOP
#define NVAR_LOOP(n)
Definition: pluto.h:618
NIONS
#define NIONS
Definition: cooling.h:28
j
int j
Definition: analysis.c:2
if
if(divB==NULL)
Definition: analysis.c:10
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
g_minCoolingTemp
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
Definition: globals.h:106
INT_LIST
Definition: structs.h:328
i
int i
Definition: analysis.c:2
INT_LIST::indx
int indx[2046]
Array of integers containg variables indices.
Definition: structs.h:329
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
Fe_IONS
#define Fe_IONS
Definition: cooling.h:18
X_H2
#define X_H2
Definition: cooling.h:30
JDIR
#define JDIR
Definition: pluto.h:194
g_maxCoolingRate
double g_maxCoolingRate
The maximum fractional variation due to cooling from one step to the next.
Definition: globals.h:104
NVAR
#define NVAR
Definition: pluto.h:609
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:

void Numerical_Jacobian ( double *  v,
double **  J 
)

Compute Jacobian matrix numerically and get eigenvector and eigenvalue decomposition

The purpose of this function is to detect stiffness.

Note: This function is EXTREMELY time consuming.

Definition at line 227 of file cooling_source.c.

238 {
239  int n, nv, k, l;
240  double eps;
241  double vpp[NVAR], vp[NVAR], vm[NVAR], vmm[NVAR];
242  double Rpp[NVAR], Rp[NVAR], Rm[NVAR], Rmm[NVAR];
243 
244  eps = 1.e-5;
245  n = NIONS + 1;
246 
247 /* -- partial derivs with respect to pressure -- */
248 
249  for (nv = 0; nv < NVAR; nv++){
250  vp[nv] = vm[nv] = v[nv];
251  }
252  vp[PRS] *= 1.0 + eps;
253  vm[PRS] *= 1.0 - eps;
254 
255  Radiat (vp, Rp);
256  Radiat (vm, Rm);
257 
258  for (k = 0; k < n - 1; k++){
259  J[k][n - 1] = (Rp[k + NFLX] - Rm[k + NFLX])/(2.0*v[PRS]*eps);
260  }
261  J[n - 1][n - 1] = (Rp[PRS] - Rm[PRS])/(2.0*v[PRS]*eps);
262 
263 /* -- partial deriv with respect to ions -- */
264 
265  for (l = 0; l < n - 1; l++){
266 
267  for (nv = 0; nv < NVAR; nv++){
268  vp[nv] = vm[nv] = v[nv];
269  }
270  vp [l + NFLX] = v[l + NFLX] + eps;
271  vm [l + NFLX] = v[l + NFLX] - eps;
272 
273  vp [l + NFLX] = MIN(vp [l + NFLX], 1.0);
274  vm [l + NFLX] = MAX(vm [l + NFLX], 0.0);
275 
276  Radiat (vp , Rp);
277  Radiat (vm , Rm);
278 
279  for (k = 0; k < n - 1; k++){
280  J[k][l] = (Rp[k + NFLX] - Rm[k + NFLX])/(vp[l + NFLX] - vm[l + NFLX]);
281  }
282  J[n - 1][l] = (Rp[PRS] - Rm[PRS])/(vp[l + NFLX] - vm[l + NFLX]);
283  }
284 }
MAX
#define MAX(a, b)
Definition: macros.h:101
n
static int n
Definition: analysis.c:3
Radiat
void Radiat(double *, double *)
Definition: radiat.c:94
MIN
#define MIN(a, b)
Definition: macros.h:104
eps
#define eps
NFLX
#define NFLX
Definition: mod_defs.h:32
NIONS
#define NIONS
Definition: cooling.h:28
k
int k
Definition: analysis.c:2
NVAR
#define NVAR
Definition: pluto.h:609

Here is the call graph for this function: