PLUTO  4.4-patch2
Functions
cooling.h File Reference

Cooling main header file. More...

Go to the source code of this file.

Functions

double CompEquil (double, double, double *)
 
void CoolingSource (const Data *, double, timeStep *, Grid *)
 
double GetMaxRate (double *, double *, double)
 
void Jacobian (double *, double *, double **)
 
void Numerical_Jacobian (double *, double **)
 
void PowerLawCooling (Data_Arr, double, timeStep *, Grid *)
 
void Radiat (double *, double *)
 
double SolveODE_CK45 (double *, double *, double *, double, double, intList *)
 
double SolveODE_RKF23 (double *, double *, double *, double, intList *)
 
double SolveODE_RKF12 (double *, double *, double *, double, intList *)
 
double SolveODE_RK4 (double *, double *, double *, double, intList *)
 
void H2RateTables (double, double *)
 

Detailed Description

Contains function prototypes and basic definitions for the cooling modules.

Author
A. Mignone (migno.nosp@m.ne@t.nosp@m.o.inf.nosp@m.n.it)
Date
Apr 08, 2019

Function Documentation

◆ CompEquil()

double CompEquil ( double  N,
double  T,
double *  v0 
)
Parameters
[in]nthe particle number density (not needed, but kept for compatibility)
[in]Tthe temperature (in K) for which equilibrium must be found.
[in,out]van array of primitive variables. On input, only density needs to be defined. On output, fractions will be updated to the equilibrium values.

Compute the equilibrium ionization balance for (rho,T)

◆ CoolingSource()

void CoolingSource ( const Data d,
double  dt,
timeStep Dts,
Grid grid 
)

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]gridpointer to an array of Grid structures

◆ GetMaxRate()

double GetMaxRate ( double *  v0,
double *  k1,
double  T0 
)

Return an estimate of the maximum rate (dimension 1/time) in the chemical network. This will serve as a "stiffness" detector in the main ode integrator.

For integration to be carried explicitly all the time, return a small value (1.e-12).

◆ H2RateTables()

void H2RateTables ( double  T,
double *  krvals 
)

On first call, compute and store rate coefficient arrays needed for the H2_COOL module. On subsequent calls, use linear interpolation between adjacent tabulated values to obtain the coefficients at the desired temperature T.

◆ Jacobian()

void Jacobian ( double *  v,
double *  rhs,
double **  dfdy 
)

Compute the jacobian J(k,l) = dfdy

k = row index l = column index

J(0,0) J(0,1) ... J(0, n-1) J(1,0) J(1,1) .... J(1, n-1) . . . . . . . . . J(n-1,0) .... J(n-1, n-1)

or,

+--------------------—+

  • | |
  • | |
  • | |
  • dX'/dX | dX'/dp |
  • (JXX) | (JXp) |
  • | |
  • | | +-----------—+-----—+
  • dp'/dX | dp'/dp |
  • (JpX) | Jpp | +--------------------—+

◆ Numerical_Jacobian()

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.

◆ PowerLawCooling()

void PowerLawCooling ( Data_Arr  VV,
double  dt,
timeStep Dts,
Grid grid 
)
Parameters
[in,out]VVa pointer to the PLUTO 3D data array containing pimitive variables.
[in]dtthe current integration time step
[in]Dtsa pointer to the timeStep structure
[in]gridpointer to an array of Grid structures

◆ Radiat()

void Radiat ( double *  v,
double *  rhs 
)

Cooling for optically thin plasma up to about 200,000 K Plasma composition: H, HeI-II, CI-V, NI-V, OI-V, NeI-V, SI-V Assumed abundances in elem_ab Uses S : Array = Variables vector x line points rhs : output for the system of ODE ibeg, iend : begin and end points of the current line

Cooling for neutral or singly ionized gas: good up to about 35,000 K in equilibrium or shocks in neutral gas up to about 80 km/s. Assumed abundances in ab Uses t : Kelvin dene : electron density cm*-3 fneut : hydrogen neutral fraction (adimensionale) ci,cr : H ionization and recombination rate coefficients

