Take a source step using power-law cooling. More...

#include "pluto.h"

Include dependency graph for cooling.c:

Functions
void	PowerLawCooling (Data_Arr VV, double dt, Time_Step Dts, Grid grid)

Detailed Description

Take a source step using power-law cooling.

Integrate the ODE

$dp_{\rm cgs}/dt_{\rm cgs} = -(\Gamma - 1) \Lambda(\rho_{\rm cgs}, T_{\rm cgs}) \qquad{\rm where}\qquad \Lambda = \frac{a_{br}}{(\mu m_H)^2} \rho_{\rm cgs}^2 \sqrt{T_{\rm cgs}} \,,\qquad a_{br} = 2.e-27 c.g.s$

which accounts for bremmstrahlung cooling. Here the subscript 'cgs' means that the corresponding quantity is given in c.g.s units. We denote with mu the molecular weight, mH the hydrogen mass (in c.g.s).

The previous integration is carried out analytically since the density does not change during this step. In non-dimensional form:

$\frac{dp}{dt} = -{\rm cost} \;\rho (\rho p)^{\HALF}$

[notice that since p/rho=T/KELVIN this is equivalent to: dp/dt=-cost rho^2 (T/KELVIN)^(1/2) ]

The quantity cost is determined by transforming the dimensional equation into the non-dimensional one. If p, rho and t are in code (non-dimensional) units and if L_0, rho_0, and V_0 are the unit length, density and velocity, then cost is found to be:

                 a_br * (gamma - 1) * L_0 * rho_0
      cost = -------------------------------------------
              sqrt(kB * mu * mH) * kB * mu * mH * V_0^2

where a_{br} = 2.e-27 (in c.g.s), kB is the Boltmann constant

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

Date: July 28, 2014

Definition in file cooling.c.

Function Documentation

void PowerLawCooling	(	Data_Arr	VV,
		double	dt,
		Time_Step *	Dts,
		Grid *	grid
	)

Parameters

[in,out]	VV	a pointer to the PLUTO 3D data array containing pimitive variables.
[in]	dt	the current integration time step
[in]	Dts	a pointer to the Time_Step structure
[in]	grid	pointer to an array of Grid structures

Definition at line 50 of file cooling.c.

 {
   int     i, j, k;
   double  cost, dE;
   double  rho, p, T, p_f, T_f;
 
   cost  = UNIT_LENGTH*UNIT_DENSITY/(UNIT_VELOCITY*UNIT_VELOCITY);
   cost *= 2.e-27*(g_gamma-1.0)/ (0.5*CONST_mH*sqrt(0.5*CONST_mH*CONST_kB));
 
 /*  -------------------------------------------------------------
                 Integrate analytically
     -------------------------------------------------------------  */
 
   dE = 1.e-18;
   DOM_LOOP(k,j,i){
 
 /*  ----  Find initial temperature in Kelvin  ----  */
 
     rho = VV[RHO][k][j][i];
     p   = VV[PRS][k][j][i];
 
     T   = (p/rho*KELVIN);
 
     if (T < g_minCoolingTemp) continue;
 
 /*  ----  Find final energy  ----  */
 
     p_f = sqrt(p) - 0.5*cost*rho*sqrt(rho)*dt;
     p_f = MAX(p_f, 0.0);
     p_f = p_f*p_f;
     T_f = p_f/rho*KELVIN;
     T_f = MAX (T_f, g_minCoolingTemp);
 
 /*  ----  Update Energy  ----  */
 
     p_f = T_f*rho/KELVIN;
 
     VV[PRS][k][j][i] = p_f;
     dE = fabs(1.0 - p_f/p) + 1.e-18;
 
     Dts->dt_cool = MIN(Dts->dt_cool, dt*g_maxCoolingRate/dE);
   }
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

Functions

Detailed Description

Function Documentation