Thermodynamic relations for the caloric EoS, e=e(T,rho) and T=T(e,rho). More...

#include "pluto.h"

Include dependency graph for internal_energy.c:

Macros
#define	TV_ENERGY_TABLE_NX 1024

#define	TV_ENERGY_TABLE_NY 512

Functions
static double	rhoeFunc (double, void *)

static double	InternalEnergyDeriv (double rho, double T)

void	MakeInternalEnergyTable ()

void	MakeEV_TemperatureTable ()

double	InternalEnergy (double *v, double T)

int	GetEV_Temperature (double rhoe, double v, double T)

void	SoundSpeed2 (double *v, double cs2, double h, int beg, int end, int pos, Grid grid)

Variables
static Table2D	rhoe_tab
	A 2D table containing pre-computed values of internal energy stored at equally spaced node values of `Log(T)` and `Log(rho)` . More...

static Table2D	Trhoe_tab
	A 2D table containing pre-computed values of temperature by inverting e=e(T,rho) at equally spaced node values of `Log(rhoe)` and `Log(rho)`. More...

Detailed Description

Thermodynamic relations for the caloric EoS, e=e(T,rho) and T=T(e,rho).

Collect various functions involving computations between internal energy, density and temperature,

$e = e(T,\vec{X})\,,\qquad T = T(e,\vec{X})$

where e is the specific internal energy, T is the temperature and X are the ionization fractions. The first relation is used to retrieve e as a function of T and X by calling InternalEnergy() while the second one computes T from internal energy and fractions and it is handled by the function GetEV_Temperature():

