diff --git a/GravityParticle.h b/GravityParticle.h index 417730f2..41d13c37 100644 --- a/GravityParticle.h +++ b/GravityParticle.h @@ -60,7 +60,10 @@ class extraSPHData double _fESNrate; /* SN energy rate */ double _fTimeCoolIsOffUntil;/* time cooling is turned back on */ Vector3D _vPred; /* Predicted velocities for velocity - dependent forces */ + 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 */ @@ -124,6 +127,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;} @@ -186,6 +192,9 @@ class extraSPHData p | _fESNrate; p | _fTimeCoolIsOffUntil; p | _vPred; +#ifdef SLIDING_PATCH + p | _dPyPred; +#endif p | _uPred; p | _divv; p | _curlv; @@ -317,6 +326,10 @@ class GravityParticle : public ExternalGravityParticle { public: SFC::Key key; Vector3D velocity; +#ifdef SLIDING_PATCH + double dPy; ///< Canonical momentum used to update y-velocity +#endif + // inline Vector3D& vPred() { return _vPred; } Vector3D treeAcceleration; cosmoType potential; cosmoType dtGrav; ///< timestep from gravity; N.B., this @@ -364,6 +377,9 @@ class GravityParticle : public ExternalGravityParticle { ExternalGravityParticle::pup(p); p | key; p | velocity; +#ifdef SLIDING_PATCH + p | dPy; +#endif p | treeAcceleration; p | dtGrav; p | fDensity; @@ -413,6 +429,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 bb017ab4..a6498431 100644 --- a/Makefile.in +++ b/Makefile.in @@ -126,6 +126,7 @@ defines := $(strip @HEXADECAPOLE@ @FLAG_GDFORCE@ @FLAG_ARCH@ \ @FLAG_JEANSSOFTONLY@ \ @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 bf2f97fb..a5af0a51 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); @@ -440,9 +441,9 @@ mainmodule ParallelGravity { entry void liveVizDumpFrameInit(liveVizRequestMsg * msg); entry void setProjections(int bOn); #ifdef PUSH_GRAVITY - entry void buildTree(int bucketSize, const CkCallback& cb, bool merge); + entry void buildTree(int bucketSize, double dTime, const CkCallback& cb, bool merge); #else - entry void buildTree(int bucketSize, const CkCallback& cb); + entry void buildTree(int bucketSize, double dTime, const CkCallback& cb); #endif entry void startOctTreeBuild(CkReductionMsg* m); diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index d1cfcd21..4df7382b 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -1695,9 +1695,9 @@ void Main::buildTree(int iPhase) CkPrintf("Building trees ... "); double startTime = CkWallTimer(); #ifdef PUSH_GRAVITY - treeProxy.buildTree(bucketSize, CkCallbackResumeThread(),!bDoPush); + treeProxy.buildTree(bucketSize, dTime, CkCallbackResumeThread(),!bDoPush); #else - treeProxy.buildTree(bucketSize, CkCallbackResumeThread()); + treeProxy.buildTree(bucketSize, dTime, CkCallbackResumeThread()); #endif double tTB = CkWallTimer()-startTime; timings[iPhase].tTBuild += tTB; @@ -1970,7 +1970,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; @@ -2160,10 +2160,13 @@ 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 ..."); @@ -2300,8 +2303,6 @@ void Main::setupICs() { if(param.iStartStep) restartNSIDM(); } - param.externalGravity.CheckParams(prm, param); - string achLogFileName = string(param.achOutName) + ".log"; ofstream ofsLog; if(bIsRestarting) @@ -2433,6 +2434,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 " @@ -2496,7 +2503,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(); @@ -2634,7 +2641,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(); @@ -2849,7 +2856,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 @@ -3039,7 +3046,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); @@ -3613,8 +3620,8 @@ 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, - CkCallbackResumeThread()); + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, + CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3658,7 +3665,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); @@ -3681,11 +3688,10 @@ 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); } - if(param.iBinaryOut == 6) writeNCXML(achFile); } diff --git a/ParallelGravity.h b/ParallelGravity.h index b6b810c5..ee319d7f 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -1203,6 +1203,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; @@ -1330,10 +1331,9 @@ class TreePiece : public CBase_TreePiece { // called when a chunk has been used completely (chunkRemaining[chunk] == 0) void finishedChunk(int chunk); - ///Array of comparison function pointers + /// Array of comparison function pointers for ORB tree build bool (*compFuncPtr[3])(GravityParticle,GravityParticle); - double tmpTime; - double totalTime; + double dTimeTree; ///< Simulation time tree is built public: #ifdef SPCUDA @@ -1468,8 +1468,7 @@ class TreePiece : public CBase_TreePiece { #endif #endif - tmpTime=0.0; - totalTime=0.0; + dTimeTree = 0.0; // temporarely set to -1, it will updated after the tree is built numChunks=-1; prefetchRoots = NULL; @@ -1581,7 +1580,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); @@ -1706,7 +1705,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); @@ -1825,10 +1824,13 @@ class TreePiece : public CBase_TreePiece { void setProjections(int bOn); /// \brief Charm entry point to build the tree (called by Main). + /// \param bucketSize Maximum particles per bucket + /// \param dTime Simulation time tree is built + /// \param cb Callback after tree is built #ifdef PUSH_GRAVITY - void buildTree(int bucketSize, const CkCallback& cb, bool merge); + void buildTree(int bucketSize, double dTime, const CkCallback& cb, bool merge); #else - void buildTree(int bucketSize, const CkCallback& cb); + void buildTree(int bucketSize, double dTime, const CkCallback& cb); #endif /// \brief Real tree build, independent of other TreePieces. @@ -2005,15 +2007,32 @@ class TreePiece : public CBase_TreePiece { // jetley // need this in TreeWalk GenericTreeNode *getRoot() {return root;} - // need this in Compute + /// @brief adjust y (azimuthal) position if we are crossing an + /// x (radial) boundary + /// @return change in y position + inline double SHEAR(int ix, ///< Interior or exterior? + double t, ///< Current simulation time + Vector3D fPeriod, ///< patch size + double dOrbFreq) ///< Orbital frequency + { +#ifdef SLIDING_PATCH + return + (ix < 0) ? fmod(0.5 * fPeriod.y - 1.5 * ix * dOrbFreq * fPeriod.x * t, fPeriod.y) - 0.5 * fPeriod.y : + (ix > 0) ? 0.5 * fPeriod.y - fmod(0.5 * fPeriod.y + 1.5 * ix * dOrbFreq * fPeriod.x * t, fPeriod.y) : 0.0; +#else + return 0.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, dTimeTree, fPeriod.y, dOrbFreq), + z*fPeriod.z); return offset; } diff --git a/Sph.cpp b/Sph.cpp index 51d06de4..cec4d12e 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -30,7 +30,9 @@ 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, @@ -718,7 +720,9 @@ 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()); @@ -738,7 +742,9 @@ 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()); @@ -751,7 +757,9 @@ 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()); @@ -785,7 +793,9 @@ 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." @@ -1145,6 +1155,16 @@ 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 + // Detect boundary crossings and adjust shear velocity + // accordingly. + 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; @@ -1616,6 +1636,16 @@ 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 + // Detect radial boundary crossings and adjust azimuthal velocity + // accordingly. + if (params.dx.x < 0.0 && (p->position.x - q->position.x > 0.0)) { + dv.y += 1.5 * dOrbFreq * fPeriod.x; + } + else if (params.dx.x > 0.0 && (p->position.x - q->position.x < 0.0)) { + dv.y -= 1.5 * dOrbFreq * fPeriod.x; + } +#endif params.dvdotdr = vFac*dot(dv, params.dx) + fDist2*H; #ifdef GDFORCE pParams.PoverRho2 = p->PoverRho2()*p->fDensity/q->fDensity; diff --git a/Sph.h b/Sph.h index 872010a6..6a97c8f1 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); @@ -44,9 +46,12 @@ class DenDvDxSmoothParams : public SmoothParams /// @param _bStarting Simulation is starting /// @param _bHaveAlpha No need to calculate alpha /// @param _dAlphaMax Maximum SPH alpha + /// @param _dOrbFreq Angular frequency of rotating patch + /// @param _fPeriod Size of patch 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 +60,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 +83,8 @@ class DenDvDxSmoothParams : public SmoothParams p|bStarting; p|bHaveAlpha; p|dAlphaMax; + p|dOrbFreq; + p|fPeriod; } }; @@ -104,12 +113,15 @@ class DenDvDxNeighborSmParams : public DenDvDxSmoothParams /// @param dTime Current time /// @param bConstantDiffusion Fixed diffusion constant /// @param dAlphaMax Maximum SPH alpha + /// @param dOrbFreq Angular frequency of rotating patch + /// @param fPeriod Size of patch /// 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 +168,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); @@ -174,10 +188,13 @@ class PressureSmoothParams : public SmoothParams /// @param dTime Current time /// @param _alpha Artificial viscosity parameter /// @param _beta Artificial viscosity parameter + /// @param _dOrbFreq Angular frequency of rotating patch + /// @param _fPeriod Size of patch 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 +212,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 +228,8 @@ class PressureSmoothParams : public SmoothParams p|dMetalDiffusionCoeff; p|dtFacCourant; p|dtFacDiffusion; + p|dOrbFreq; + p|fPeriod; } }; diff --git a/TreePiece.cpp b/TreePiece.cpp index f13c326f..4254a11d 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,49 @@ 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 @@ -1622,7 +1667,7 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], p->fMFracIronPred() = p->fMFracIron(); #endif } - else { // predicted quantities are at the beginning + else { // predicted quantities are at the beginning // of step p->vPred() = p->velocity; if(!bGasIsothermal) { @@ -1681,8 +1726,15 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], 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); } } @@ -2008,6 +2060,7 @@ void TreePiece::drift(double dDelta, // time step in x containing bool buildTree, // is a treebuild happening before the // next drift? double dMaxEnergy, // Maximum internal energy of gas. + double dTime, // Absolute simulation time const CkCallback& cb) { callback = cb; // called by assignKeys() deleteTree(); @@ -2022,15 +2075,34 @@ 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.x; + p->position.y += 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 +2119,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 @@ -2100,7 +2179,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 + } CkMustAssert(bInBox, "binbox2 failed\n"); if(buildTree) contribute(sizeof(OrientedBox), &boundingBox, @@ -2578,12 +2668,12 @@ void TreePiece::DumpFrame(InDumpFrame in, const CkCallback& cb, int liveVizDump) */ #ifdef PUSH_GRAVITY -void TreePiece::buildTree(int bucketSize, const CkCallback& cb, bool _merge) +void TreePiece::buildTree(int bucketSize, double dTime, const CkCallback& cb, bool _merge) #else -void TreePiece::buildTree(int bucketSize, const CkCallback& cb) +void TreePiece::buildTree(int bucketSize, double dTime, const CkCallback& cb) #endif { - + dTimeTree = dTime; #if COSMO_DEBUG > 1 auto file_name = make_formatted_string("tree.%d.%d.after",thisIndex,iterationNo); char const* fout = file_name.c_str(); @@ -5742,6 +5832,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 15adbdff..f5217865 100755 --- a/configure +++ b/configure @@ -633,6 +633,7 @@ LIBOBJS FLAG_SANITIZER PROJECTIONS DEFAULT_LB +FLAG_SLIDING_PATCH FLAG_SIDMINTERACT FLAG_JEANSSOFTONLY FLAG_GDFORCE @@ -752,6 +753,7 @@ enable_vsigvisc enable_gdforce enable_jeanssoft enable_sidminter +enable_slidingpatch enable_default_lb enable_projections enable_opts @@ -1409,6 +1411,7 @@ Optional Features: --enable-gdforce Geometric Density forces --enable-jeanssoft gravitational softening jeans floor --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 @@ -4895,6 +4898,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 -------------------------------------------------------- # ----------------------------------------------------------------------------- @@ -6944,7 +6962,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_GDFORCE $FLAG_JEANSSOFTONLY $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 aee0baae..ef4407fb 100644 --- a/configure.ac +++ b/configure.ac @@ -299,6 +299,9 @@ ARG_ENABLE([jeanssoft], [gravitational softening jeans floor], [FLAG_JEANSSOFTON # 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 -------------------------------------------------------- # ----------------------------------------------------------------------------- @@ -361,7 +364,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_GDFORCE $FLAG_JEANSSOFTONLY $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..a2241b86 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -53,8 +53,12 @@ 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)); + } + } } /* @@ -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;