https://github.com/N-BodyShop/changa
Tip revision: 98c212cc2ac931baa082740ce37b03f09eef89a4 authored by Tom Quinn on 05 August 2022, 19:44:28 UTC
Update version info for v3.5 release.
Update version info for v3.5 release.
Tip revision: 98c212c
Sph.h
/* Classes to describe smooth functions needed for SPH */
#ifndef __SPH_H
#define __SPH_H
#ifdef DIFFUSION
#if defined(FEEDBACKDIFFLIMIT) && !defined(DIFFUSIONHARMONIC)
#define DIFFUSIONHARMONIC
#endif
#endif
/// @brief Parameters and functions for the first SPH smooth: density
/// and velocity derivatives.
class DenDvDxSmoothParams : public SmoothParams
{
protected:
double a, H; // Cosmological parameters
int bActiveOnly;
int bConstantDiffusion;
double dTime;
double dAlphaMax; ///< Maximum SPH alpha
int bStarting; ///< We are starting (or restarting)
/// the simulation
int bHaveAlpha; ///< Alpha has been read in.
virtual void fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nList);
virtual int isSmoothActive(GravityParticle *p);
virtual void initTreeParticle(GravityParticle *p);
virtual void postTreeParticle(GravityParticle *p);
virtual void initSmoothParticle(GravityParticle *p);
virtual void initSmoothCache(GravityParticle *p);
virtual void combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2);
public:
DenDvDxSmoothParams() {}
/// @param _iType Type of particle to operate on
/// @param am Active rung
/// @param csm Cosmology information
/// @param _dTime Current time
/// @param _bActiveOnly Only operate on active particles.
/// @param _bConstantDiffusion Fixed diffusion constant
/// @param _bStarting Simulation is starting
/// @param _bHaveAlpha No need to calculate alpha
/// @param _dAlphaMax Maximum SPH alpha
DenDvDxSmoothParams(int _iType, int am, CSM csm, double _dTime,
int _bActiveOnly, int _bConstantDiffusion,
int _bStarting, int _bHaveAlpha, double _dAlphaMax) {
iType = _iType;
activeRung = am;
bActiveOnly = _bActiveOnly;
bConstantDiffusion = _bConstantDiffusion;
bStarting = _bStarting;
bHaveAlpha = _bHaveAlpha;
dAlphaMax = _dAlphaMax;
dTime = _dTime;
if(csm->bComove) {
H = csmTime2Hub(csm,dTime);
a = csmTime2Exp(csm,dTime);
}
else {
H = 0.0;
a = 1.0;
}
}
PUPable_decl(DenDvDxSmoothParams);
DenDvDxSmoothParams(CkMigrateMessage *m) : SmoothParams(m) {}
virtual void pup(PUP::er &p) {
SmoothParams::pup(p);//Call base class
p|dTime;
p|a;
p|H;
p|bActiveOnly;
p|bConstantDiffusion;
p|bStarting;
p|bHaveAlpha;
p|dAlphaMax;
}
};
/// @brief Get density and velocity derivatives of "Neighbor of
/// Active" particles
///
/// Like the above class, but only does non active particles marked
/// "Neighbor of Active" for the "fast gas" option.
/// Also, marking of neighbors is not complete. This is used in the
/// "FastGas" step. It uses the same fcnSmooth() as DenDvDxSmoothParams.
class DenDvDxNeighborSmParams : public DenDvDxSmoothParams
{
virtual int isSmoothActive(GravityParticle *p);
virtual void initTreeParticle(GravityParticle *p) {}
virtual void postTreeParticle(GravityParticle *p);
virtual void initSmoothParticle(GravityParticle *p) {}
virtual void initSmoothCache(GravityParticle *p) {}
virtual void combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2) {}
public:
DenDvDxNeighborSmParams() {}
/// @param _iType Type of particle to operate on
/// @param am Active rung
/// @param csm Cosmology information
/// @param dTime Current time
/// @param bConstantDiffusion Fixed diffusion constant
/// @param dAlphaMax Maximum SPH alpha
/// This calls the DenDvDx constructor, and we assume bActiveOnly,
/// bStarting, and bHaveAlpha are not set.
DenDvDxNeighborSmParams(int _iType, int am, CSM csm, double dTime,
int bConstDiffusion, double dAlphaMax)
: DenDvDxSmoothParams(_iType, am, csm, dTime, 0, bConstDiffusion,
0, 0, dAlphaMax) {}
PUPable_decl(DenDvDxNeighborSmParams);
DenDvDxNeighborSmParams(CkMigrateMessage *m) : DenDvDxSmoothParams(m) {}
virtual void pup(PUP::er &p) {
DenDvDxSmoothParams::pup(p);//Call base class
}
};
/// @brief Parameters for "Mark Smooth", used to find inverse nearest
/// neighbors.
class MarkSmoothParams : public SmoothParams
{
virtual void fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nList) {}
virtual int isSmoothActive(GravityParticle *p);
virtual void initTreeParticle(GravityParticle *p) {}
virtual void postTreeParticle(GravityParticle *p) {}
virtual void initSmoothParticle(GravityParticle *p) {}
virtual void initSmoothCache(GravityParticle *p) {}
virtual void combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2);
public:
MarkSmoothParams() {}
/// @param _iType Type of particle to operate on
/// @param am Active rung
MarkSmoothParams(int _iType, int am) {
iType = _iType;
activeRung = am;
}
PUPable_decl(MarkSmoothParams);
MarkSmoothParams(CkMigrateMessage *m) : SmoothParams(m) {}
virtual void pup(PUP::er &p) {
SmoothParams::pup(p);//Call base class
}
};
/// @brief Second pass in SPH: calculate pressure forces.
class PressureSmoothParams : public SmoothParams
{
double dTime;
double a, H; // Cosmological parameters
double alpha, beta; // SPH viscosity parameters
double dThermalDiffusionCoeff;
double dMetalDiffusionCoeff;
double dtFacCourant; // Courant timestep factor
double dtFacDiffusion; // Diffusion timestep factor
virtual void fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nList);
virtual int isSmoothActive(GravityParticle *p);
virtual void initTreeParticle(GravityParticle *p) {}
virtual void postTreeParticle(GravityParticle *p) {}
virtual void initSmoothParticle(GravityParticle *p);
virtual void initSmoothCache(GravityParticle *p);
virtual void combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2);
public:
PressureSmoothParams() {}
/// @param _iType Type of particles to smooth
/// @param am Active rung
/// @param csm Cosmological parameters
/// @param dTime Current time
/// @param _alpha Artificial viscosity parameter
/// @param _beta Artificial viscosity parameter
PressureSmoothParams(int _iType, int am, CSM csm, double _dTime,
double _alpha, double _beta,
double _dThermalDiff, double _dMetalDiff,
double dEtaCourant, double dEtaDiffusion) {
iType = _iType;
activeRung = am;
dTime = _dTime;
if(csm->bComove) {
H = csmTime2Hub(csm,dTime);
a = csmTime2Exp(csm,dTime);
}
else {
H = 0.0;
a = 1.0;
}
alpha = _alpha;
beta = _beta;
dThermalDiffusionCoeff = _dThermalDiff;
dMetalDiffusionCoeff = _dMetalDiff;
dtFacCourant = dEtaCourant*a*2.0/1.6;
dtFacDiffusion = 2.0*dEtaDiffusion;
}
PUPable_decl(PressureSmoothParams);
PressureSmoothParams(CkMigrateMessage *m) : SmoothParams(m) {}
virtual void pup(PUP::er &p) {
SmoothParams::pup(p);//Call base class
p|dTime;
p|a;
p|H;
p|alpha;
p|beta;
p|dThermalDiffusionCoeff;
p|dMetalDiffusionCoeff;
p|dtFacCourant;
p|dtFacDiffusion;
}
};
///
/// @brief SmoothParams class for distributing deleted gas to neighboring
/// particles.
///
class DistDeletedGasSmoothParams : public SmoothParams
{
virtual void fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nList);
virtual int isSmoothActive(GravityParticle *p);
virtual void initSmoothParticle(GravityParticle *p) {};
virtual void initTreeParticle(GravityParticle *p) {}
virtual void postTreeParticle(GravityParticle *p) {}
virtual void initSmoothCache(GravityParticle *p);
virtual void combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2);
public:
DistDeletedGasSmoothParams() {}
/// @param _iType Type of particles to smooth
/// @param am Active rung
DistDeletedGasSmoothParams(int _iType, int am) {
iType = _iType;
activeRung = am;
bUseBallMax = 0;
}
PUPable_decl(DistDeletedGasSmoothParams);
DistDeletedGasSmoothParams(CkMigrateMessage *m) : SmoothParams(m) {}
virtual void pup(PUP::er &p) {
SmoothParams::pup(p);//Call base class
}
};
#ifdef SUPERBUBBLE
///
/// @brief SmoothParams class for the sub-grid evaporation,
/// promoting cold neighbours to hot.
///
/// This method selects which particles cold particles are
/// to be evaporated.
///
class PromoteToHotGasSmoothParams : public SmoothParams
{
double dEvapCoeff;
double dEvapMinTemp;
double dErgPerGmUnit;
double dGmPerCcUnit;
double dDeltaStarForm;
double dTime;
virtual void fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nList);
virtual int isSmoothActive(GravityParticle *p);
virtual void initSmoothParticle(GravityParticle *p) ;
virtual void initTreeParticle(GravityParticle *p) ;
virtual void postTreeParticle(GravityParticle *p) {}
virtual void initSmoothCache(GravityParticle *p);
virtual void combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2) ;
public:
PromoteToHotGasSmoothParams() {}
/// @param _iType Type of particles to smooth
/// @param am Active rung
/// @param _dEvapCoeff evaporation/conduction rate coefficient
/// @param _dEvapMinTemp lowest temperature to allow evaporation
/// @param _dErgPerGmUnit specific energy in code units
/// @param _dGmPerCcUnit density in code units
/// @param _dDeltaStarForm starformation timestep
/// @param _dTime Current time
PromoteToHotGasSmoothParams(int _iType, int am, double _dEvapCoeff, double _dEvapMinTemp,
double _dErgPerGmUnit, double _dGmPerCcUnit, double _dDeltaStarForm, double _dTime) {
iType = _iType;
activeRung = am;
bUseBallMax = 0;
dEvapCoeff = _dEvapCoeff;
dEvapMinTemp = _dEvapMinTemp;
dErgPerGmUnit = _dErgPerGmUnit;
dGmPerCcUnit = _dGmPerCcUnit;
dDeltaStarForm = _dDeltaStarForm;
dTime = _dTime;
}
PUPable_decl(PromoteToHotGasSmoothParams);
PromoteToHotGasSmoothParams(CkMigrateMessage *m) : SmoothParams(m) {}
virtual void pup(PUP::er &p) {
SmoothParams::pup(p);//Call base class
p|dEvapCoeff;
p|dEvapMinTemp;
p|dErgPerGmUnit;
p|dGmPerCcUnit;
p|dDeltaStarForm;
p|dTime;
}
};
///
/// @brief SmoothParams class for the sub-grid evaporation,
/// promoting cold neighbours to hot.
///
/// This method takes cold particles marked to be promoted,
/// and averages their internal energy with the hot promoter.
///
class ShareWithHotGasSmoothParams : public SmoothParams
{
double dEvapMinTemp;
double dErgPerGmUnit;
double dGmPerCcUnit;
virtual void fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nList);
virtual int isSmoothActive(GravityParticle *p);
virtual void initSmoothParticle(GravityParticle *p) {}
virtual void initTreeParticle(GravityParticle *p) {}
virtual void postTreeParticle(GravityParticle *p) {}
virtual void initSmoothCache(GravityParticle *p);
virtual void combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2) ;
public:
ShareWithHotGasSmoothParams() {}
/// @param _iType Type of particles to smooth
/// @param am Active rung
/// @param _dEvapMinTemp lowest temperature to allow evaporation
/// @param _dErgPerGmUnit specific energy in code units
/// @param _dGmPerCcUnit density in code units
ShareWithHotGasSmoothParams(int _iType, int am, double _dEvapMinTemp,
double _dErgPerGmUnit, double _dGmPerCcUnit ){
iType = _iType;
activeRung = am;
bUseBallMax = 0;
dEvapMinTemp = _dEvapMinTemp;
dErgPerGmUnit = _dErgPerGmUnit;
dGmPerCcUnit = _dGmPerCcUnit;
}
PUPable_decl(ShareWithHotGasSmoothParams);
ShareWithHotGasSmoothParams(CkMigrateMessage *m) : SmoothParams(m) {}
virtual void pup(PUP::er &p) {
SmoothParams::pup(p);//Call base class
p|dEvapMinTemp;
p|dErgPerGmUnit;
p|dGmPerCcUnit;
}
};
#endif
#endif