If no chemistry is used (NIONS==0 ), ionization fractions are computed assuming equilibrium conditions and X=X(T,rho). This requires the numerical inversion of the second equation T=(e,X(T,rho)) which can be carried out using either a root finder algorithm (typically Brent's method) or, alternatively, lookup table together with cubic/linear (default) or bilinear interpolation.

For the root-finder version, a 2D table (Trhoe_tab) giving T=T(rhoe,rho) is pre-computed in the function MakeEV_TemperatureTable() and lookup table is used to narrow down the root interval using the values stored inside this table.

In the tabulated EOS version, a different table (rhoe_tab) rhoe=rhoe(T,rho) giving the internal energy at predefined values of temperature and density is pre-calculated in the function MakeInternalEnergyTable() so that internal energy and temperature can be found by a combination of lookup table and bilinear (direct or inverse) interpolation.

Author: A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
B. Vaidya

Date: 7 Jan, 2015

Definition in file internal_energy.c.

Macro Definition Documentation

#define TV_ENERGY_TABLE_NX 1024

Definition at line 47 of file internal_energy.c.

#define TV_ENERGY_TABLE_NY 512

Definition at line 50 of file internal_energy.c.

Function Documentation

int GetEV_Temperature	(	double	rhoe,
		double *	v,
		double *	T
	)

Compute temperature from internal energy, density and fractions by solving either

$\begin{array}{ll} \rho e - \rho e(T,\vec{X}) = 0 & \quad{(\rm with\; chemistry)} \\ \noalign{\medskip} \rho e - \rho e(T,\vec{X}(T,\rho)) = 0 & \quad{(\rm without\; chemistry)} \end{array}$

The previous equations are solved using a root finder algorithm (default is Brent or Ridder) or by lookup table / interpolation methods.

Parameters

[in]	rhoe	Internal energy of gas in code units.
[in]	v	1D array of primitive variables (only density and species need to be set).
[out]	T	Temperature in K

Returns: 0 on success, 1 on failure.

Definition at line 290 of file internal_energy.c.

 {
 #if TV_ENERGY_TABLE ==  YES
   int status;
   double rho;
   
   if (rhoe < 0.0) return 1;
 
   rho    = v[RHO];
   status = InverseLookupTable2D(&rhoe_tab, rho, rhoe, T);
   if (status != 0){
     WARNING(
       print ("! GetEV_Temperature(): table interpolation failure");
       print ("  [rho = %12.6e, rhoe = %12.6e]\n",rho,rhoe);
     )
     return 1;
   }  
 
   if (*T < T_CUT_RHOE) return 1;
   return 0;
 #else
 
   int    nv, status, i, j;
   double Tmin, Tmax,rho, lnx, lny;
   struct func_param param;
 
   if (rhoe < 0.0) return 1;
 
 /* ---------------------------------------------------
     Constrain species to lie between [0,1]
    --------------------------------------------------- */
  
 #if NIONS > 0
   NIONS_LOOP(nv) {
     v[nv] = MAX(v[nv],0.0);
     v[nv] = MIN(v[nv],1.0);
   }
   #if COOLING == H2_COOL
     v[X_H2] = MIN(v[X_H2], 0.5);
   #endif
 #endif
  
   param.v[RHO] = v[RHO];  /* maybe set param.v to point to v ?? */
   param.rhoe   = rhoe;
   NSCL_LOOP(nv) param.v[nv] = v[nv];
 
   #if NIONS == 0   /* No chemistry */
 
    rho = v[RHO];
    lnx = log10(rhoe);
    lny = log10(rho);
    i   = INT_FLOOR((lnx - Trhoe_tab.lnxmin)*Trhoe_tab.dlnx_1);
    j   = INT_FLOOR((lny - Trhoe_tab.lnymin)*Trhoe_tab.dlny_1);
 
 /* ------------------------------------------------------------------
     Since Trhoe_tab[j][i] increases with i but decrease with j,
     the min and max between four adjacent node values are the two
     opposite points of the secondary diagonal (i,j+1) and (i+1,j),
     respectively.
    ------------------------------------------------------------------ */
 
    Tmin = Trhoe_tab.f[j+1][i];
    Tmax = Trhoe_tab.f[j][i+1];
 
   #else           /* Chemistry */
    InternalEnergyBracket(rhoe, v, &Tmin, &Tmax);
   #endif
 
   if (Tmin < 0.0 || Tmax < 0.0) return 1;
 
 /*  status = Ridder(rhoeFunc, &param, Tmin, Tmax,-1, 1.0e-7, &T);  */
   status = Brent(rhoeFunc, &param, Tmin, Tmax, -1,   1.e-7,  T);
 
   if (*T < T_CUT_RHOE || status != 0) return 1;
   return 0;
 
 #endif /* TV_ENERGY_TABLE == NO */
 }

Here is the call graph for this function:

Here is the caller graph for this function:

double InternalEnergy	(	double *	v,
		double	T
	)

Compute and return the gas internal energy as a function of temperature and density or fractions: rhoe = rhoe(T,rho) or rhoe = rhoe(T,X) .

Parameters

[in]	v	1D Array of primitive variables containing density and species. Other variables are ignored.
[in]	T	Gas temperature

Returns: The gas internal energy (rhoe) in code units.

Definition at line 255 of file internal_energy.c.

 {
 #if TV_ENERGY_TABLE == YES
 
   int    status;
   double rho, rhoe;
 
   rho    = v[RHO];
   status = Table2DInterpolate(&rhoe_tab, T, rho, &rhoe);
   if (status != 0){
     print ("! InternalEnergy(): table interpolation failure (bound exceeded)\n");
     QUIT_PLUTO(1);
   }
   return rhoe;
 
 #else 
 
   return InternalEnergyFunc (v,T);
 
 #endif /* TV_ENERGY_TABLE == NO */
 }

Here is the call graph for this function:

Here is the caller graph for this function:

double InternalEnergyDeriv	(	double	rho,
		double	T
	)

static

Compute derivative of internal energy using numerical differentiation

Definition at line 147 of file internal_energy.c.

 {
   double Epp, Emm, Ep, Em, dEdT;
   double dT = 1.e-3, v[256];
 
   v[RHO] = rho;
   Epp   = InternalEnergyFunc(v,T*(1.0 + 2.0*dT));
   Ep    = InternalEnergyFunc(v,T*(1.0 + 1.0*dT));
   Em    = InternalEnergyFunc(v,T*(1.0 - 1.0*dT));
   Emm   = InternalEnergyFunc(v,T*(1.0 - 2.0*dT));
 
   dEdT = (-Epp + 8.0*Ep - 8.0*Em + Emm)/(12.0*T*dT);
   return dEdT;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void MakeEV_TemperatureTable ( )

Pre-calculate a 2D table of temperature T(rhoe,rho) by inverting, at specified values of p and rho, the nonlinear equation

$f(T) = \rho e - \rho e(T,\rho) = 0 \qquad\Longrightarrow\qquad T_{ij} = T(x_i, y_j) \qquad \left(x=\log_{10}\rho e\,,\quad y=\log_{10}\rho\right)$

For convenience, the table is constructed using equally spaced bins of x=Log(rhoe) and y=Log(rho) where rhoe and rho are in code units. This function must be called only once to initialize the table Trhoe_tab which is private to this file only.

This table is required only by the root-finder methods and it is useless for the tabulated version (default).

Definition at line 178 of file internal_energy.c.

 {
   int i,j, status;
   double rhoe, rho, Tlo, Thi, T;
   struct func_param par;
 
   print1 ("> MakeEV_TemperatureTable(): Generating table...\n");
 
   InitializeTable2D(&Trhoe_tab, 1.e-9, 1.e9, 1200, 1.e-12, 1.e12, 1200);
   Tlo = 1.0;
   Thi = 1.e12; 
 
   for (j = 0; j < Trhoe_tab.ny; j++){
   for (i = 0; i < Trhoe_tab.nx; i++){
     
     par.rhoe   = rhoe = Trhoe_tab.x[i];
     par.v[RHO] = rho  = Trhoe_tab.y[j];
 
 /*  status = Ridder(rhoeFunc, &param, Tlo, Thi, -1, 1.0e-12, &T);  */
     status = Brent(rhoeFunc, &par, Tlo, Thi, -1, 1.e-12, &T);
       
     if (status != 0) T = -1.0;
 
     Trhoe_tab.f[j][i] = T;
   }}
 
 
 /* let's plot the exact solution and the fundamental derivative  */
 #if 0 
 int ii;
 double G, lnT, v[256];
 FILE *fp;
 
 /* -- write exact internal energy, rhoe(T) -- */
 fp  = fopen("rhoe.dat","w");
 v[RHO] = rho = 1.245197084735e+00;
 for (lnT = 2.0; lnT < 4.5; lnT += 1.e-3){
   T    = pow(10.0, lnT);
   rhoe = InternalEnergyFunc(v,T);
   G    = FundamentalDerivative(v,T); 
   fprintf (fp,"%12.6e  %12.6e  %12.6e\n",T,rhoe,G);
 }
 fclose(fp);
 #endif
 
 /* -------------------------------------------------------
     Compute forward differences 
    ------------------------------------------------------- */
 
   FinalizeTable2D(&Trhoe_tab);
   if (prank == 0)  WriteBinaryTable2D("Trhoe_tab.bin",&Trhoe_tab);  
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void MakeInternalEnergyTable ( )

Compute a 2D table of the internal energy rhoe(T,rho) as function of temperature and density. The grid is equally spaced in log10(T) (=x) and log10(rho) (=y). Values are stored in the static table rhoe_tab.

Definition at line 62 of file internal_energy.c.

 {
   int i,j;
   double x, y, q;
   double T, rho, v[NVAR];
 
   print1 ("> MakeInternalEnergyTable(): Generating table (%d x %d points)\n",
            TV_ENERGY_TABLE_NX, TV_ENERGY_TABLE_NY);
 
   InitializeTable2D(&rhoe_tab, 1.0, 1.e8, TV_ENERGY_TABLE_NX, 
                                1.e-6, 1.e6, TV_ENERGY_TABLE_NY);  
   for (j = 0; j < rhoe_tab.ny; j++){
   for (i = 0; i < rhoe_tab.nx; i++){
     T   = rhoe_tab.x[i];
     rho = rhoe_tab.y[j];
     v[RHO] = rho;
     rhoe_tab.f[j][i] = InternalEnergyFunc(v,T);
   }}
 
   rhoe_tab.interpolation = SPLINE1;
 
 /* ----------------------------------------------------------------
     Compute cubic spline coefficients
    ---------------------------------------------------------------- */
 
   static double *dfdx;
 
   if (dfdx == NULL) dfdx = ARRAY_1D(TV_ENERGY_TABLE_NX, double);
 
   for (j = 0; j < rhoe_tab.ny; j++){
     rho = rhoe_tab.y[j]; 
     for (i = 0; i < rhoe_tab.nx; i++){
       T       = rhoe_tab.x[i]; 
       dfdx[i] = InternalEnergyDeriv(rho, T);
     }
     if (rhoe_tab.interpolation == SPLINE1){
       MonotoneSplineCoeffs(rhoe_tab.x, rhoe_tab.f[j], dfdx, rhoe_tab.nx,
                            rhoe_tab.a[j], rhoe_tab.b[j], rhoe_tab.c[j],
                            rhoe_tab.d[j]);
     }else if (rhoe_tab.interpolation == SPLINE2){
       SplineCoeffs(rhoe_tab.x, rhoe_tab.f[j],
                    dfdx[0], dfdx[rhoe_tab.nx-1], rhoe_tab.nx,
                    rhoe_tab.a[j], rhoe_tab.b[j], rhoe_tab.c[j], rhoe_tab.d[j]);
     }
   }
 
 /* let's plot the exact solutin and the spline using 20 points between
    adjacent nodes */
 #if 0 
 int ii;
 double G, lnT, rhoe, rhoe_ex;
 FILE *fp;
 
 j = TV_ENERGY_TABLE_NY/2;
 
 /* -- write exact internal energy, rhoe(T) -- */
 
 fp  = fopen("rhoe.dat","w");
 v[RHO] = rho = rhoe_tab.y[j];
 for (i = 10; i < rhoe_tab.nx-11; i++){
   for (ii = 0; ii < 20; ii++){   // Sub grid 
     T   = rhoe_tab.x[i] + ii*rhoe_tab.dx[i]/20.0;
     rhoe_ex = InternalEnergyFunc(v,T);
     G       = FundamentalDerivative(v,T); 
     Table2DInterpolate(&rhoe_tab, T, rho, &rhoe);
     fprintf (fp, "%12.6e  %12.6e  %12.6e  %12.6e\n",T,rhoe_ex, rhoe, G);
   }
 }
 fclose(fp);
 #endif
 
 
   FinalizeTable2D(&rhoe_tab);
   if (prank == 0) WriteBinaryTable2D("rhoe_tab.bin",&rhoe_tab);  
 }

Here is the call graph for this function:

Here is the caller graph for this function:

double rhoeFunc	(	double	T,
		void *	par
	)

static

Return the nonlinear function f(T)=rhoe(T,rho)-rhoe used by the root finder to find temperature T.

Parameters

[in]	T	The temperature in Kelvin
[in]	par	A pointer to a `func_param` structure containing quantities required for the computation.

Returns: f(T)

Definition at line 390 of file internal_energy.c.

 {
   double func, rhoe;
   struct func_param *p = (struct func_param *) par;
 
   if (T < 0.0){
     printf ("! rhoeFunc: T = %12.6e\n",T);
     exit(1);
   }
   rhoe = InternalEnergyFunc(p->v, T);
   func = rhoe - p->rhoe;
   return func;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void SoundSpeed2	(	double **	v,
		double *	cs2,
		double *	h,
		int	beg,
		int	end,
		int	pos,
		Grid *	grid
	)

Define the square of the sound speed.

Parameters

[in]	v	1D array of primitive quantities
[out]	cs2	1D array containing the square of the sound speed
[in]	h	1D array of enthalpy values
[in]	beg	initial index of computation
[in]	end	final index of computation
[in]	pos	an integer specifying the spatial position inside the cell (only for spatially-dependent EOS)
[in]	grid	pointer to an array of Grid structures

Returns: This function has no return value.

Definition at line 414 of file internal_energy.c.

 {
   int  i;
 
   for (i = beg; i <= end; i++){
      cs2[i] = Gamma1(v[i])*v[i][PRS]/v[i][RHO];
 /*     cs2[i] = 1.6667*v[i][PRS]/v[i][RHO];  */
   }
 }

Here is the call graph for this function:

Variable Documentation

Table2D rhoe_tab

static

A 2D table containing pre-computed values of internal energy stored at equally spaced node values of Log(T) and Log(rho) .

Definition at line 57 of file internal_energy.c.

Table2D Trhoe_tab

static

A 2D table containing pre-computed values of temperature by inverting e=e(T,rho) at equally spaced node values of Log(rhoe) and Log(rho).

This is used only in conjunction with a root finder to bracket the zero in a narrower interval.

Definition at line 171 of file internal_energy.c.

Macros

Functions

Variables

Detailed Description

Macro Definition Documentation

Function Documentation

Variable Documentation