Skip to content

Commit

Permalink
Merge pull request #184 from N-BodyShop/public_planetesimal
Browse files Browse the repository at this point in the history 
Public planetesimal
  • Loading branch information
@trquinn
trquinn authored Oct 7, 2024
2 parents b0736f5 + 946cd58 commit a05be57
Show file tree
Hide file tree
Showing 27 changed files with 2,248 additions and 124 deletions.
51 changes: 50 additions & 1 deletion GravityParticle.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,10 @@ class extraSPHData
double _fMFracIron; /* Iron mass fraction */
double _fESNrate; /* SN energy rate */
double _fTimeCoolIsOffUntil;/* time cooling is turned back on */
#ifndef COLLISION
Vector3D<cosmoType> _vPred; /* Predicted velocities for velocity
dependent forces */
#endif
double _uPred; /* Predicted internal energy */
double _divv; /* Diverence of the velocity */
Vector3D<double> _curlv; /* Curl of the velocity */
Expand Down Expand Up @@ -123,7 +125,9 @@ class extraSPHData
inline double& fMFracIron() {return _fMFracIron;}
inline double& fESNrate() {return _fESNrate;}
inline double& fTimeCoolIsOffUntil() {return _fTimeCoolIsOffUntil;}
#ifndef COLLISION
inline Vector3D<cosmoType>& vPred() {return _vPred;}
#endif
inline double& uPred() {return _uPred;}
inline double& divv() {return _divv;}
inline Vector3D<double>& curlv() {return _curlv;}
Expand Down Expand Up @@ -185,7 +189,9 @@ class extraSPHData
p | _fMFracOxygen;
p | _fESNrate;
p | _fTimeCoolIsOffUntil;
#ifndef COLLISION
p | _vPred;
#endif
p | _uPred;
p | _divv;
p | _curlv;
Expand Down Expand Up @@ -314,9 +320,23 @@ class ExternalSmoothParticle;
/// Other classes of particles require this plus an "extra data" class.