em(1) = TOTAL EMISSIVITY : (ergs cm**3 s**-1) em(2) = Ly alpha + two photon continuum: Aggarwal MNRAS 202, 10**4.3 K em(3) = H alpha: Aggarwal, Case B em(4) = He I 584 + two photon + 623 (all n=2 excitations): Berrington &Kingston,JPB 20:

em(5) = C I 9850 + 9823: Mendoza, IAU 103, 5000 K em(6) = C II, 156 micron: Mendoza, 10,000 K em(7) = C II] 2325 A: Mendoza, 15,000 K em(8) = N I 5200 A: Mendoza, 7500 K em(9) = N II 6584 + 6548 A: Mendoza em(10) = O I 63 micron: Mendoza,2500 K em(11) = O I 6300 A + 6363 A: Mendoza, 7500 K em(12) = O II 3727: Mendoza em(13) = Mg II 2800: Mendoza em(14) = Si II 35 micron: Dufton&Kingston, MNRAS 248 em(15) = S II 6717+6727: Mendoza em(16) = Fe II 25 micron: Nussbaumer&Storey em(17) = Fe II 1.6 micron em(18) = thermal energy lost by ionization em(19) = thermal energy lost by recombination (2/3 kT per recombination. The ionization energy lost is not included here.

Provide r.h.s. for tabulated cooling.

◆ SolveODE_CK45()

double SolveODE_CK45 ( double *  v0,
double *  k1,
double *  v5th,
double  dt,
double  tol,
intList vars 
)

Explicit 5th order adaptive step size Cash-Karp integrator. Attempt to integrate the cooling network with a single step if the error is within the given tolerance. Otherwise, reduce the time step and repeat the integration as many times as necessary.

Parameters
[in]v_ininitial condition array
[in]k_inright hand side, R(v_in)
[out]v5thfinal solution array
[in]dtthe time step we want to achieve
[in]tolthe relative tolerance
[in]*varsthe list of dependent variables
Returns
The suggested time step for the next iteration.

◆ SolveODE_RK4()

double SolveODE_RK4 ( double *  v0,
double *  k1,
double *  v4th,
double  dt,
intList var_list 
)

Solve the system of ODE with a standard RK4 integrator (no adaptive step).

◆ SolveODE_RKF12()

double SolveODE_RKF12 ( double *  v0,
double *  k1,
double *  v2nd,
double  dt,
intList vars 
)

Explicit 2nd order integrator with 1st-order embedded solution. Attempt to integrate the cooling network with a single step if the error is within the given tolerance. The error is computed from the difference between 1st and 2nd order solutions.

Parameters
[in]v0initial condition array
[in]k1right hand side, R(v0)
[out]v2ndfinal solution array
[in]dtthe time step we want to achieve
[in]*varsthe list of dependent variables
Returns
The (scaled) error between 1st and 2nd-order solutions

◆ SolveODE_RKF23()

double SolveODE_RKF23 ( double *  v0,
double *  k1,
double *  v3rd,
double  dt,
intList vars 
)

Explicit 3rd order integrator with 2nd-order embedded solution. Attempt to integrate the cooling network with a single step if the error is within the given tolerance.

\[ y^{\rm 3rd} = y_0 + dt*(k1 + k2 + 4*k3)/6 k1 = RHS(y(0)) k2 = RHS(y(0) + dt*k1) k3 = RHS(y(0) + 0.25*dt*(k1 + k2)) y(2nd) = y(0) + dt*(k1 + k2)/2 \]

The error is computed from the difference between 1st and 2nd order solutions.

Parameters
[in]v0initial condition array
[in]k1right hand side, R(v0)
[out]v2ndfinal solution array
[in]dtthe time step we want to achieve
[in]*varsthe list of dependent variables
Returns
The (scaled) error between 2nd and 3rd-order solutions