Revision d334fbbb65c972b40e6431562405f85416011757 authored by Tom Quinn on 30 January 2016, 05:07:15 UTC, committed by Tom Quinn on 30 January 2016, 05:07:15 UTC
Change-Id: Iabe4025b68f4b9355e20d71ac6bce1e72af28b65
1 parent e572b13
cooling_planet.c
#ifndef NOCOOLING
/*
* Cooling code originally written by James Wadsley, McMaster
* University for GASOLINE.
*
* Updated for ChaNGa by Isaac Backus, University of Washington
*/
/*
* Cooling functions for a planetary disk where the cooling time is
* proportional to the orbital time.
*
* ALGORITHM:
*
* Integrating over a time step (see CoolIntegrateEnergyCode) gives:
* E(t+dt) = exp(-dt/tcool) (E(t) - PdV*tcool) + PdV*tcool
* where:
* tcool = beta/omega
* beta is a dimensionless cooling time (see CoolAddParams) and omega is the
* approximate Keplerian angular velocity, calculated as:
* omega = sqrt(G*M/R**3)
* where M is the total mass of all star particles and R is the distance to
* the center of mass of star particles.
*
* NOTES:
*
* The convention is that Cool...() functions are public and cl..()
* functions are private (which could be enforced if we move this to
* C++). The COOL class will contain data which is constant across
* all processors, and will be contained in a Charm++ NodeGroup.
* The clDerivsData structure will contain data that changes from
* particle to particle.
*
* By convention, all cooling is done in CGS. The header file defines macros
* which convert between code and CGS units [cooling_planet.h]. For example,
* to convert energy into CGS, do :
* E = CoolCodeEnergyToErgPerGm(cl, E);
* Conversion factors are stored in a COOL struct (see clInitConstants)
*
* For a list of parameters to set, see CoolAddParams in cooling_planet.c
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#include "cooling.h"
/**
* @brief CoolInit initializes a COOL data structure for storing
* cooling constants and parameters
* @return a blank COOL data structure with nTableRead = 0
*/
COOL *CoolInit( )
{
COOL *cl;
cl = (COOL *) malloc(sizeof(COOL));
assert(cl!=NULL);
/* Internal Tables read from Files */
cl->nTableRead = 0;
return cl;
}
/**
* Per thread initialization of cooling data.
* @param cl Initialized COOL structure.
*/
clDerivsData *CoolDerivsInit(COOL *cl)
{
clDerivsData *Data;
assert(cl != NULL);
Data = malloc(sizeof(clDerivsData));
assert(Data != NULL);
return Data;
}
///Frees memory and deletes cl
void CoolFinalize(COOL *cl )
{
free(cl);
}
/**
* Deallocate memory for per-thread data.
*/
void CoolDerivsFinalize(clDerivsData *clData)
{
free(clData);
}
/**
* @brief clInitConstants sets physical constants and parameters
* for cooling
* @param cl a COOL data structure (see CoolInit)
* @param CoolParam an initialized COOLPARAM structure (see CoolAddParams)
*/
void clInitConstants(COOL *cl, double dGmPerCcUnit,
double dComovingGmPerCcUnit, double dErgPerGmUnit,
double dSecUnit, double dKpcUnit, COOLPARAM CoolParam)
{
assert(cl!=NULL);
/* Units */
cl->dGmPerCcUnit = dGmPerCcUnit;
cl->dComovingGmPerCcUnit = dComovingGmPerCcUnit;
cl->dErgPerGmUnit = dErgPerGmUnit;
cl->dSecUnit = dSecUnit;
cl->dErgPerGmPerSecUnit = cl->dErgPerGmUnit / cl->dSecUnit;
cl->diErgPerGmUnit = 1./dErgPerGmUnit;
cl->dKpcUnit = dKpcUnit;
/* config parameters */
cl->Y_Total = CoolParam.Y_Total;
cl->Tmin = CoolParam.dCoolingTmin;
cl->Tmax = CoolParam.dCoolingTmax;
cl->beta = CoolParam.dBetaCooling;
}
/**
* @brief Calculates thermal energy from temperature and Y (deg. freedom/3)
* @param Y_Total (deg. freedom)/3
* @param T Temperature
* @return Thermal (internal) energy per mass (CGS)
*/
double clThermalEnergy( double Y_Total, double T ) {
return Y_Total*CL_Eerg_gm_degK3_2*T;
}
/**
* @brief Calculates temperature from internal energy and Y (deg.freedom/3)
* @param Y_Total (deg. freedom)/3
* @param E Internal energy per mass (CGS units)
* @return Temperature
*/
double clTemperature( double Y_Total, double E ) {
return E/(Y_Total*CL_Eerg_gm_degK3_2);
}
/* Module Interface routines */
/**
* @brief CoolAddParams Parses parameters and adds them to the COOLPARAM
* struct
* @param CoolParam A COOLPARAM struct
* @param prm The param struct
*
* The available parameters to be set are:
* dBetaCooling : Effective cooling time = t_cool * omega. Also
* equal to 2*pi*(t_cool/T_orbit)
* dY_Total : (degrees of freedom)/3 for the gas
* dCoolingTmin : Minimum allowed temperature (in Kelvin)
* dCoolingTmax : Maximum allowed temperature (in Kelvin)
*/
void CoolAddParams( COOLPARAM *CoolParam, PRM prm ) {
CoolParam->dBetaCooling = 1.0;
prmAddParam(prm, "dBetaCooling", paramDouble, &CoolParam->dBetaCooling,
sizeof(double), "betacool",
"<Effective cooling time (tCool*omega)> = 1");
CoolParam->Y_Total = 2;
prmAddParam(prm,"dY_Total",paramDouble,&CoolParam->Y_Total,
sizeof(double),"Y_Total",
"<Y_Total> = 2");
CoolParam->dCoolingTmin = 0;
prmAddParam(prm,"dCoolingTmin",paramDouble,&CoolParam->dCoolingTmin,
sizeof(double),"ctmin",
"<Minimum Temperature for Cooling> = 0K");
CoolParam->dCoolingTmax = 1e9;
prmAddParam(prm,"dCoolingTmax",paramDouble,&CoolParam->dCoolingTmax,
sizeof(double),"ctmax",
"<Maximum Temperature for Cooling> = 1e9K");
}
/* Placeholder functions which just need to be defined to make
* Sph.C and the compiler happy */
void CoolTableReadInfo( COOLPARAM *CoolParam, int cntTable, int *nTableColumns,
char *suffix )
{
*nTableColumns = 0;
}
void CoolTableRead( COOL *Cool, int nData, void *vData) {}
void CoolInitRatesTable( COOL *cl, COOLPARAM CoolParam ) {}
double CoolEdotInstantCode(COOL *cl, COOLPARTICLE *cp, double ECode,
double rhoCode, double ZMetal, double *posCode ) {}
/* Initialization Routines */
/**
* @brief Initializes data and sets defaults for a cooling particle
* @param cp a COOLPARTICLE structure (see cooling_planet.h)
*/
void CoolDefaultParticleData( COOLPARTICLE *cp )
{
cp->Y_Total = 2;
}
/**
* @brief Initializes data for a cooling particle, given a COOL struct
* @param cl a COOL struct which has been initialized with constants defined
* @param cp a cooling particle COOLPARTICLE struct
* @param E Pointer to the particle's internal energy
*
* Initializes data for a cooling particle, given a COOL struct.
* Sets Y_total and E in the cooling particle
*/
void CoolInitEnergyAndParticleData( COOL *cl, COOLPARTICLE *cp, double *E,
double dDensity, double dTemp, double fMetal)
{
cp->Y_Total = cl->Y_Total;
*E = clThermalEnergy(cp->Y_Total,dTemp)*cl->diErgPerGmUnit;
}
///A hold over from cosmology. Sets the time which is available to all threads
void CoolSetTime( COOL *cl, double dTime, double z ) {
cl->z = z;
cl->dTime = dTime;
}
/**
* @brief CoolSetStarCM saves the center of mass of the star(s) to the COOL
* struct cl
* @param cl COOL struct, which stores cooling data/info
* @param dCenterOfMass Array(length 4) which contains the star(s) center of
* mass as the first 3 entries and the total star mass as the final entry
*/
void CoolSetStarCM(COOL *cl, double dCenterOfMass[4]) {
int i;
for(i=0; i<4; ++i) {
cl->dStarCenterOfMass[i] = dCenterOfMass[i];
}
}
///Calculates the temperature of a COOLPARTICLE given its energy (in Code units)
double CoolCodeEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E,
double fMetal) {
return clTemperature(cp->Y_Total, E*Cool->dErgPerGmUnit);
}
/**
* Updates a particle's internal energy based on cooling and PdV work.
* Cooling is done according to the equation:
* Edot = -E/tcool + PdV
* where PdV is the PdV work per time done on the particle and tcool
* is given by:
* tcool = beta/omega
* beta is a dimensionless cooling set by the runtime parameter dBetaCooling
* and omega is the approximate keplerian orbital angular velocity
* Solving the cooling equation over a time interval from 0 to t gives:
* E(t) = exp(-t/tcool) (E0 - PdV*tcool) + PdV*tcool
*
* The Energy (Ecode) is updated
*
* @brief CoolIntegrateEnergyCode Integrates a particle's internal energy
* due to cooling and PdV work
* @param cl The COOL structure containing constants/params for the simulation
* @param clData [unused]
* @param cp Cooling particle
* @param ECode Particle internal energy in Code units
* @param PdVCode PdV work per time in code units
* @param rhoCode [unused]
* @param ZMetal [Unused]
* @param posCode particle position in code units
* @param tStep Time step to integrate over, in seconds
*/
void CoolIntegrateEnergyCode(COOL *cl, clDerivsData *clData, COOLPARTICLE *cp,
double *ECode, double PdVCode, double rhoCode,
double ZMetal, double *posCode, double tStep )
{
double radius = 0;
double omega, tcool, PdV, T, dECGS;
/* This is used as a switch to disable cooling before a certain time*/
if (tStep <= 0) return;
/* Calculate the radius (in simulation units)
* Defined as the distance to the center of mass of all the star
* particles */
int i=0;
for(i=0; i<3; ++i) {
radius += pow(posCode[i] - cl->dStarCenterOfMass[i], 2);
}
radius = sqrt(radius);
/* Calculate approximate keplerian angular velocity */
omega = sqrt(cl->dStarCenterOfMass[3] / pow(radius,3));
/* Calculate cooling time (in seconds)*/
tcool = cl->dSecUnit * cl->beta/omega;
/* Convert the particle's internal energy to CGS */
dECGS = CoolCodeEnergyToErgPerGm( cl, *ECode );
/* Convert PdV work to CGS */
PdV = CoolCodeWorkToErgPerGmPerSec(cl, PdVCode);
/* Integrate energy */
dECGS = exp(-tStep/tcool) * (dECGS - PdV*tcool) + PdV*tcool;
/* apply minimum temperature cutoff */
T = clTemperature(cp->Y_Total, dECGS);
if (T < cl->Tmin) {
T = cl->Tmin;
dECGS = clThermalEnergy(cp->Y_Total, T);
}
/*Convert energy back to simulation units */
*ECode = CoolErgPerGmToCodeEnergy(cl, dECGS);
}
#endif /* NOCOOLING */
Computing file changes ...