class GravityParticle : public ExternalGravityParticle {
#ifdef COLLISION
private:
// Necessary for gas drag forces on solid bodies
Vector3D<cosmoType> _vPred;
#endif

public:
Vector3D<cosmoType> velocity;
Vector3D<cosmoType> treeAcceleration;
#ifdef COLLISION
cosmoType dtCol;
int64_t iOrderCol;
Vector3D<cosmoType> w;
cosmoType dtKep;
// For consistency with non-collision version
inline Vector3D<cosmoType>& vPred() {return _vPred;}
#endif
cosmoType potential;
cosmoType dtGrav; ///< timestep from gravity; N.B., this
/// is actually stored as (1/time^2)
Expand Down Expand Up @@ -365,6 +385,13 @@ class GravityParticle : public ExternalGravityParticle {
p | key;
p | velocity;
p | treeAcceleration;
#ifdef COLLISION
p | dtCol;
p | iOrderCol;
p | _vPred;
p | dtKep;
p | w;
#endif
p | dtGrav;
p | fDensity;
p | fBall;
Expand Down Expand Up @@ -412,7 +439,9 @@ class GravityParticle : public ExternalGravityParticle {
inline double& fMFracIron() {IMAGAS; return (((extraSPHData*)extraData)->fMFracIron());}
inline double& fESNrate() {IMAGAS; return (((extraSPHData*)extraData)->fESNrate());}
inline double& fTimeCoolIsOffUntil() {IMAGAS; return (((extraSPHData*)extraData)->fTimeCoolIsOffUntil());}
#ifndef COLLISION
inline Vector3D<cosmoType>& vPred() { IMAGAS; return (((extraSPHData*)extraData)->vPred());}
#endif
inline double& uPred() {IMAGAS; return (((extraSPHData*)extraData)->uPred());}
inline double& divv() { IMAGAS; return (((extraSPHData*)extraData)->divv());}
inline Vector3D<double>& curlv() { IMAGAS; return (((extraSPHData*)extraData)->curlv());}
Expand Down Expand Up @@ -579,6 +608,11 @@ class ExternalSmoothParticle {
#ifdef DTADJUST
double dt;
double dtNew;
#endif
#ifdef COLLISION
double soft;
double dtCol;
int64_t iOrderCol;
#endif
Vector3D<cosmoType> vPred;
Vector3D<cosmoType> treeAcceleration;
Expand Down Expand Up @@ -642,6 +676,11 @@ class ExternalSmoothParticle {
iType = p->iType;
rung = p->rung;
treeAcceleration = p->treeAcceleration;
#ifdef COLLISION
soft = p->soft;
dtCol = p->dtCol;
iOrderCol = p->iOrderCol;
#endif
if(TYPETest(p, TYPE_GAS)) {
vPred = p->vPred();
mumax = p->mumax();
Expand Down Expand Up @@ -711,6 +750,11 @@ class ExternalSmoothParticle {
tmp->iType = iType;
tmp->rung = rung;
tmp->treeAcceleration = treeAcceleration;
#ifdef COLLISION
tmp->soft = soft;
tmp->dtCol = dtCol;
tmp->iOrderCol = iOrderCol;
#endif
if(TYPETest(tmp, TYPE_GAS)) {
tmp->vPred() = vPred;
tmp->mumax() = mumax;
Expand Down Expand Up @@ -773,6 +817,7 @@ class ExternalSmoothParticle {
void pup(PUP::er &p) {
p | position;
p | velocity;
p | vPred;
p | mass;
p | fBall;
p | fDensity;
Expand All @@ -784,7 +829,6 @@ class ExternalSmoothParticle {
p | dtNew;
#endif
p | treeAcceleration;
p | vPred;
p | mumax;
p | PdV;
p | uDotPdV;
Expand Down Expand Up @@ -831,6 +875,11 @@ class ExternalSmoothParticle {
p | iEaterOrder;
p | dTimeFB;
p | iBucketOff;
#ifdef COLLISION
p | soft;
p | dtCol;
p | iOrderCol;
#endif
}
#endif
};
Expand Down
34 changes: 34 additions & 0 deletions InOutput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,13 @@ void load_tipsy_dark(Tipsy::TipsyReader &r, GravityParticle &p)
#endif
p.fDensity = 0.0;
p.iType = TYPE_DARK;
#ifdef COLLISION
p.dtKep = 0;
p.vPred() = dp.vel;
p.dtCol = DBL_MAX;
p.iOrderCol = -1;
p.w = Vector3D<double>(0.);
#endif
}

template <typename TPos, typename TVel>
Expand All @@ -106,6 +113,11 @@ void load_tipsy_star(Tipsy::TipsyReader &r, GravityParticle &p)
p.soft = sp.eps;
#ifdef CHANGESOFT
p.fSoft0 = sp.eps;
#endif
#ifdef COLLISION
p.dtCol = DBL_MAX;
p.iOrderCol = -1;
p.w = Vector3D<double>(0.);
#endif
p.fDensity = 0.0;
p.iType = TYPE_STAR;
Expand Down Expand Up @@ -713,6 +725,28 @@ static void load_NC_dark(std::string filename, int64_t startParticle,

load_NC_base(filename, startParticle, myNumDark, myParts);

#ifdef COLLISION
if(verbosity && startParticle == 0)
CkPrintf("loading spins\n");

FieldHeader fh;
void *data = readFieldData(filename + "/spin", fh, 3, myNumDark,
startParticle);
for(int i = 0; i < myNumDark; ++i) {
switch(fh.code) {
case float32:
myParts[i].w = static_cast<Vector3D<float> *>(data)[i];
break;
case float64:
myParts[i].w = static_cast<Vector3D<double> *>(data)[i];
break;
default:
throw XDRException("I don't recognize the type of this field!");
}
}
deleteField(fh, data);
#endif

for(int i = 0; i < myNumDark; ++i) {
myParts[i].fDensity = 0.0;
myParts[i].iType = TYPE_DARK;
Expand Down
29 changes: 29 additions & 0 deletions InOutput.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,35 @@ class PotOutputParams : public OutputParams
}
};


#ifdef COLLISION
/// @brief Output particle spin
class SpinOutputParams : public OutputParams
{
public:
virtual double dValue(GravityParticle *p) {CkAssert(0); return 0.0;}
virtual Vector3D<double> vValue(GravityParticle *p)
{return p->w;}
virtual void setDValue(GravityParticle *p, double val) {CkAssert(0);}
virtual void setVValue(GravityParticle *p, Vector3D<double> val)
{p->position = val;}
virtual int64_t iValue(GravityParticle *p) {CkAssert(0); return 0.0;}
virtual void setIValue(GravityParticle *p, int64_t iValue) {CkAssert(0);}
SpinOutputParams() {}
SpinOutputParams(std::string _fileName, int _iBinaryOut, double _dTime) {
bFloat = 1;
bVector = 1; fileName = _fileName; iBinaryOut = _iBinaryOut;
sTipsyExt = "spin"; sNChilExt = "spin";
dTime = _dTime;
iType = TYPE_DARK; }
PUPable_decl(SpinOutputParams);
SpinOutputParams(CkMigrateMessage *m) {}
virtual void pup(PUP::er &p) {
OutputParams::pup(p);//Call base class
}
};
#endif

/// @brief Output particle gas density
class GasDenOutputParams : public OutputParams
{
Expand Down
2 changes: 1 addition & 1 deletion Makefile.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ defines := $(strip @HEXADECAPOLE@ @FLAG_GDFORCE@ @FLAG_ARCH@ \
@FLAG_SPH_KERNEL@ @FLAG_COOLING@ @FLAG_DIFFUSION@ \
@FLAG_FEEDBACKDIFFLIMIT@ @FLAG_CULLENALPHA@ \
@FLAG_VSIGVISC@ @FLAG_DAMPING@ @FLAG_DIFFHARMONIC@ \
@FLAG_JEANSSOFTONLY@ @FLAG_FLOAT@ \
@FLAG_JEANSSOFTONLY@ @FLAG_FLOAT@ @FLAG_COLLISION@ \
@FLAG_TREE_BUILD@ $(debug_defines) @FLAG_INTERLIST@ \
@FLAG_NSMOOTHINNER@ @FLAG_SPLITGAS@ @FLAG_SIDMINTERACT@ \
@FLAG_SUPERBUBBLE@ $(cuda_defines) -DREDUCTION_HELPER)
Expand Down
50 changes: 42 additions & 8 deletions ParallelGravity.ci
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,10 +145,17 @@ mainmodule ParallelGravity {
#endif
//SIDM
PUPable SIDMSmoothParams;
#ifdef COLLISION
PUPable CollisionSmoothParams;
#endif

PUPable MassOutputParams;
PUPable PosOutputParams;
PUPable VelOutputParams;
PUPable PotOutputParams;
#ifdef COLLISION
PUPable SpinOutputParams;
#endif
PUPable GasDenOutputParams;
PUPable TempOutputParams;
PUPable AccOutputParams;
Expand Down Expand Up @@ -375,20 +382,32 @@ mainmodule ParallelGravity {
double dMultiPhaseMaxTime, double dMultiPhaseMinTemp, double dEvapCoeff, const CkCallback& cb);
entry void initAccel(int iKickRung, const CkCallback& cb);
entry void applyFrameAcc(int iKickRung, Vector3D<double> frameAcc, const CkCallback& cb);
entry void externalGravity(int activeRung, const ExternalGravity exGrav,
entry void externalForce(int activeRung, const ExternalForce& exForce, int bKepStep,
const CkCallback& cb);
#ifdef COOLING_MOLECULARH
entry void distribLymanWerner(const CkCallback& cb);
#endif
entry void adjust(int iKickRung, int bEpsAccStep, int bGravStep,
int bSphStep, int bViscosityLimitdt,
double dEta, double dEtaCourant, double dEtauDot,
double dDiffCoeff, double dEtaDiffusion,
double dDelta, double dAccFac,
double dCosmoFac, double dhMinOverSoft,
double dResolveJeans,
#ifdef COLLISION
entry void adjust(int iKickRung, int bCollStep, int bEpsAccStep,
int bGravStep, int bKepStep, int bSphStep,
int bViscosityLimitdt, double dEta, double dEtaCourant,
double dEtauDot, double dDiffCoeff, double dEtaDiffusion,
double dDelta, double dAccFac,
double dCosmoFac, double dhMinOverSoft,
double dResolveJeans,
int bDoGas,
const CkCallback& cb);
#else
entry void adjust(int iKickRung, int bEpsAccStep,
int bGravStep, int bSphStep,
int bViscosityLimitdt, double dEta, double dEtaCourant,
double dEtauDot, double dDiffCoeff, double dEtaDiffusion,
double dDelta, double dAccFac,
double dCosmoFac, double dhMinOverSoft,
double dResolveJeans,
int bDoGas,
const CkCallback& cb);
#endif
entry void truncateRung(int iCurrMaxRung, const CkCallback& cb);
entry void rungStats(const CkCallback& cb);
entry void countActive(int activeRung, const CkCallback& cb);
Expand Down Expand Up @@ -441,6 +460,21 @@ mainmodule ParallelGravity {
entry void Feedback(Fdbk fb, double dTime, double dDelta,
const CkCallback& cb);
entry void massMetalsEnergyCheck(int bPreDist, const CkCallback& cb);
#ifdef COLLISION
entry void delEjected(double dDelDist, const CkCallback& cb);
entry void getNearCollPartners(const CkCallback& cb);
entry void logOverlaps(const CkCallback& cb);
entry void getCollInfo(const CkCallback& cb);
entry void getCollInfo(int64_t iOrder, const CkCallback& cb);
entry void resolveCollision(Collision coll, const ColliderInfo &c1, const ColliderInfo &c2,
double baseStep, double timeNow, double dCentMass, const CkCallback& cb);
entry void sameHigherRung(int64_t iord1, int rung1, int64_t iord2, int rung2, const CkCallback& cb);
entry void resolveWallCollision(Collision coll, const ColliderInfo &c1, const CkCallback& cb);
entry void unKickCollStep(int iKickRung, double dDeltaBase, const CkCallback& cb);
entry void placeOnCollRung(int64_t iOrder, int collStepRung, const CkCallback& cb);
entry void resetRungs(const CkCallback& cb);
entry void getNeedCollStep(int collStepRung, const CkCallback& cb);
#endif
entry void setTypeFromFile(int iSetMask, char file[256],
const CkCallback& cb);
entry void getCOM(const CkCallback& cb, int bLiveViz);
Expand Down
Loading

0 comments on commit a05be57

Please sign in to comment.