diff --git a/GravityParticle.h b/GravityParticle.h index b9faa04b..698ae57f 100644 --- a/GravityParticle.h +++ b/GravityParticle.h @@ -61,6 +61,9 @@ class extraSPHData double _fTimeCoolIsOffUntil;/* time cooling is turned back on */ Vector3D _vPred; /* Predicted velocities for velocity dependent forces */ +#ifdef SLIDING_PATCH + double _dPyPred; ///< Predicted canonical momentum +#endif double _uPred; /* Predicted internal energy */ double _divv; /* Diverence of the velocity */ Vector3D _curlv; /* Curl of the velocity */ @@ -122,6 +125,9 @@ class extraSPHData inline double& fESNrate() {return _fESNrate;} inline double& fTimeCoolIsOffUntil() {return _fTimeCoolIsOffUntil;} inline Vector3D& vPred() {return _vPred;} +#ifdef SLIDING_PATCH + inline double& dPyPred() {return _dPyPred;} +#endif inline double& uPred() {return _uPred;} inline double& divv() {return _divv;} inline Vector3D& curlv() {return _curlv;} @@ -182,6 +188,9 @@ class extraSPHData p | _fESNrate; p | _fTimeCoolIsOffUntil; p | _vPred; +#ifdef SLIDING_PATCH + p | _dPyPred; +#endif p | _uPred; p | _divv; p | _curlv; @@ -311,6 +320,9 @@ class GravityParticle : public ExternalGravityParticle { public: SFC::Key key; Vector3D velocity; +#ifdef SLIDING_PATCH + double dPy; ///< Canonical momentum used to update y-velocity +#endif Vector3D treeAcceleration; cosmoType potential; cosmoType dtGrav; ///< timestep from gravity; N.B., this @@ -358,6 +370,9 @@ class GravityParticle : public ExternalGravityParticle { ExternalGravityParticle::pup(p); p | key; p | velocity; +#ifdef SLIDING_PATCH + p | dPy; +#endif p | treeAcceleration; p | dtGrav; p | fDensity; @@ -407,6 +422,9 @@ class GravityParticle : public ExternalGravityParticle { inline double& fESNrate() {IMAGAS; return (((extraSPHData*)extraData)->fESNrate());} inline double& fTimeCoolIsOffUntil() {IMAGAS; return (((extraSPHData*)extraData)->fTimeCoolIsOffUntil());} inline Vector3D& vPred() { IMAGAS; return (((extraSPHData*)extraData)->vPred());} +#ifdef SLIDING_PATCH + inline double& dPyPred() { IMAGAS; return (((extraSPHData*)extraData)->dPyPred());} +#endif inline double& uPred() {IMAGAS; return (((extraSPHData*)extraData)->uPred());} inline double& divv() { IMAGAS; return (((extraSPHData*)extraData)->divv());} inline Vector3D& curlv() { IMAGAS; return (((extraSPHData*)extraData)->curlv());} diff --git a/Makefile.in b/Makefile.in index 39caa3bf..3d1bc2ed 100644 --- a/Makefile.in +++ b/Makefile.in @@ -125,6 +125,7 @@ defines := $(strip @HEXADECAPOLE@ @FLAG_RTFORCE@ @FLAG_ARCH@ \ @FLAG_VSIGVISC@ @FLAG_DAMPING@ @FLAG_DIFFHARMONIC@ \ @FLAG_TREE_BUILD@ $(debug_defines) @FLAG_INTERLIST@ \ @FLAG_NSMOOTHINNER@ @FLAG_SPLITGAS@ @FLAG_SIDMINTERACT@ \ + @FLAG_SLIDING_PATCH@ \ @FLAG_SUPERBUBBLE@ $(cuda_defines) -DREDUCTION_HELPER) modules := $(strip -language charm++ -balancer @DEFAULT_LB@ \ diff --git a/ParallelGravity.ci b/ParallelGravity.ci index afeb7500..6fa81501 100644 --- a/ParallelGravity.ci +++ b/ParallelGravity.ci @@ -271,7 +271,7 @@ mainmodule ParallelGravity { entry void resetObjectLoad(const CkCallback& cb); entry void setPeriodic(int nReplicas, Vector3D fPeriod, int bEwald, double fEwCut, double fEwhCut, int bPeriod, - int bComove, double dRhoFac); + int bComove, double dRhoFac, double dOrbFreq); entry [notrace] void EwaldInit(); entry [notrace] void initCoolingData(const CkCallback& cb); entry void calculateEwald(dummyMsg *m); @@ -396,7 +396,8 @@ mainmodule ParallelGravity { entry void assignDomain(const CkCallback &cb); entry void drift(double dDelta, int bNeedVPred, int bGasIsoThermal, double dvDelta, double duDelta, int nGrowMass, - bool buildTree, double dMaxEnergy, const CkCallback& cb); + bool buildTree, double dMaxEnergy, + double dTime, const CkCallback& cb); entry void starCenterOfMass(const CkCallback& cb); entry void calcEnergy(const CkCallback& cb); entry void colNParts(const CkCallback &cb); diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index c46504e6..b33411ff 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -1971,7 +1971,7 @@ void Main::advanceBigStep(int iStep) { bool bBuildTree = (iSub + 1 == driftSteps); treeProxy.drift(dDriftFac, param.bDoGas, param.bGasIsothermal, dKickFac, dTimeSub, nGrowMassDrift, bBuildTree, - param.dMaxEnergy, + param.dMaxEnergy, dTime, CkCallbackResumeThread()); double tDrift = CkWallTimer() - startTime; timings[activeRung].tDrift += tDrift; @@ -2156,10 +2156,12 @@ void Main::advanceBigStep(int iStep) { void Main::setupICs() { double startTime; + param.externalGravity.CheckParams(prm, param); + treeProxy.setPeriodic(param.nReplicas, param.vPeriod, param.bEwald, param.dEwCut, param.dEwhCut, param.bPeriodic, param.csm->bComove, - 0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0); + 0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0, param.externalGravity.dOrbFreq); /******** Particles Loading ********/ CkPrintf("Loading particles ..."); @@ -2289,8 +2291,6 @@ void Main::setupICs() { } } - param.externalGravity.CheckParams(prm, param); - string achLogFileName = string(param.achOutName) + ".log"; ofstream ofsLog; if(bIsRestarting) @@ -2423,6 +2423,12 @@ void Main::setupICs() { #endif #ifdef DAMPING ofsLog << " DAMPING"; +#endif +#ifdef SLIDING_PATCH + ofsLog << " SLIDING_PATCH"; +#endif +#ifdef NO_HILL + ofsLog << " NO_HILL"; #endif ofsLog << endl; ofsLog << "# Key sizes: " << sizeof(KeyType) << " bytes particle " @@ -2487,7 +2493,7 @@ void Main::setupICs() { // for periodic, puts all particles within the boundary // Also assigns keys and sorts. - treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime, CkCallbackResumeThread()); initialForces(); @@ -2625,7 +2631,7 @@ Main::restart(CkCheckpointStatusMsg *msg) } else { CkPrintf("Not Using CkLoop %d\n", param.bUseCkLoopPar); } - treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime, CkCallbackResumeThread()); if(param.bGasCooling || param.bStarForm) initCooling(); @@ -2840,7 +2846,7 @@ Main::doSimulation() } // The following drift is called because it deletes the tree // so it won't be saved on disk. - treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, false, param.dMaxEnergy, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, false, param.dMaxEnergy, dTime, CkCallbackResumeThread()); treeProxy[0].flushStarLog(CkCallbackResumeThread()); param.iStartStep = iStep; // update so that restart continues on @@ -3033,7 +3039,7 @@ Main::doSimulation() if(param.bDoGas && param.bDoDensity) { // The following call is to get the particles in key order // before the sort. - treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3607,7 +3613,7 @@ void Main::writeOutput(int iStep) if(param.nSteps != 0 && param.bDoDensity) { // The following call is to get the particles in key order // before the sort. - treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3652,7 +3658,7 @@ void Main::writeOutput(int iStep) startTime = CkWallTimer(); // The following call is to get the particles in key order // before the sort. - treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3675,7 +3681,7 @@ void Main::writeOutput(int iStep) // processors for continuing the simulation. // The following call is to get the particles in key order // before the sort. - treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime, CkCallbackResumeThread()); domainDecomp(0); } diff --git a/ParallelGravity.h b/ParallelGravity.h index 3966a469..690837af 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -1192,6 +1192,7 @@ class TreePiece : public CBase_TreePiece { /// Background density of the Universe double dRhoFac; Vector3D fPeriod; + double dOrbFreq; ///< Orbital frequency of patch int nReplicas; int bEwald; /* Perform Ewald */ double fEwCut; @@ -1570,7 +1571,7 @@ class TreePiece : public CBase_TreePiece { void setPeriodic(int nReplicas, Vector3D fPeriod, int bEwald, double fEwCut, double fEwhCut, int bPeriod, - int bComove, double dRhoFac); + int bComove, double dRhoFac, double dOrbFreq); void BucketEwald(GenericTreeNode *req, int nReps,double fEwCut); void EwaldInit(); void calculateEwald(dummyMsg *m); @@ -1695,7 +1696,7 @@ class TreePiece : public CBase_TreePiece { double dMultiPhaseMaxTime, double dMultiPhaseMinTemp, double dEvapCoeff, const CkCallback& cb); void drift(double dDelta, int bNeedVPred, int bGasIsothermal, double dvDelta, double duDelta, int nGrowMass, bool buildTree, double dMaxEnergy, - const CkCallback& cb); + double dTime, const CkCallback& cb); void initAccel(int iKickRung, const CkCallback& cb); #ifdef COOLING_MOLECULARH void distribLymanWerner(const CkCallback& cb); @@ -1995,14 +1996,26 @@ class TreePiece : public CBase_TreePiece { // need this in TreeWalk GenericTreeNode *getRoot() {return root;} // need this in Compute + inline double SHEAR(int ix, ///< Interior or exterior? + double t, ///< Current simulation time + Vector3D fPeriod, ///< vector for dxPeriod and dyPeriod + double dOrbFreq) ///< Orbital frequency + { +#ifdef SLIDING_PATCH + return (ix < 0) ? fmod(0.5 * fPeriod[1] - 1.5 * ix * dOrbFreq * fPeriod[0] * t, fPeriod[1]) - 0.5 * fPeriod[1] : + (ix > 0) ? 0.5 * fPeriod[1] - fmod(0.5 * fPeriod[1] + 1.5 * ix * dOrbFreq * fPeriod[0] * t, fPeriod[1]) : 0.0; +#else + return 0; +#endif + } + inline Vector3D decodeOffset(int reqID) { int offsetcode = reqID >> 22; int x = (offsetcode & 0x7) - 3; int y = ((offsetcode >> 3) & 0x7) - 3; int z = ((offsetcode >> 6) & 0x7) - 3; - Vector3D offset(x*fPeriod.x, y*fPeriod.y, z*fPeriod.z); - + Vector3D offset(x*fPeriod.x, y*fPeriod.y + SHEAR(x, totalTime, fPeriod.y, dOrbFreq), z*fPeriod.z); return offset; } diff --git a/Sph.cpp b/Sph.cpp index 8c511252..1d3dd586 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -30,7 +30,7 @@ Main::initSph() // Starting is true DenDvDxSmoothParams pDen(TYPE_GAS, 0, param.csm, dTime, 0, param.bConstantDiffusion, 1, bHaveAlpha, - param.dConstAlphaMax); + param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod); double startTime = CkWallTimer(); double dfBall2OverSoft2 = 4.0*param.dhMinOverSoft*param.dhMinOverSoft; treeProxy.startSmooth(&pDen, 1, param.nSmooth, dfBall2OverSoft2, @@ -712,7 +712,7 @@ Main::doSph(int activeRung, int bNeedDensity) // This also marks neighbors of actives DenDvDxSmoothParams pDen(TYPE_GAS, activeRung, param.csm, dTime, 1, param.bConstantDiffusion, 0, 0, - param.dConstAlphaMax); + param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod); double startTime = CkWallTimer(); treeProxy.startSmooth(&pDen, 1, param.nSmooth, dfBall2OverSoft2, CkCallbackResumeThread()); @@ -732,7 +732,7 @@ Main::doSph(int activeRung, int bNeedDensity) // additional marking DenDvDxNeighborSmParams pDenN(TYPE_GAS, activeRung, param.csm, dTime, param.bConstantDiffusion, - param.dConstAlphaMax); + param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod); startTime = CkWallTimer(); treeProxy.startSmooth(&pDenN, 1, param.nSmooth, dfBall2OverSoft2, CkCallbackResumeThread()); @@ -745,7 +745,7 @@ Main::doSph(int activeRung, int bNeedDensity) // actives, and those who have actives as neighbors. DenDvDxSmoothParams pDen(TYPE_GAS, activeRung, param.csm, dTime, 0, param.bConstantDiffusion, 0, 0, - param.dConstAlphaMax); + param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod); double startTime = CkWallTimer(); treeProxy.startSmooth(&pDen, 1, param.nSmooth, dfBall2OverSoft2, CkCallbackResumeThread()); @@ -779,7 +779,7 @@ Main::doSph(int activeRung, int bNeedDensity) PressureSmoothParams pPressure(TYPE_GAS, activeRung, param.csm, dTime, param.dConstAlpha, param.dConstBeta, param.dThermalDiffusionCoeff, param.dMetalDiffusionCoeff, - param.dEtaCourant, param.dEtaDiffusion); + param.dEtaCourant, param.dEtaDiffusion, param.externalGravity.dOrbFreq, param.fPeriod); double startTime = CkWallTimer(); treeProxy.startReSmooth(&pPressure, CkCallbackResumeThread()); ckout << " took " << (CkWallTimer() - startTime) << " seconds." @@ -1119,6 +1119,14 @@ void DenDvDxSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, dvx = (-p->vPred().x + q->vPred().x)*vFac; dvy = (-p->vPred().y + q->vPred().y)*vFac; dvz = (-p->vPred().z + q->vPred().z)*vFac; +#ifdef SLIDING_PATCH + if (dx < 0.0 && (p->position.x - q->position.x > 0.0)) { + dvy -= 1.5 * dOrbFreq * fPeriod.x; + } + else if (dx > 0.0 && (p->position.x - q->position.x < 0.0)) { + dvy += 1.5 * dOrbFreq * fPeriod.x; + } +#endif dvxdx += dvx*dx*rs1; dvxdy += dvx*dy*rs1; dvxdz += dvx*dz*rs1; @@ -1575,6 +1583,14 @@ void PressureSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, qParams.rNorm = rs1 * q->mass; params.dx = nnList[i].dx; dv = p->vPred() - q->vPred(); +#ifdef SLIDING_PATCH + if (params.dx.x < 0.0 && (p->position.x - q->position.x > 0.0)) { + dv[1] += 1.5 * dOrbFreq * fPeriod.x; + } + else if (params.dx.x > 0.0 && (p->position.x - q->position.x < 0.0)) { + dv[1] -= 1.5 * dOrbFreq * fPeriod.x; + } +#endif params.dvdotdr = vFac*dot(dv, params.dx) + fDist2*H; #ifdef RTFORCE pParams.PoverRho2 = p->PoverRho2()*p->fDensity/q->fDensity; diff --git a/Sph.h b/Sph.h index 872010a6..41e11eaf 100644 --- a/Sph.h +++ b/Sph.h @@ -23,6 +23,8 @@ class DenDvDxSmoothParams : public SmoothParams int bStarting; ///< We are starting (or restarting) /// the simulation int bHaveAlpha; ///< Alpha has been read in. + double dOrbFreq; ///< Orbital Frequency + Vector3D fPeriod;///< Dimensions of patch virtual void fcnSmooth(GravityParticle *p, int nSmooth, pqSmoothNode *nList); @@ -46,7 +48,8 @@ class DenDvDxSmoothParams : public SmoothParams /// @param _dAlphaMax Maximum SPH alpha DenDvDxSmoothParams(int _iType, int am, CSM csm, double _dTime, int _bActiveOnly, int _bConstantDiffusion, - int _bStarting, int _bHaveAlpha, double _dAlphaMax) { + int _bStarting, int _bHaveAlpha, double _dAlphaMax, + double _dOrbFreq, Vector3D _fPeriod) { iType = _iType; activeRung = am; bActiveOnly = _bActiveOnly; @@ -55,6 +58,8 @@ class DenDvDxSmoothParams : public SmoothParams bHaveAlpha = _bHaveAlpha; dAlphaMax = _dAlphaMax; dTime = _dTime; + dOrbFreq = _dOrbFreq; + fPeriod = _fPeriod; if(csm->bComove) { H = csmTime2Hub(csm,dTime); a = csmTime2Exp(csm,dTime); @@ -76,6 +81,8 @@ class DenDvDxSmoothParams : public SmoothParams p|bStarting; p|bHaveAlpha; p|dAlphaMax; + p|dOrbFreq; + p|fPeriod; } }; @@ -107,9 +114,9 @@ class DenDvDxNeighborSmParams : public DenDvDxSmoothParams /// 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) + int bConstDiffusion, double dAlphaMax, double dOrbFreq, Vector3D _fPeriod) : DenDvDxSmoothParams(_iType, am, csm, dTime, 0, bConstDiffusion, - 0, 0, dAlphaMax) {} + 0, 0, dAlphaMax, dOrbFreq, fPeriod) {} PUPable_decl(DenDvDxNeighborSmParams); DenDvDxNeighborSmParams(CkMigrateMessage *m) : DenDvDxSmoothParams(m) {} virtual void pup(PUP::er &p) { @@ -156,6 +163,8 @@ class PressureSmoothParams : public SmoothParams double dMetalDiffusionCoeff; double dtFacCourant; // Courant timestep factor double dtFacDiffusion; // Diffusion timestep factor + double dOrbFreq; // Orbital Frequency + Vector3D fPeriod; // Dimensions of patch virtual void fcnSmooth(GravityParticle *p, int nSmooth, pqSmoothNode *nList); @@ -177,7 +186,8 @@ class PressureSmoothParams : public SmoothParams PressureSmoothParams(int _iType, int am, CSM csm, double _dTime, double _alpha, double _beta, double _dThermalDiff, double _dMetalDiff, - double dEtaCourant, double dEtaDiffusion) { + double dEtaCourant, double dEtaDiffusion, + double _dOrbFreq, Vector3D _fPeriod) { iType = _iType; activeRung = am; dTime = _dTime; @@ -195,6 +205,8 @@ class PressureSmoothParams : public SmoothParams dMetalDiffusionCoeff = _dMetalDiff; dtFacCourant = dEtaCourant*a*2.0/1.6; dtFacDiffusion = 2.0*dEtaDiffusion; + dOrbFreq = _dOrbFreq; + fPeriod = _fPeriod; } PUPable_decl(PressureSmoothParams); PressureSmoothParams(CkMigrateMessage *m) : SmoothParams(m) {} @@ -209,6 +221,8 @@ class PressureSmoothParams : public SmoothParams p|dMetalDiffusionCoeff; p|dtFacCourant; p|dtFacDiffusion; + p|dOrbFreq; + p|fPeriod; } }; diff --git a/TreePiece.cpp b/TreePiece.cpp index 3e69a044..a9dc7934 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -93,7 +93,8 @@ void TreePiece::setPeriodic(int nRepsPar, // Number of replicas in double dEwhCutPar, // Cutoff on Fourier summation int bPeriodPar, // Periodic boundaries int bComovePar, // Comoving coordinates - double dRhoFacPar // Background density + double dRhoFacPar, // Background density + double dOrbFreqPar // Orbital frequency ) { nReplicas = nRepsPar; @@ -104,6 +105,7 @@ void TreePiece::setPeriodic(int nRepsPar, // Number of replicas in bPeriodic = bPeriodPar; bComove = bComovePar; dRhoFac = dRhoFacPar; + dOrbFreq = dOrbFreqPar; if(ewt == NULL) { ewt = new EWT[nMaxEwhLoop]; } @@ -1465,6 +1467,46 @@ void TreePiece::calcEnergy(const CkCallback& cb) { #include "physconst.h" +/** +* @brief closePatch performs a velocity update as part of the kick-cross-drift-cross-kick +* algorithm described in Quinn et al. 2010 +* @param bClosing bool indicating whether the patch is opening or closing +* @param dDelta timestep +* @param p particle to update +* @param dOrbFreq orbital frequency of rotating patch +*/ +inline void closePatch(int bClosing, double dDelta, GravityParticle *p, double dOrbFreq) { +#ifdef SLIDING_PATCH +#ifndef NO_HILL + if (bClosing) { + p->velocity.x += 2.0 * dDelta * dOrbFreq * p->dPy; + p->velocity.y = p->dPy - 2 * dOrbFreq * p->position.x; + } +#endif +#endif +} +/** +* @brief openPatch performs a velocity update as part of the kick-cross-drift-cross-kick +* algorithm described in Quinn et al. 2010 +* @param bClosing bool indicating whether the patch is opening or closing +* @param dDelta timestep +* @param p particle to update +* @param dOrbFreq orbital frequency of rotating patch +*/ +inline void openPatch(int bClosing, double dDelta, GravityParticle* p, double dOrbFreq) { +#ifdef SLIDING_PATCH +#ifndef NO_HILL + if (!bClosing) { + p->dPy = p->velocity.x + 2.0 * dOrbFreq * p->position.x; + + // Cross hamiltonian + p->velocity.x += 2.0 * dDelta * dOrbFreq * p->dPy; + p->velocity.y = p->dPy - dOrbFreq * p->position.x - dOrbFreq * (p->position.x + 2.0 * dDelta * p->velocity.x); + } +#endif +#endif +} + void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], int bClosing, // Are we at the end of a timestep int bNeedVPred, // do we need to update vpred @@ -1679,8 +1721,15 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], CkAssert(p->u() < LIGHTSPEED*LIGHTSPEED/dm->Cool->dErgPerGmUnit); CkAssert(p->uPred() < LIGHTSPEED*LIGHTSPEED/dm->Cool->dErgPerGmUnit); #endif - } +#ifdef SLIDING_PATCH +#ifndef NO_HILL + p->dPyPred() = p->vPred().y + 2.0*dOrbFreq*p->position.x; +#endif +#endif + } + closePatch(bClosing, dDelta[p->rung], p, dOrbFreq); p->velocity += dDelta[p->rung]*p->treeAcceleration; + openPatch(bClosing, dDelta[p->rung], p, dOrbFreq); glassDamping(p->velocity, dDelta[p->rung], dGlassDamper); } } @@ -1994,7 +2043,7 @@ void TreePiece::assignDomain(const CkCallback &cb) myParticles[i].interMass = thisIndex; } contribute(cb); - } +} void TreePiece::drift(double dDelta, // time step in x containing // cosmo scaling @@ -2007,8 +2056,10 @@ void TreePiece::drift(double dDelta, // time step in x containing // in place bool buildTree, // is a treebuild happening before the // next drift? - double dMaxEnergy, // Maximum internal energy of gas. + double dMaxEnergy, // Maximum internal energy of gas. + double dTime, const CkCallback& cb) { + totalTime = dTime; callback = cb; // called by assignKeys() deleteTree(); @@ -2022,15 +2073,35 @@ void TreePiece::drift(double dDelta, // time step in x containing for(unsigned int i = 1; i <= myNumParticles; ++i) { GravityParticle *p = &myParticles[i]; +#ifdef SLIDING_PATCH + cosmoType fShear = 0.0; // Shear effect when crossing radial boundary + cosmoType fDx = 0.0; // Acceleration adjustment when + // crossing radial boundary +#endif if (p->iOrder >= nGrowMass) p->position += dDelta*p->velocity; if(bPeriodic) { for(int j = 0; j < 3; j++) { if(p->position[j] >= 0.5*fPeriod[j]){ p->position[j] -= fPeriod[j]; +#ifdef SLIDING_PATCH + if (j == 0) { /* radial wrap */ + fShear = 1.5 * dOrbFreq * fPeriod[0]; + p->position[1] += SHEAR(-1, dTime + dDelta, fPeriod, dOrbFreq); + fDx = -fPeriod.x; + } +#endif } + if(p->position[j] < -0.5*fPeriod[j]){ p->position[j] += fPeriod[j]; +#ifdef SLIDING_PATCH + if (j == 0) { /* radial wrap */ + fShear = -1.5 * dOrbFreq * fPeriod.x; + p->position.y += SHEAR(1, dTime + dDelta, fPeriod, dOrbFreq); + fDx = fPeriod.x; + } +#endif } bool a = (p->position[j] >= -0.5*fPeriod[j]); @@ -2047,7 +2118,14 @@ void TreePiece::drift(double dDelta, // time step in x containing } boundingBox.grow(p->position); if(bNeedVpred && TYPETest(p, TYPE_GAS)) { +#if !defined(SLIDING_PATCH) || defined(NO_HILL) p->vPred() += dvDelta*p->treeAcceleration; +#else + p->vPred().x += dvDelta*(2.0*dOrbFreq*p->dPy + p->treeAcceleration.x); + p->dPyPred() += dvDelta*p->treeAcceleration.y; + p->vPred().y = p->dPyPred() - 2.0*dOrbFreq*p->position.x; + p->vPred().z += dvDelta*p->treeAcceleration.z; +#endif glassDamping(p->vPred(), dvDelta, dGlassDamper); if(!bGasIsothermal) { #ifndef COOLING_NONE @@ -2098,7 +2176,18 @@ void TreePiece::drift(double dDelta, // time step in x containing p->fMFracIronPred() += p->fMFracIronDot()*duDelta; #endif /* DIFFUSION */ } +#ifdef SLIDING_PATCH + p->velocity.y += fShear; + p->dPy -= fShear / 3.0; /* Angular momentum is also changed. */ + if(p->isGas()) { + p->vPred().y += fShear; +#ifndef NO_HILL + p->dPyPred() -= fShear/3.0; + p->treeAcceleration.y -= dOrbFreq*dOrbFreq*fDx; +#endif } +#endif + } CkAssert(bInBox); if(!bInBox){ CkAbort("binbox2 failed\n"); @@ -5748,6 +5837,7 @@ void TreePiece::pup(PUP::er& p) { // Periodic variables p | nReplicas; p | fPeriod; + p | dOrbFreq; p | bEwald; p | fEwCut; p | dEwhCut; diff --git a/configure b/configure index 37ff1d66..eddeabbc 100755 --- a/configure +++ b/configure @@ -633,6 +633,7 @@ LIBOBJS FLAG_SANITIZER PROJECTIONS DEFAULT_LB +FLAG_SLIDING_PATCH FLAG_SIDMINTERACT FLAG_RTFORCE FLAG_VSIGVISC @@ -750,6 +751,7 @@ enable_cullenalpha enable_vsigvisc enable_rtforce enable_sidminter +enable_slidingpatch enable_default_lb enable_projections enable_opts @@ -1406,6 +1408,7 @@ Optional Features: --enable-vsigvisc alternative Monahan artificial viscosity --enable-rtforce Richie-Thomas forces --enable-sidminter SIDM interactions + --enable-slidingpatch Sliding Patch EOM --enable-default-lb default load balancer --enable-projections Charm++ projections --enable-opts DEPRECATED - Do not use @@ -4875,6 +4878,21 @@ fi +# Sliding patch (Hill) approximation + # Check whether --enable-slidingpatch was given. +if test "${enable_slidingpatch+set}" = set; then : + enableval=$enable_slidingpatch; case "$enableval" in + yes | no ) val=$enableval;; + *) as_fn_error $? "invalid argument for '--enable-slidingpatch': $enableval" "$LINENO" 5;; + esac +else + val=no +fi + + if test x$val = xyes; then FLAG_SLIDING_PATCH=-DSLIDING_PATCH; else FLAG_SLIDING_PATCH=""; fi + + + # ----------------------------------------------------------------------------- # ---- Charm++ Options -------------------------------------------------------- # ----------------------------------------------------------------------------- @@ -6924,7 +6942,7 @@ echo " "Charm path " " $CHARM_PATH echo " "Charm compiler" " $CHARMC echo " "C++ compiler" " $CXX echo -echo " "Gravity Flags " " $HEXADECAPOLE $FLAG_FLOAT $FLAG_CHANGESOFT $FLAG_DTADJUST $FLAG_TREE_BUILD $FLAG_SIDMINTERACT +echo " "Gravity Flags " " $HEXADECAPOLE $FLAG_FLOAT $FLAG_CHANGESOFT $FLAG_DTADJUST $FLAG_TREE_BUILD $FLAG_SIDMINTERACT $FLAG_SLIDING_PATCH echo " "SPH flags " " $FLAG_SPH_KERNEL $FLAG_DAMPING $FLAG_COOLING $FLAG_DIFFUSION $FLAG_FEEDBACKDIFFLIMIT $FLAG_DIFFHARMONIC $FLAG_CULLENALPHA $FLAG_VSIGVISC $FLAG_RTFORCE $FLAG_SUPERBUBBLE echo " "Misc Flags " " projections=$PROJECTIONS $FLAG_BIGKEYS sanitizer=$FLAG_SANITIZER echo " "Load balancer " " $DEFAULT_LB diff --git a/configure.ac b/configure.ac index b3daf6af..0ec72eff 100644 --- a/configure.ac +++ b/configure.ac @@ -294,6 +294,9 @@ ARG_ENABLE([rtforce], [Richie-Thomas forces], [FLAG_RTFORCE], [-DRTFORCE], [yes] # SIDM interactions ARG_ENABLE([sidminter], [SIDM interactions], [FLAG_SIDMINTERACT], [-DSIDMINTERACT], [no]) +# Sliding patch (Hill) approximation +ARG_ENABLE([slidingpatch], [Sliding Patch EOM], [FLAG_SLIDING_PATCH], [-DSLIDING_PATCH], [no]) + # ----------------------------------------------------------------------------- # ---- Charm++ Options -------------------------------------------------------- # ----------------------------------------------------------------------------- @@ -356,7 +359,7 @@ echo " "Charm path " " $CHARM_PATH echo " "Charm compiler" " $CHARMC echo " "C++ compiler" " $CXX echo -echo " "Gravity Flags " " $HEXADECAPOLE $FLAG_FLOAT $FLAG_CHANGESOFT $FLAG_DTADJUST $FLAG_TREE_BUILD $FLAG_SIDMINTERACT +echo " "Gravity Flags " " $HEXADECAPOLE $FLAG_FLOAT $FLAG_CHANGESOFT $FLAG_DTADJUST $FLAG_TREE_BUILD $FLAG_SIDMINTERACT $FLAG_SLIDING_PATCH echo " "SPH flags " " $FLAG_SPH_KERNEL $FLAG_DAMPING $FLAG_COOLING $FLAG_DIFFUSION $FLAG_FEEDBACKDIFFLIMIT $FLAG_DIFFHARMONIC $FLAG_CULLENALPHA $FLAG_VSIGVISC $FLAG_RTFORCE $FLAG_SUPERBUBBLE echo " "Misc Flags " " projections=$PROJECTIONS $FLAG_BIGKEYS sanitizer=$FLAG_SANITIZER echo " "Load balancer " " $DEFAULT_LB diff --git a/externalGravity.cpp b/externalGravity.cpp index e4133cef..06fd4d41 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -53,9 +53,13 @@ void ExternalGravity::AddParams(PRM prm) void ExternalGravity::CheckParams(PRM prm, struct parameters ¶m) { // Enable external gravity if any of the flags are set - if (bBodyForce || bPatch || bCentralBody) + if (bBodyForce || bPatch || bCentralBody) { param.bDoExternalGravity = 1; + if (bPatch) { + param.externalGravity.dOrbFreq = sqrt(param.externalGravity.dCentMass / pow(param.externalGravity.dOrbDist, 3)); + } } +} /* * @brief This function applies the external potential force to every applicable @@ -126,14 +130,19 @@ Vector3D ExternalGravity::applyPotential(GravityParticle *p) const } if (bPatch) { +#ifndef NO_HILL double r2 = dOrbDist*dOrbDist + p->position.z*p->position.z; double idt2 = dCentMass*pow(r2, -1.5); - p->treeAcceleration.z -= dCentMass*p->position.z - *pow(r2, -1.5); + p->treeAcceleration.x -= dOrbFreq*dOrbFreq*p->position.x; + // TODO: for self gravity there is a vertical frequency + // enhancement here. + p->treeAcceleration.z -= dOrbFreq*dOrbFreq*p->position.z; + p->potential += dCentMass/sqrt(r2); if(idt2 > p->dtGrav) p->dtGrav = idt2; +#endif } if(bCentralBody) { diff --git a/externalGravity.h b/externalGravity.h index 11da5e18..dc09a388 100644 --- a/externalGravity.h +++ b/externalGravity.h @@ -11,6 +11,7 @@ class ExternalGravity { int bPatch; ///< Patch in a disk double dCentMass; ///< Central mass double dOrbDist; ///< Distance of the patch from the center + double dOrbFreq; ///< Orbital frequency int bCentralBody; ///< Mass at the origin double dEqRad; ///< Equatorial radius of central body double dJ2; ///< Oblateness coefficients of central body @@ -30,6 +31,7 @@ inline void ExternalGravity::pup(PUP::er &p) { p | bPatch; p | dCentMass; p | dOrbDist; + p | dOrbFreq; p | bCentralBody; p | dEqRad; p | dJ2;