From ead977b4079e1d135fd416bfee98dd6a4938ba4f Mon Sep 17 00:00:00 2001 From: brianxchen Date: Fri, 29 Jan 2021 10:30:06 -0800 Subject: [PATCH 01/20] undid stylistic changes, adding shear --- GravityParticle.h | 23 +++++++++++++---------- ParallelGravity.ci | 3 ++- ParallelGravity.cpp | 5 ++++- ParallelGravity.h | 2 +- TreePiece.cpp | 44 ++++++++++++++++++++++++++++++++++---------- externalGravity.cpp | 2 -- 6 files changed, 54 insertions(+), 25 deletions(-) diff --git a/GravityParticle.h b/GravityParticle.h index 4d870bfb..05bf6b7b 100644 --- a/GravityParticle.h +++ b/GravityParticle.h @@ -59,7 +59,8 @@ class extraSPHData double _fMFracIron; /* Iron mass fraction */ double _fESNrate; /* SN energy rate */ double _fTimeCoolIsOffUntil;/* time cooling is turned back on */ - Vector3D _vPred; /* Predicted velocities for velocity dependent forces */ + Vector3D _vPred; /* Predicted velocities for velocity + dependent forces */ double _uPred; /* Predicted internal energy */ double _divv; /* Diverence of the velocity */ Vector3D _curlv; /* Curl of the velocity */ @@ -120,7 +121,7 @@ class extraSPHData inline double& fMFracIron() {return _fMFracIron;} inline double& fESNrate() {return _fESNrate;} inline double& fTimeCoolIsOffUntil() {return _fTimeCoolIsOffUntil;} - // inline Vector3D& vPred() {return _vPred;} + inline Vector3D& vPred() {return _vPred;} inline double& uPred() {return _uPred;} inline double& divv() {return _divv;} inline Vector3D& curlv() {return _curlv;} @@ -180,7 +181,7 @@ class extraSPHData p | _fMFracOxygen; p | _fESNrate; p | _fTimeCoolIsOffUntil; - // p | _vPred; + p | _vPred; p | _uPred; p | _divv; p | _curlv; @@ -310,9 +311,11 @@ class GravityParticle : public ExternalGravityParticle { public: SFC::Key key; Vector3D velocity; - double dPy; +#ifdef SLIDING_PATCH + double dPy; ///< Canonical momentum used to update y-velocity +#endif Vector3D _vPred; - inline Vector3D& vPred() { return _vPred;} + // inline Vector3D& vPred() { return _vPred; } Vector3D treeAcceleration; cosmoType potential; cosmoType dtGrav; ///< timestep from gravity; N.B., this @@ -360,9 +363,11 @@ class GravityParticle : public ExternalGravityParticle { ExternalGravityParticle::pup(p); p | key; p | velocity; +#ifdef SLIDING_PATCH p | dPy; - p | _vPred; +#endif p | treeAcceleration; + p | _vPred; p | dtGrav; p | fDensity; p | fBall; @@ -410,7 +415,7 @@ 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());} - // inline Vector3D& vPred() { IMAGAS; return (((extraSPHData*)extraData)->vPred());} + inline Vector3D& vPred() { IMAGAS; return (((extraSPHData*)extraData)->vPred());} inline double& uPred() {IMAGAS; return (((extraSPHData*)extraData)->uPred());} inline double& divv() { IMAGAS; return (((extraSPHData*)extraData)->divv());} inline Vector3D& curlv() { IMAGAS; return (((extraSPHData*)extraData)->curlv());} @@ -636,7 +641,6 @@ class ExternalSmoothParticle { iType = p->iType; rung = p->rung; treeAcceleration = p->treeAcceleration; - vPred = p->vPred(); if(TYPETest(p, TYPE_GAS)) { mumax = p->mumax(); PdV = p->PdV(); @@ -703,9 +707,8 @@ class ExternalSmoothParticle { tmp->iType = iType; tmp->rung = rung; tmp->treeAcceleration = treeAcceleration; - tmp->vPred() = vPred; if(TYPETest(tmp, TYPE_GAS)) { - tmp->mumax() = mumax; + tmp->mumax() = mumax; tmp->PdV() = PdV; tmp->c() = c; tmp->PoverRho2() = PoverRho2; diff --git a/ParallelGravity.ci b/ParallelGravity.ci index 3692bd5f..fc7ee85a 100644 --- a/ParallelGravity.ci +++ b/ParallelGravity.ci @@ -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 dOrbFreq, 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 46e82054..273d011e 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -1877,7 +1877,10 @@ void Main::kick(bool bClosing, int iActiveRung, int nextMaxRung, double a = csmTime2Exp(param.csm,dTime); double startTime = CkWallTimer(); - double dOrbFreq = sqrt(param.ExternalGravity.dCentMass / pow(param.ExternalGravity.dOrbDist, 3)); + double dOrbFreq = sqrt(param.ExternalGravity.dCentMass / + pow(param.ExternalGravity.dOrbDist, 3)); // Orbital frequency + // sqrt(G*M/r^3) + // with G=1 treeProxy.kick(iActiveRung, dKickFac, bClosing, param.bDoGas, param.bGasIsothermal, dOrbFreq, param.dMaxEnergy, duKick, (param.dConstGamma-1), param.dThermalCondSatCoeff/a, diff --git a/ParallelGravity.h b/ParallelGravity.h index a24306f3..251743af 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -1695,7 +1695,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 dOrbFreq, double dTime, const CkCallback& cb); void initAccel(int iKickRung, const CkCallback& cb); #ifdef COOLING_MOLECULARH void distribLymanWerner(const CkCallback& cb); diff --git a/TreePiece.cpp b/TreePiece.cpp index ef034902..72c6b919 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1469,7 +1469,7 @@ 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 int bGasIsothermal, // Isothermal EOS - double dOrbFreq, // Orbital frequency + double dOrbFreq, ///< Orbital frequency double dMaxEnergy, // Maximum internal energy of gas. double duDelta[MAXRUNG+1], // dts for energy double gammam1, // Adiabatic index - 1 @@ -1482,12 +1482,9 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], for(unsigned int i = 1; i <= myNumParticles; ++i) { GravityParticle *p = &myParticles[i]; if(p->rung >= iKickRung) { - if(bNeedVPred && TYPETest(p, TYPE_GAS)) { if(bClosing) { // update predicted quantities to end of step - p->vPred() = p->velocity + dDelta[p->rung]*p->treeAcceleration; - glassDamping(p->vPred(), dDelta[p->rung], dGlassDamper); if(!bGasIsothermal) { #ifndef COOLING_NONE @@ -1617,12 +1614,9 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], p->fMFracIronPred() = p->fMFracIron(); #endif } - if(!bClosing) { // predicted quantities are at the beginning // of step p->vPred() = p->velocity; - - if(!bGasIsothermal) { p->uPred() = p->u(); #ifndef COOLING_NONE @@ -2001,6 +1995,11 @@ void TreePiece::assignDomain(const CkCallback &cb) contribute(cb); } +inline double SHEAR(int ix, double t, Vector3D fPeriod, double dOrbFreq, GravityParticle *p) { + ((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); +} + void TreePiece::drift(double dDelta, // time step in x containing // cosmo scaling int bNeedVpred, // Update predicted velocities @@ -2012,7 +2011,9 @@ 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 dOrbFreq, // Orbital frequency + double dTime, const CkCallback& cb) { callback = cb; // called by assignKeys() deleteTree(); @@ -2027,15 +2028,32 @@ 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 + FLOAT fShear = 0.0; + p->bAzWrap = 0; +#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 * p->dOrbFreq * p->fPeriod[0]; + p->position[1] += SHEAR(-1, dTime + dDelta, fPeriod, dOrbFreq, *p); + } +#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 * p->dOrbFreq * p->fPeriod[0]; + p->position[1] += SHEAR(1, dTime + dDelta, fPeriod, dOrbFreq, *p); + } +#endif } bool a = (p->position[j] >= -0.5*fPeriod[j]); @@ -2051,9 +2069,9 @@ void TreePiece::drift(double dDelta, // time step in x containing } } boundingBox.grow(p->position); - p->vPred() += dvDelta * p->treeAcceleration; if(bNeedVpred && TYPETest(p, TYPE_GAS)) { - glassDamping(p->vPred(), dvDelta, dGlassDamper); + p->vPred() += p->treeAcceleration * dvDelta; + glassDamping(p->vPred(), dvDelta, dGlassDamper); if(!bGasIsothermal) { #ifndef COOLING_NONE p->uPred() += p->uDot()*duDelta; @@ -2104,6 +2122,12 @@ void TreePiece::drift(double dDelta, // time step in x containing #endif /* DIFFUSION */ } } +#ifdef SLIDING_PATCH + p->velocity[1] += fShear; + p->dPy -= fShear / 3.0; /* Angular momentum is + also changed. */ +#endif + CkAssert(bInBox); if(!bInBox){ CkAbort("binbox2 failed\n"); diff --git a/externalGravity.cpp b/externalGravity.cpp index 841085ff..72bec35e 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -55,8 +55,6 @@ void ExternalGravity::CheckParams(PRM prm, struct parameters ¶m) // Enable external gravity if any of the flags are set if (bBodyForce || bPatch || bCentralBody) param.bDoExternalGravity = 1; - - // Calculate dOrbFreq } /* From 8cc79045c40323f96dc1adcef8a40ba9e53bb12c Mon Sep 17 00:00:00 2001 From: brianxchen Date: Fri, 5 Feb 2021 09:36:42 -0800 Subject: [PATCH 02/20] Working shear --- GravityParticle.h | 8 +++---- ParallelGravity.cpp | 20 +++++++++-------- TreePiece.cpp | 55 +++++++++++++++++++++++++++------------------ externalGravity.cpp | 2 +- 4 files changed, 49 insertions(+), 36 deletions(-) diff --git a/GravityParticle.h b/GravityParticle.h index 05bf6b7b..43aa3545 100644 --- a/GravityParticle.h +++ b/GravityParticle.h @@ -59,8 +59,8 @@ class extraSPHData double _fMFracIron; /* Iron mass fraction */ double _fESNrate; /* SN energy rate */ double _fTimeCoolIsOffUntil;/* time cooling is turned back on */ - Vector3D _vPred; /* Predicted velocities for velocity - dependent forces */ + Vector3D _vPred; /* Predicted velocities for velocity + dependent forces */ double _uPred; /* Predicted internal energy */ double _divv; /* Diverence of the velocity */ Vector3D _curlv; /* Curl of the velocity */ @@ -314,7 +314,6 @@ class GravityParticle : public ExternalGravityParticle { #ifdef SLIDING_PATCH double dPy; ///< Canonical momentum used to update y-velocity #endif - Vector3D _vPred; // inline Vector3D& vPred() { return _vPred; } Vector3D treeAcceleration; cosmoType potential; @@ -367,7 +366,6 @@ class GravityParticle : public ExternalGravityParticle { p | dPy; #endif p | treeAcceleration; - p | _vPred; p | dtGrav; p | fDensity; p | fBall; @@ -642,6 +640,7 @@ class ExternalSmoothParticle { rung = p->rung; treeAcceleration = p->treeAcceleration; if(TYPETest(p, TYPE_GAS)) { + vPred = p->vPred(); mumax = p->mumax(); PdV = p->PdV(); c = p->c(); @@ -708,6 +707,7 @@ class ExternalSmoothParticle { tmp->rung = rung; tmp->treeAcceleration = treeAcceleration; if(TYPETest(tmp, TYPE_GAS)) { + tmp->vPred() = vPred; tmp->mumax() = mumax; tmp->PdV() = PdV; tmp->c() = c; diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index 273d011e..802ec9b5 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -1877,8 +1877,8 @@ void Main::kick(bool bClosing, int iActiveRung, int nextMaxRung, double a = csmTime2Exp(param.csm,dTime); double startTime = CkWallTimer(); - double dOrbFreq = sqrt(param.ExternalGravity.dCentMass / - pow(param.ExternalGravity.dOrbDist, 3)); // Orbital frequency + double dOrbFreq = sqrt(param.externalGravity.dCentMass / + pow(param.externalGravity.dOrbDist, 3)); // Orbital frequency // sqrt(G*M/r^3) // with G=1 treeProxy.kick(iActiveRung, dKickFac, bClosing, param.bDoGas, @@ -1973,9 +1973,11 @@ void Main::advanceBigStep(int iStep) { double dDriftFac = csmComoveDriftFac(param.csm, dTime, dTimeSub); double dKickFac = csmComoveKickFac(param.csm, dTime, dTimeSub); bool bBuildTree = (iSub + 1 == driftSteps); + double dOrbFreq = sqrt(param.externalGravity.dCentMass / + pow(param.externalGravity.dOrbDist, 3)); treeProxy.drift(dDriftFac, param.bDoGas, param.bGasIsothermal, dKickFac, dTimeSub, nGrowMassDrift, bBuildTree, - param.dMaxEnergy, + param.dMaxEnergy, dOrbFreq, dTime, CkCallbackResumeThread()); double tDrift = CkWallTimer() - startTime; timings[activeRung].tDrift += tDrift; @@ -2491,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, 0, 0, CkCallbackResumeThread()); initialForces(); @@ -2629,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, 0, 0, CkCallbackResumeThread()); if(param.bGasCooling || param.bStarForm) initCooling(); @@ -2844,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, 0, 0, CkCallbackResumeThread()); treeProxy[0].flushStarLog(CkCallbackResumeThread()); param.iStartStep = iStep; // update so that restart continues on @@ -3037,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, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3610,7 +3612,7 @@ void Main::writeOutput(int iStep) << endl; // 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, 0, CkCallbackResumeThread()); domainDecomp(0); if(param.nSteps != 0 && param.bDoDensity) { @@ -3656,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, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); diff --git a/TreePiece.cpp b/TreePiece.cpp index 72c6b919..4ba01fb6 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1465,6 +1465,22 @@ void TreePiece::calcEnergy(const CkCallback& cb) { #include "physconst.h" +inline void updatePatch(int bClosing, double dDelta, GravityParticle& p, double dOrbFreq) { +#ifdef SLIDING_PATCH + if (bClosing) { + p->velocity[0] += 2.0 * dDelta[p->rung] * dOrbFreq * p->dPy; + p->velocity[1] = p->dPy - 2 * dOrbFreq * p->position[0]; + } + else { + p->dPy = p->velocity[1] + 2.0 * dOrbFreq * p->position[0]; + + // Cross hamiltonian + p->velocity[0] += 2.0 * dDelta[p->rung] * dOrbFreq * p->dPy; + p->velocity[1] = p->dPy - dOrbFreq * p->position[0] - dOrbFreq * (p->position[0] + 2.0 * dDelta[p->rung] * p->velocity[0]); + } +#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 @@ -1484,7 +1500,8 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], if(p->rung >= iKickRung) { if(bNeedVPred && TYPETest(p, TYPE_GAS)) { if(bClosing) { // update predicted quantities to end of step - p->vPred() = p->velocity + dDelta[p->rung]*p->treeAcceleration; + p->vPred() = p->velocity + + dDelta[p->rung]*p->treeAcceleration; glassDamping(p->vPred(), dDelta[p->rung], dGlassDamper); if(!bGasIsothermal) { #ifndef COOLING_NONE @@ -1614,7 +1631,7 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], p->fMFracIronPred() = p->fMFracIron(); #endif } - if(!bClosing) { // predicted quantities are at the beginning + else { // predicted quantities are at the beginning // of step p->vPred() = p->velocity; if(!bGasIsothermal) { @@ -1667,18 +1684,9 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], CkAssert(p->uPred() < LIGHTSPEED*LIGHTSPEED/dm->Cool->dErgPerGmUnit); #endif } - if (bClosing) { - p->velocity[0] += 2.0 * dDelta[p->rung] * dOrbFreq * p->dPy; - p->velocity[1] = p->dPy - 2 * dOrbFreq * p->position[0]; - } + updatePatch(bClosing, dDelta[p->rung], *p, dOrbFreq); p->velocity += dDelta[p->rung]*p->treeAcceleration; - if (!bClosing) { - p->dPy = p->velocity[1] + 2.0 * dOrbFreq * p->position[0]; - - // Cross hamiltonian - p->velocity[0] += 2.0 * dDelta[p->rung] * dOrbFreq * p->dPy; - p->velocity[1] = p->dPy - dOrbFreq * p->position[0] - dOrbFreq * (p->position[0] + 2.0 * dDelta[p->rung] * p->velocity[0]); - } + updatePatch(bClosing, dDelta[p->rung], *p, dOrbFreq); glassDamping(p->velocity, dDelta[p->rung], dGlassDamper); } @@ -1993,11 +2001,15 @@ void TreePiece::assignDomain(const CkCallback &cb) myParticles[i].interMass = thisIndex; } contribute(cb); - } +} -inline double SHEAR(int ix, double t, Vector3D fPeriod, double dOrbFreq, GravityParticle *p) { - ((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); +inline double SHEAR(int ix, ///< Interior or exterior? + double t, ///< time, in dTime+dDelta + Vector3D fPeriod, ///< vector for dxPeriod and dyPeriod + double dOrbFreq) ///< Orbital frequency +{ + 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; } void TreePiece::drift(double dDelta, // time step in x containing @@ -2012,7 +2024,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 dOrbFreq, // Orbital frequency + double dOrbFreq, ///< Orbital frequency of patch double dTime, const CkCallback& cb) { callback = cb; // called by assignKeys() @@ -2029,8 +2041,7 @@ 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 - FLOAT fShear = 0.0; - p->bAzWrap = 0; + cosmoType fShear = 0.0; #endif if (p->iOrder >= nGrowMass) p->position += dDelta*p->velocity; @@ -2070,8 +2081,8 @@ void TreePiece::drift(double dDelta, // time step in x containing } boundingBox.grow(p->position); if(bNeedVpred && TYPETest(p, TYPE_GAS)) { - p->vPred() += p->treeAcceleration * dvDelta; - glassDamping(p->vPred(), dvDelta, dGlassDamper); + p->vPred() += dvDelta * p->treeAcceleration; + glassDamping(p->vPred(), dvDelta, dGlassDamper); if(!bGasIsothermal) { #ifndef COOLING_NONE p->uPred() += p->uDot()*duDelta; diff --git a/externalGravity.cpp b/externalGravity.cpp index 72bec35e..2144a35b 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -55,7 +55,7 @@ void ExternalGravity::CheckParams(PRM prm, struct parameters ¶m) // Enable external gravity if any of the flags are set if (bBodyForce || bPatch || bCentralBody) param.bDoExternalGravity = 1; -} + } /* * @brief This function applies the external potential force to every applicable From fdbb78fe124df3d13943b21ddd7009a42ea49f6a Mon Sep 17 00:00:00 2001 From: brianxchen Date: Fri, 5 Feb 2021 11:25:24 -0800 Subject: [PATCH 03/20] Actual working shear update --- TreePiece.cpp | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/TreePiece.cpp b/TreePiece.cpp index 4ba01fb6..f4f0b6aa 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1465,18 +1465,18 @@ void TreePiece::calcEnergy(const CkCallback& cb) { #include "physconst.h" -inline void updatePatch(int bClosing, double dDelta, GravityParticle& p, double dOrbFreq) { +inline void updatePatch(int bClosing, double dDelta, GravityParticle *p, double dOrbFreq) { #ifdef SLIDING_PATCH if (bClosing) { - p->velocity[0] += 2.0 * dDelta[p->rung] * dOrbFreq * p->dPy; + p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; p->velocity[1] = p->dPy - 2 * dOrbFreq * p->position[0]; } else { p->dPy = p->velocity[1] + 2.0 * dOrbFreq * p->position[0]; // Cross hamiltonian - p->velocity[0] += 2.0 * dDelta[p->rung] * dOrbFreq * p->dPy; - p->velocity[1] = p->dPy - dOrbFreq * p->position[0] - dOrbFreq * (p->position[0] + 2.0 * dDelta[p->rung] * p->velocity[0]); + p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; + p->velocity[1] = p->dPy - dOrbFreq * p->position[0] - dOrbFreq * (p->position[0] + 2.0 * dDelta * p->velocity[0]); } #endif } @@ -1684,9 +1684,9 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], CkAssert(p->uPred() < LIGHTSPEED*LIGHTSPEED/dm->Cool->dErgPerGmUnit); #endif } - updatePatch(bClosing, dDelta[p->rung], *p, dOrbFreq); + updatePatch(bClosing, dDelta[p->rung], p, dOrbFreq); p->velocity += dDelta[p->rung]*p->treeAcceleration; - updatePatch(bClosing, dDelta[p->rung], *p, dOrbFreq); + updatePatch(bClosing, dDelta[p->rung], p, dOrbFreq); glassDamping(p->velocity, dDelta[p->rung], dGlassDamper); } @@ -2051,8 +2051,8 @@ void TreePiece::drift(double dDelta, // time step in x containing p->position[j] -= fPeriod[j]; #ifdef SLIDING_PATCH if (j == 0) { /* radial wrap */ - fShear = 1.5 * p->dOrbFreq * p->fPeriod[0]; - p->position[1] += SHEAR(-1, dTime + dDelta, fPeriod, dOrbFreq, *p); + fShear = 1.5 * dOrbFreq * fPeriod[0]; + p->position[1] += SHEAR(-1, dTime + dDelta, fPeriod, dOrbFreq); } #endif } @@ -2061,8 +2061,8 @@ void TreePiece::drift(double dDelta, // time step in x containing p->position[j] += fPeriod[j]; #ifdef SLIDING_PATCH if (j == 0) { /* radial wrap */ - fShear = -1.5 * p->dOrbFreq * p->fPeriod[0]; - p->position[1] += SHEAR(1, dTime + dDelta, fPeriod, dOrbFreq, *p); + fShear = -1.5 * dOrbFreq * fPeriod[0]; + p->position[1] += SHEAR(1, dTime + dDelta, fPeriod, dOrbFreq); } #endif } @@ -2132,12 +2132,13 @@ void TreePiece::drift(double dDelta, // time step in x containing p->fMFracIronPred() += p->fMFracIronDot()*duDelta; #endif /* DIFFUSION */ } - } #ifdef SLIDING_PATCH - p->velocity[1] += fShear; - p->dPy -= fShear / 3.0; /* Angular momentum is - also changed. */ + p->velocity[1] += fShear; + p->dPy -= fShear / 3.0; /* Angular momentum is + also changed. */ #endif + } + CkAssert(bInBox); if(!bInBox){ From 91b54b0ce38822eff4f6aed111cc1b24a6ea3372 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Thu, 11 Feb 2021 10:01:41 -0800 Subject: [PATCH 04/20] Split updatePatch into two inlines - now consistent with pkdgrav --- TreePiece.cpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/TreePiece.cpp b/TreePiece.cpp index f4f0b6aa..42a4158e 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1465,17 +1465,21 @@ void TreePiece::calcEnergy(const CkCallback& cb) { #include "physconst.h" -inline void updatePatch(int bClosing, double dDelta, GravityParticle *p, double dOrbFreq) { +inline void openPatch(int bClosing, double dDelta, GravityParticle *p, double dOrbFreq) { #ifdef SLIDING_PATCH if (bClosing) { p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; p->velocity[1] = p->dPy - 2 * dOrbFreq * p->position[0]; } - else { +#endif +} +inline void closePatch(int bClosing, double dDelta, GravityParticle* p, double dOrbFreq) { +#ifdef SLIDING_PATCH + if (!bClosing) { p->dPy = p->velocity[1] + 2.0 * dOrbFreq * p->position[0]; // Cross hamiltonian - p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; + p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; p->velocity[1] = p->dPy - dOrbFreq * p->position[0] - dOrbFreq * (p->position[0] + 2.0 * dDelta * p->velocity[0]); } #endif @@ -1684,9 +1688,9 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], CkAssert(p->uPred() < LIGHTSPEED*LIGHTSPEED/dm->Cool->dErgPerGmUnit); #endif } - updatePatch(bClosing, dDelta[p->rung], p, dOrbFreq); + openPatch(bClosing, dDelta[p->rung], p, dOrbFreq); p->velocity += dDelta[p->rung]*p->treeAcceleration; - updatePatch(bClosing, dDelta[p->rung], p, dOrbFreq); + closePatch(bClosing, dDelta[p->rung], p, dOrbFreq); glassDamping(p->velocity, dDelta[p->rung], dGlassDamper); } @@ -2138,8 +2142,6 @@ void TreePiece::drift(double dDelta, // time step in x containing also changed. */ #endif } - - CkAssert(bInBox); if(!bInBox){ CkAbort("binbox2 failed\n"); From ff3416f96f0f1f4bb2d04853fa034d13170c1c22 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Fri, 26 Feb 2021 10:13:15 -0800 Subject: [PATCH 05/20] changes to shearing patch --- ParallelGravity.h | 12 +++++++++++- Sph.cpp | 24 ++++++++++++++++++++---- Sph.h | 18 ++++++++++++++++-- TreePiece.cpp | 17 ++++------------- parameters.h | 2 ++ 5 files changed, 53 insertions(+), 20 deletions(-) diff --git a/ParallelGravity.h b/ParallelGravity.h index 251743af..e37e09be 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; int nReplicas; int bEwald; /* Perform Ewald */ double fEwCut; @@ -1995,13 +1996,22 @@ 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, ///< time, in dTime+dDelta + Vector3D fPeriod, ///< vector for dxPeriod and dyPeriod + double dOrbFreq) ///< Orbital frequency + { + 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; + } + 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[1], dOrbFreq), z*fPeriod.z); return offset; } diff --git a/Sph.cpp b/Sph.cpp index efc723f7..7cc49c42 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.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.dOrbFreq, param.fPeriod); double startTime = CkWallTimer(); treeProxy.startSmooth(&pDen, 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.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.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, dvzdx += dvz*dx*rs1; dvzdy += dvz*dy*rs1; dvzdz += dvz*dz*rs1; +#ifdef SLIDING_PATCH + if (dx < 0.0 && (p->position[0] - q->position[0] > 0.0)) { + dvy += 1.5 * dOrbFreq * fPeriod[0]; + } + else if (dx > 0.0 && (p->position[0] - q->position[0] < 0.0)) { + dvy -= 1.5 * dOrbFreq * fPeriod[0]; + } +#endif divvnorm += (dx*dx+dy*dy+dz*dz)*rs1; /* Grad P estimate */ /* This used to be: @@ -1562,6 +1570,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[0] - q->position[0] > 0.0)) { + dv[1] += 1.5 * dOrbFreq * fPeriod[0]; + } + else if (params.dx.x > 0.0 && (p->position[0] - q->position[0] < 0.0)) { + dv[1] -= 1.5 * dOrbFreq * fPeriod[0]; + } +#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..9efd87ff 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; } }; @@ -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 42a4158e..83d2c522 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1465,7 +1465,7 @@ void TreePiece::calcEnergy(const CkCallback& cb) { #include "physconst.h" -inline void openPatch(int bClosing, double dDelta, GravityParticle *p, double dOrbFreq) { +inline void closePatch(int bClosing, double dDelta, GravityParticle *p, double dOrbFreq) { #ifdef SLIDING_PATCH if (bClosing) { p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; @@ -1473,7 +1473,7 @@ inline void openPatch(int bClosing, double dDelta, GravityParticle *p, double dO } #endif } -inline void closePatch(int bClosing, double dDelta, GravityParticle* p, double dOrbFreq) { +inline void openPatch(int bClosing, double dDelta, GravityParticle* p, double dOrbFreq) { #ifdef SLIDING_PATCH if (!bClosing) { p->dPy = p->velocity[1] + 2.0 * dOrbFreq * p->position[0]; @@ -1688,9 +1688,9 @@ void TreePiece::kick(int iKickRung, double dDelta[MAXRUNG+1], CkAssert(p->uPred() < LIGHTSPEED*LIGHTSPEED/dm->Cool->dErgPerGmUnit); #endif } - openPatch(bClosing, dDelta[p->rung], p, dOrbFreq); - p->velocity += dDelta[p->rung]*p->treeAcceleration; 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); } @@ -2007,15 +2007,6 @@ void TreePiece::assignDomain(const CkCallback &cb) contribute(cb); } -inline double SHEAR(int ix, ///< Interior or exterior? - double t, ///< time, in dTime+dDelta - Vector3D fPeriod, ///< vector for dxPeriod and dyPeriod - double dOrbFreq) ///< Orbital frequency -{ - 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; -} - void TreePiece::drift(double dDelta, // time step in x containing // cosmo scaling int bNeedVpred, // Update predicted velocities diff --git a/parameters.h b/parameters.h index 7f8c6da2..30adeea6 100644 --- a/parameters.h +++ b/parameters.h @@ -47,6 +47,7 @@ typedef struct parameters { int iOrder; int bConcurrentSph; double dFracNoDomainDecomp; + double dOrbFreq; #ifdef PUSH_GRAVITY double dFracPushParticles; #endif @@ -193,6 +194,7 @@ inline void operator|(PUP::er &p, Parameters ¶m) { p|param.iOrder; p|param.bConcurrentSph; p|param.dFracNoDomainDecomp; + p|param.dOrbFreq; #ifdef PUSH_GRAVITY p|param.dFracPushParticles; #endif From 47155b897d55e6d4c273f1e7565ebd66c6f20bb0 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Tue, 2 Mar 2021 14:20:24 -0800 Subject: [PATCH 06/20] Shearing patch with too wide spread --- Sph.cpp | 8 ++++---- externalGravity.cpp | 1 + externalGravity.h | 2 ++ parameters.h | 2 -- 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/Sph.cpp b/Sph.cpp index 7cc49c42..5d2583be 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.dOrbFreq, param.fPeriod); + 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.dOrbFreq, param.fPeriod); + param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod); double startTime = CkWallTimer(); treeProxy.startSmooth(&pDen, 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.dOrbFreq, param.fPeriod); + 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.dOrbFreq, param.fPeriod); + param.dEtaCourant, param.dEtaDiffusion, param.externalGravity.dOrbFreq, param.fPeriod); double startTime = CkWallTimer(); treeProxy.startReSmooth(&pPressure, CkCallbackResumeThread()); ckout << " took " << (CkWallTimer() - startTime) << " seconds." diff --git a/externalGravity.cpp b/externalGravity.cpp index 2144a35b..3aa88643 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -55,6 +55,7 @@ void ExternalGravity::CheckParams(PRM prm, struct parameters ¶m) // Enable external gravity if any of the flags are set if (bBodyForce || bPatch || bCentralBody) param.bDoExternalGravity = 1; + param.externalGravity.dOrbFreq = sqrt(param.externalGravity.dCentMass / pow(param.externalGravity.dOrbDist, 3)); } /* 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; diff --git a/parameters.h b/parameters.h index 30adeea6..7f8c6da2 100644 --- a/parameters.h +++ b/parameters.h @@ -47,7 +47,6 @@ typedef struct parameters { int iOrder; int bConcurrentSph; double dFracNoDomainDecomp; - double dOrbFreq; #ifdef PUSH_GRAVITY double dFracPushParticles; #endif @@ -194,7 +193,6 @@ inline void operator|(PUP::er &p, Parameters ¶m) { p|param.iOrder; p|param.bConcurrentSph; p|param.dFracNoDomainDecomp; - p|param.dOrbFreq; #ifdef PUSH_GRAVITY p|param.dFracPushParticles; #endif From fda28940451b08662ba4e3178e7e4c770c4ecc5a Mon Sep 17 00:00:00 2001 From: brianxchen Date: Fri, 5 Mar 2021 10:40:59 -0800 Subject: [PATCH 07/20] Working shearing patch - tends to ideal line after a long time --- Sph.cpp | 2 +- Sph.h | 4 ++-- externalGravity.cpp | 3 ++- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Sph.cpp b/Sph.cpp index 5d2583be..564af778 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -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()); diff --git a/Sph.h b/Sph.h index 9efd87ff..41e11eaf 100644 --- a/Sph.h +++ b/Sph.h @@ -114,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) { diff --git a/externalGravity.cpp b/externalGravity.cpp index 3aa88643..f941fdad 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -53,10 +53,11 @@ 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; param.externalGravity.dOrbFreq = sqrt(param.externalGravity.dCentMass / pow(param.externalGravity.dOrbDist, 3)); } + } /* * @brief This function applies the external potential force to every applicable From 1b68a18d95918d282ce6ab4d4769648a62216b34 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Fri, 19 Mar 2021 10:31:57 -0700 Subject: [PATCH 08/20] moved velocity update earlier --- Sph.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/Sph.cpp b/Sph.cpp index 564af778..3b4d1eb6 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -1110,6 +1110,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[0] - q->position[0] > 0.0)) { + dvy += 1.5 * dOrbFreq * fPeriod[0]; + } + else if (dx > 0.0 && (p->position[0] - q->position[0] < 0.0)) { + dvy -= 1.5 * dOrbFreq * fPeriod[0]; + } +#endif dvxdx += dvx*dx*rs1; dvxdy += dvx*dy*rs1; dvxdz += dvx*dz*rs1; @@ -1119,14 +1127,7 @@ void DenDvDxSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, dvzdx += dvz*dx*rs1; dvzdy += dvz*dy*rs1; dvzdz += dvz*dz*rs1; -#ifdef SLIDING_PATCH - if (dx < 0.0 && (p->position[0] - q->position[0] > 0.0)) { - dvy += 1.5 * dOrbFreq * fPeriod[0]; - } - else if (dx > 0.0 && (p->position[0] - q->position[0] < 0.0)) { - dvy -= 1.5 * dOrbFreq * fPeriod[0]; - } -#endif + divvnorm += (dx*dx+dy*dy+dz*dz)*rs1; /* Grad P estimate */ /* This used to be: From 8ed66ee1a2a2c7cf77e883234a7264e3df64722a Mon Sep 17 00:00:00 2001 From: brianxchen Date: Mon, 3 May 2021 15:05:04 -0700 Subject: [PATCH 09/20] added orbital frequency to treepiece calls, shears correctly now --- ParallelGravity.ci | 2 +- ParallelGravity.cpp | 35 +++++++++++++++++++---------------- ParallelGravity.h | 2 +- Sph.cpp | 4 ++-- TreePiece.cpp | 9 +++++++-- 5 files changed, 30 insertions(+), 22 deletions(-) diff --git a/ParallelGravity.ci b/ParallelGravity.ci index fc7ee85a..921ed9a6 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); diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index 802ec9b5..dee11611 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -1877,12 +1877,8 @@ void Main::kick(bool bClosing, int iActiveRung, int nextMaxRung, double a = csmTime2Exp(param.csm,dTime); double startTime = CkWallTimer(); - double dOrbFreq = sqrt(param.externalGravity.dCentMass / - pow(param.externalGravity.dOrbDist, 3)); // Orbital frequency - // sqrt(G*M/r^3) - // with G=1 treeProxy.kick(iActiveRung, dKickFac, bClosing, param.bDoGas, - param.bGasIsothermal, dOrbFreq, param.dMaxEnergy, duKick, + param.bGasIsothermal, param.externalGravity.dOrbFreq, param.dMaxEnergy, duKick, (param.dConstGamma-1), param.dThermalCondSatCoeff/a, param.feedback->dMultiPhaseMaxTime, param.feedback->dMultiPhaseMinTemp, @@ -1973,11 +1969,9 @@ void Main::advanceBigStep(int iStep) { double dDriftFac = csmComoveDriftFac(param.csm, dTime, dTimeSub); double dKickFac = csmComoveKickFac(param.csm, dTime, dTimeSub); bool bBuildTree = (iSub + 1 == driftSteps); - double dOrbFreq = sqrt(param.externalGravity.dCentMass / - pow(param.externalGravity.dOrbDist, 3)); treeProxy.drift(dDriftFac, param.bDoGas, param.bGasIsothermal, dKickFac, dTimeSub, nGrowMassDrift, bBuildTree, - param.dMaxEnergy, dOrbFreq, dTime, + param.dMaxEnergy, param.externalGravity.dOrbFreq, dTime, CkCallbackResumeThread()); double tDrift = CkWallTimer() - startTime; timings[activeRung].tDrift += tDrift; @@ -2165,7 +2159,7 @@ void Main::setupICs() { 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 ..."); @@ -2493,7 +2487,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, 0, 0, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, dTime, CkCallbackResumeThread()); initialForces(); @@ -2631,7 +2625,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, 0, 0, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, dTime, CkCallbackResumeThread()); if(param.bGasCooling || param.bStarForm) initCooling(); @@ -3610,12 +3604,12 @@ void Main::writeOutput(int iStep) if(verbosity) ckout << " took " << (CkWallTimer() - startTime) << " seconds." << endl; - // 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, 0, 0, - CkCallbackResumeThread()); - domainDecomp(0); 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, 0, 0, + CkCallbackResumeThread()); + domainDecomp(0); buildTree(0); if(verbosity) @@ -3677,6 +3671,15 @@ void Main::writeOutput(int iStep) << endl; } } + if (param.nSteps != 0) { // Get particles back to home + // 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, param.externalGravity.dOrbFreq, dTime, + CkCallbackResumeThread()); + domainDecomp(0); + } if(param.iBinaryOut == 6) writeNCXML(achFile); } diff --git a/ParallelGravity.h b/ParallelGravity.h index e37e09be..9d098f2c 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -1571,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); diff --git a/Sph.cpp b/Sph.cpp index 3b4d1eb6..55a1ef25 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -1112,10 +1112,10 @@ void DenDvDxSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, dvz = (-p->vPred().z + q->vPred().z)*vFac; #ifdef SLIDING_PATCH if (dx < 0.0 && (p->position[0] - q->position[0] > 0.0)) { - dvy += 1.5 * dOrbFreq * fPeriod[0]; + dvy -= 1.5 * dOrbFreq * fPeriod[0]; } else if (dx > 0.0 && (p->position[0] - q->position[0] < 0.0)) { - dvy -= 1.5 * dOrbFreq * fPeriod[0]; + dvy += 1.5 * dOrbFreq * fPeriod[0]; } #endif dvxdx += dvx*dx*rs1; diff --git a/TreePiece.cpp b/TreePiece.cpp index 83d2c522..b550bd8c 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]; } @@ -2019,9 +2021,11 @@ 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 dOrbFreq, ///< Orbital frequency of patch + double _dOrbFreq, ///< Orbital frequency of patch double dTime, const CkCallback& cb) { + totalTime = dTime; + dOrbFreq = _dOrbFreq; callback = cb; // called by assignKeys() deleteTree(); @@ -5782,6 +5786,7 @@ void TreePiece::pup(PUP::er& p) { // Periodic variables p | nReplicas; p | fPeriod; + p | dOrbFreq; p | bEwald; p | fEwCut; p | dEwhCut; From 6eb73d101add51f014677941203bd1aedfdab904 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Mon, 3 May 2021 15:16:28 -0700 Subject: [PATCH 10/20] added dorbfreq to remaining drift calls --- ParallelGravity.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index dee11611..2e08ca61 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -2487,7 +2487,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, 0, dTime, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, param.externalGravity.dOrbFreq, dTime, CkCallbackResumeThread()); initialForces(); @@ -2625,7 +2625,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, 0, dTime, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, param.externalGravity.dOrbFreq, dTime, CkCallbackResumeThread()); if(param.bGasCooling || param.bStarForm) initCooling(); @@ -2840,7 +2840,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, 0, 0, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, false, param.dMaxEnergy, param.externalGravity.dOrbFreq, dTime, CkCallbackResumeThread()); treeProxy[0].flushStarLog(CkCallbackResumeThread()); param.iStartStep = iStep; // update so that restart continues on From 050c1ff0f3b821170cf2a7f86784170cb0c07f6f Mon Sep 17 00:00:00 2001 From: brianxchen Date: Fri, 7 May 2021 10:54:49 -0700 Subject: [PATCH 11/20] attempted removal of dorbfreq in drift calls --- ParallelGravity.ci | 2 +- ParallelGravity.cpp | 16 ++++++++-------- ParallelGravity.h | 2 +- TreePiece.cpp | 2 -- 4 files changed, 10 insertions(+), 12 deletions(-) diff --git a/ParallelGravity.ci b/ParallelGravity.ci index 921ed9a6..ea1e1348 100644 --- a/ParallelGravity.ci +++ b/ParallelGravity.ci @@ -397,7 +397,7 @@ mainmodule ParallelGravity { entry void drift(double dDelta, int bNeedVPred, int bGasIsoThermal, double dvDelta, double duDelta, int nGrowMass, bool buildTree, double dMaxEnergy, - double dOrbFreq, double dTime, const CkCallback& cb); + 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 2e08ca61..0aef730c 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.externalGravity.dOrbFreq, dTime, + param.dMaxEnergy, dTime, CkCallbackResumeThread()); double tDrift = CkWallTimer() - startTime; timings[activeRung].tDrift += tDrift; @@ -2487,7 +2487,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, param.externalGravity.dOrbFreq, dTime, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime, CkCallbackResumeThread()); initialForces(); @@ -2625,7 +2625,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, param.externalGravity.dOrbFreq, dTime, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime, CkCallbackResumeThread()); if(param.bGasCooling || param.bStarForm) initCooling(); @@ -2840,7 +2840,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, param.externalGravity.dOrbFreq, dTime, + 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 +3033,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, 0, 0, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3607,7 +3607,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, 0, 0, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3652,7 +3652,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, 0, 0, + treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0, CkCallbackResumeThread()); domainDecomp(0); buildTree(0); @@ -3676,7 +3676,7 @@ void Main::writeOutput(int iStep) // 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, param.externalGravity.dOrbFreq, dTime, + 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 9d098f2c..4f7c3411 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -1696,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, - double dOrbFreq, double dTime, const CkCallback& cb); + double dTime, const CkCallback& cb); void initAccel(int iKickRung, const CkCallback& cb); #ifdef COOLING_MOLECULARH void distribLymanWerner(const CkCallback& cb); diff --git a/TreePiece.cpp b/TreePiece.cpp index b550bd8c..9be10b79 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -2021,11 +2021,9 @@ 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 _dOrbFreq, ///< Orbital frequency of patch double dTime, const CkCallback& cb) { totalTime = dTime; - dOrbFreq = _dOrbFreq; callback = cb; // called by assignKeys() deleteTree(); From ad5d55fc1612b3d482caf6930372f702e136d6d5 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Mon, 10 May 2021 15:30:52 -0700 Subject: [PATCH 12/20] fixed shear - moved checkparams before setperiodic --- ParallelGravity.cpp | 4 ++-- ParallelGravity.h | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index 0aef730c..36e82c38 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -2156,6 +2156,8 @@ 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, @@ -2289,8 +2291,6 @@ void Main::setupICs() { } } - param.externalGravity.CheckParams(prm, param); - string achLogFileName = string(param.achOutName) + ".log"; ofstream ofsLog; if(bIsRestarting) diff --git a/ParallelGravity.h b/ParallelGravity.h index 4f7c3411..95c911b7 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -2012,7 +2012,6 @@ class TreePiece : public CBase_TreePiece { int z = ((offsetcode >> 6) & 0x7) - 3; Vector3D offset(x*fPeriod.x, y*fPeriod.y + SHEAR(x, totalTime, fPeriod[1], dOrbFreq), z*fPeriod.z); - return offset; } From 91eb7bc2cb85fa49980c3b4304a3c7e141f88e97 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Sun, 13 Jun 2021 15:33:29 -0700 Subject: [PATCH 13/20] comments, whitespace changes --- GravityParticle.h | 2 +- Sph.cpp | 1 - TreePiece.cpp | 16 ++++++++++++++++ 3 files changed, 17 insertions(+), 2 deletions(-) diff --git a/GravityParticle.h b/GravityParticle.h index ec9cf19e..f959dcda 100644 --- a/GravityParticle.h +++ b/GravityParticle.h @@ -362,7 +362,7 @@ class GravityParticle : public ExternalGravityParticle { p | key; p | velocity; #ifdef SLIDING_PATCH - p | dPy; + p | dPy; ///< Canonical momentum used to update y-velocity #endif p | treeAcceleration; p | dtGrav; diff --git a/Sph.cpp b/Sph.cpp index 565c335c..0625aa83 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -1136,7 +1136,6 @@ void DenDvDxSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, dvzdx += dvz*dx*rs1; dvzdy += dvz*dy*rs1; dvzdz += dvz*dz*rs1; - divvnorm += (dx*dx+dy*dy+dz*dz)*rs1; /* Grad P estimate */ /* This used to be: diff --git a/TreePiece.cpp b/TreePiece.cpp index 2d831591..f8c427f2 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1467,6 +1467,14 @@ void TreePiece::calcEnergy(const CkCallback& cb) { #include "physconst.h" +/** +* @brief closePatch performs a velocity update to particles +* at the edge of a sliding patch. +* @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 if (bClosing) { @@ -1475,6 +1483,14 @@ inline void closePatch(int bClosing, double dDelta, GravityParticle *p, double d } #endif } +/** +* @brief openPatch performs a velocity update to particles +* at the edge of a sliding patch. +* @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 if (!bClosing) { From 1282b04861102f1166b3bf6e5280650889d63ee3 Mon Sep 17 00:00:00 2001 From: Tom Quinn Date: Mon, 14 Jun 2021 13:45:10 -0700 Subject: [PATCH 14/20] Add --enable-slidingpatch as a configure option. --- Makefile.in | 1 + configure | 20 +++++++++++++++++++- configure.ac | 5 ++++- 3 files changed, 24 insertions(+), 2 deletions(-) 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/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 From e49553e5d3f593d1ae27733c354efc1aeed30ad5 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Thu, 17 Jun 2021 17:15:51 -0700 Subject: [PATCH 15/20] documentation, removed dOrbFreq from TreePiece::kick --- GravityParticle.h | 2 +- ParallelGravity.ci | 2 +- ParallelGravity.cpp | 2 +- ParallelGravity.h | 6 +++--- TreePiece.cpp | 8 ++++---- externalGravity.cpp | 4 +++- 6 files changed, 13 insertions(+), 11 deletions(-) diff --git a/GravityParticle.h b/GravityParticle.h index f959dcda..ec9cf19e 100644 --- a/GravityParticle.h +++ b/GravityParticle.h @@ -362,7 +362,7 @@ class GravityParticle : public ExternalGravityParticle { p | key; p | velocity; #ifdef SLIDING_PATCH - p | dPy; ///< Canonical momentum used to update y-velocity + p | dPy; #endif p | treeAcceleration; p | dtGrav; diff --git a/ParallelGravity.ci b/ParallelGravity.ci index ea1e1348..6fa81501 100644 --- a/ParallelGravity.ci +++ b/ParallelGravity.ci @@ -361,7 +361,7 @@ mainmodule ParallelGravity { entry void evaluateParticleCounts(ORBSplittersMsg *splittersMsg); entry void kick(int iKickRung, double dDelta[MAXRUNG+1], int bClosing, - int bNeedVPred, int bGasIsothermal, double dOrbFreq, double dMaxEnergy, double duDelta[MAXRUNG+1], + int bNeedVPred, int bGasIsothermal, double dMaxEnergy, double duDelta[MAXRUNG+1], double gammam1, double dThermalCondSatCoeff, double dMultiPhaseMaxTime, double dMultiPhaseMinTemp, double dEvapCoeff, const CkCallback& cb); entry void initAccel(int iKickRung, const CkCallback& cb); diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index 6c0a6c31..9b60fcd2 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -1878,7 +1878,7 @@ void Main::kick(bool bClosing, int iActiveRung, int nextMaxRung, double a = csmTime2Exp(param.csm,dTime); double startTime = CkWallTimer(); treeProxy.kick(iActiveRung, dKickFac, bClosing, param.bDoGas, - param.bGasIsothermal, param.externalGravity.dOrbFreq, param.dMaxEnergy, duKick, + param.bGasIsothermal, param.dMaxEnergy, duKick, (param.dConstGamma-1), param.dThermalCondSatCoeff/a, param.feedback->dMultiPhaseMaxTime, param.feedback->dMultiPhaseMinTemp, diff --git a/ParallelGravity.h b/ParallelGravity.h index 95c911b7..2b05c29f 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -1192,7 +1192,7 @@ class TreePiece : public CBase_TreePiece { /// Background density of the Universe double dRhoFac; Vector3D fPeriod; - double dOrbFreq; + double dOrbFreq; ///< Orbital frequency of patch int nReplicas; int bEwald; /* Perform Ewald */ double fEwCut; @@ -1691,7 +1691,7 @@ class TreePiece : public CBase_TreePiece { /*****************************/ void kick(int iKickRung, double dDelta[MAXRUNG+1], int bClosing, - int bNeedVPred, int bGasIsothermal, double dOrbFreq, double dMaxEnergy, double duDelta[MAXRUNG+1], + int bNeedVPred, int bGasIsothermal, double dMaxEnergy, double duDelta[MAXRUNG+1], double gammam1, double dThermalCondSatCoeff, double dMultiPhaseMaxTime, double dMultiPhaseMinTemp, double dEvapCoeff, const CkCallback& cb); void drift(double dDelta, int bNeedVPred, int bGasIsothermal, double dvDelta, @@ -1997,7 +1997,7 @@ class TreePiece : public CBase_TreePiece { GenericTreeNode *getRoot() {return root;} // need this in Compute inline double SHEAR(int ix, ///< Interior or exterior? - double t, ///< time, in dTime+dDelta + double t, ///< Current simulation time Vector3D fPeriod, ///< vector for dxPeriod and dyPeriod double dOrbFreq) ///< Orbital frequency { diff --git a/TreePiece.cpp b/TreePiece.cpp index f8c427f2..730f5715 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1468,7 +1468,8 @@ void TreePiece::calcEnergy(const CkCallback& cb) { #include "physconst.h" /** -* @brief closePatch performs a velocity update to particles +* @brief closePatch performs a velocity update as part of the kick-cross-drift-cross-kick +algorithm described in Quinn et al. 2010 * at the edge of a sliding patch. * @param bClosing bool indicating whether the patch is opening or closing * @param dDelta timestep @@ -1484,8 +1485,8 @@ inline void closePatch(int bClosing, double dDelta, GravityParticle *p, double d #endif } /** -* @brief openPatch performs a velocity update to particles -* at the edge of a sliding patch. +* @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 @@ -1507,7 +1508,6 @@ 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 int bGasIsothermal, // Isothermal EOS - double dOrbFreq, ///< Orbital frequency double dMaxEnergy, // Maximum internal energy of gas. double duDelta[MAXRUNG+1], // dts for energy double gammam1, // Adiabatic index - 1 diff --git a/externalGravity.cpp b/externalGravity.cpp index 562f72bd..897e49fc 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -55,7 +55,9 @@ void ExternalGravity::CheckParams(PRM prm, struct parameters ¶m) // Enable external gravity if any of the flags are set if (bBodyForce || bPatch || bCentralBody) { param.bDoExternalGravity = 1; - param.externalGravity.dOrbFreq = sqrt(param.externalGravity.dCentMass / pow(param.externalGravity.dOrbDist, 3)); + if (bPatch) { + param.externalGravity.dOrbFreq = sqrt(param.externalGravity.dCentMass / pow(param.externalGravity.dOrbDist, 3)); + } } } From 716e4fc14b1db9657b8780ead16a3528589376b9 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Mon, 9 Aug 2021 14:34:36 -0700 Subject: [PATCH 16/20] changed vector components to .x, .y etc, added sliding patch ifdef to shear --- ParallelGravity.h | 6 +++++- Sph.cpp | 16 ++++++++-------- TreePiece.cpp | 9 ++++----- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/ParallelGravity.h b/ParallelGravity.h index 2b05c29f..690837af 100644 --- a/ParallelGravity.h +++ b/ParallelGravity.h @@ -2001,8 +2001,12 @@ class TreePiece : public CBase_TreePiece { 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) { @@ -2011,7 +2015,7 @@ class TreePiece : public CBase_TreePiece { int y = ((offsetcode >> 3) & 0x7) - 3; int z = ((offsetcode >> 6) & 0x7) - 3; - Vector3D offset(x*fPeriod.x, y*fPeriod.y + SHEAR(x, totalTime, fPeriod[1], dOrbFreq), 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 0625aa83..359bb0a2 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -1120,11 +1120,11 @@ void DenDvDxSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, 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[0] - q->position[0] > 0.0)) { - dvy -= 1.5 * dOrbFreq * fPeriod[0]; + if (dx < 0.0 && (p->position.x - q->position.x > 0.0)) { + dvy -= 1.5 * dOrbFreq * fPeriod.y; } - else if (dx > 0.0 && (p->position[0] - q->position[0] < 0.0)) { - dvy += 1.5 * dOrbFreq * fPeriod[0]; + else if (dx > 0.0 && (p->position.x - q->position.x < 0.0)) { + dvy += 1.5 * dOrbFreq * fPeriod.y; } #endif dvxdx += dvx*dx*rs1; @@ -1584,11 +1584,11 @@ void PressureSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, params.dx = nnList[i].dx; dv = p->vPred() - q->vPred(); #ifdef SLIDING_PATCH - if (params.dx.x < 0.0 && (p->position[0] - q->position[0] > 0.0)) { - dv[1] += 1.5 * dOrbFreq * fPeriod[0]; + 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[0] - q->position[0] < 0.0)) { - dv[1] -= 1.5 * dOrbFreq * fPeriod[0]; + 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; diff --git a/TreePiece.cpp b/TreePiece.cpp index 730f5715..69ae6f63 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1469,8 +1469,7 @@ void TreePiece::calcEnergy(const CkCallback& cb) { /** * @brief closePatch performs a velocity update as part of the kick-cross-drift-cross-kick -algorithm described in Quinn et al. 2010 -* at the edge of a sliding patch. +* 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 @@ -1479,14 +1478,14 @@ algorithm described in Quinn et al. 2010 inline void closePatch(int bClosing, double dDelta, GravityParticle *p, double dOrbFreq) { #ifdef SLIDING_PATCH if (bClosing) { - p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; - p->velocity[1] = p->dPy - 2 * dOrbFreq * p->position[0]; + p->velocity.x += 2.0 * dDelta * dOrbFreq * p->dPy; + p->velocity.y = p->dPy - 2 * dOrbFreq * p->position.x; } #endif } /** * @brief openPatch performs a velocity update as part of the kick-cross-drift-cross-kick -algorithm described in Quinn et al. 2010 +* 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 From 56b2b3295cfb461360b99a97bd34b5e0a42b1c67 Mon Sep 17 00:00:00 2001 From: brianxchen Date: Mon, 9 Aug 2021 14:46:11 -0700 Subject: [PATCH 17/20] typo fixes --- Sph.cpp | 4 ++-- TreePiece.cpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Sph.cpp b/Sph.cpp index 359bb0a2..1d3dd586 100644 --- a/Sph.cpp +++ b/Sph.cpp @@ -1121,10 +1121,10 @@ void DenDvDxSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, 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.y; + dvy -= 1.5 * dOrbFreq * fPeriod.x; } else if (dx > 0.0 && (p->position.x - q->position.x < 0.0)) { - dvy += 1.5 * dOrbFreq * fPeriod.y; + dvy += 1.5 * dOrbFreq * fPeriod.x; } #endif dvxdx += dvx*dx*rs1; diff --git a/TreePiece.cpp b/TreePiece.cpp index 69ae6f63..92cdf0bf 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1494,11 +1494,11 @@ inline void closePatch(int bClosing, double dDelta, GravityParticle *p, double d inline void openPatch(int bClosing, double dDelta, GravityParticle* p, double dOrbFreq) { #ifdef SLIDING_PATCH if (!bClosing) { - p->dPy = p->velocity[1] + 2.0 * dOrbFreq * p->position[0]; + p->dPy = p->velocity.x + 2.0 * dOrbFreq * p->position.x; // Cross hamiltonian - p->velocity[0] += 2.0 * dDelta * dOrbFreq * p->dPy; - p->velocity[1] = p->dPy - dOrbFreq * p->position[0] - dOrbFreq * (p->position[0] + 2.0 * dDelta * p->velocity[0]); + 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 } From 60dfdd9763c9c5d53d34faaf0e4174c50686f528 Mon Sep 17 00:00:00 2001 From: Tom Quinn Date: Fri, 30 Jun 2023 17:55:13 -0700 Subject: [PATCH 18/20] ExternalGravity::applyPotential(): correct Hill acceleration. --- externalGravity.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/externalGravity.cpp b/externalGravity.cpp index 897e49fc..06fd4d41 100644 --- a/externalGravity.cpp +++ b/externalGravity.cpp @@ -130,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) { From a8b18906981eac50d2cbbc2b9d6ef35a6cd482b5 Mon Sep 17 00:00:00 2001 From: Tom Quinn Date: Fri, 30 Jun 2023 17:56:19 -0700 Subject: [PATCH 19/20] Log SLIDING_PATCH and NO_HILL compile options. --- ParallelGravity.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/ParallelGravity.cpp b/ParallelGravity.cpp index 9b60fcd2..b33411ff 100644 --- a/ParallelGravity.cpp +++ b/ParallelGravity.cpp @@ -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 " From 136628748c04a1d3c386e9e6e5394bff501f76d7 Mon Sep 17 00:00:00 2001 From: Tom Quinn Date: Fri, 30 Jun 2023 17:57:28 -0700 Subject: [PATCH 20/20] SLIDING_PATCH: add predicted canonical y momentum for gas. --- GravityParticle.h | 12 ++++++++++++ TreePiece.cpp | 40 +++++++++++++++++++++++++++++++++------- 2 files changed, 45 insertions(+), 7 deletions(-) diff --git a/GravityParticle.h b/GravityParticle.h index ec9cf19e..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; @@ -413,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/TreePiece.cpp b/TreePiece.cpp index 92cdf0bf..a9dc7934 100644 --- a/TreePiece.cpp +++ b/TreePiece.cpp @@ -1477,11 +1477,13 @@ void TreePiece::calcEnergy(const CkCallback& cb) { */ 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 @@ -1493,6 +1495,7 @@ inline void closePatch(int bClosing, double dDelta, GravityParticle *p, double d */ 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; @@ -1501,6 +1504,7 @@ inline void openPatch(int bClosing, double dDelta, GravityParticle* p, double dO 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], @@ -1717,7 +1721,12 @@ 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); @@ -2065,7 +2074,9 @@ 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; + 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; @@ -2077,6 +2088,7 @@ void TreePiece::drift(double dDelta, // time step in x containing if (j == 0) { /* radial wrap */ fShear = 1.5 * dOrbFreq * fPeriod[0]; p->position[1] += SHEAR(-1, dTime + dDelta, fPeriod, dOrbFreq); + fDx = -fPeriod.x; } #endif } @@ -2085,8 +2097,9 @@ void TreePiece::drift(double dDelta, // time step in x containing 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); + fShear = -1.5 * dOrbFreq * fPeriod.x; + p->position.y += SHEAR(1, dTime + dDelta, fPeriod, dOrbFreq); + fDx = fPeriod.x; } #endif } @@ -2105,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 @@ -2157,11 +2177,17 @@ void TreePiece::drift(double dDelta, // time step in x containing #endif /* DIFFUSION */ } #ifdef SLIDING_PATCH - p->velocity[1] += fShear; - p->dPy -= fShear / 3.0; /* Angular momentum is - also changed. */ + 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");