Shear

GravityParticle.h

@@ -311,6 +311,9 @@ class GravityParticle : public ExternalGravityParticle {
 public:
 	SFC::Key key;
 	Vector3D<double> velocity;
+#ifdef SLIDING_PATCH
+        double dPy;     ///< Canonical momentum used to update y-velocity
+#endif
 	Vector3D<cosmoType> treeAcceleration;
 	cosmoType potential;
         cosmoType dtGrav;       ///< timestep from gravity; N.B., this
@@ -358,6 +361,9 @@ class GravityParticle : public ExternalGravityParticle {
           ExternalGravityParticle::pup(p);
           p | key;
           p | velocity;
+#ifdef SLIDING_PATCH
+          p | dPy; ///< Canonical momentum used to update y-velocity
+#endif
           p | treeAcceleration;
           p | dtGrav;
           p | fDensity;

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@   \

ParallelGravity.ci

@@ -271,7 +271,7 @@ mainmodule ParallelGravity {
     entry void resetObjectLoad(const CkCallback& cb);
     entry void setPeriodic(int nReplicas, Vector3D<cosmoType> 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);
@@ -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 dMaxEnergy, double duDelta[MAXRUNG+1],
+	    int bNeedVPred, int bGasIsothermal, double dOrbFreq, 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);
@@ -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);

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.dMaxEnergy, duKick,
+                   param.bGasIsothermal, param.externalGravity.dOrbFreq, param.dMaxEnergy, duKick,
                    (param.dConstGamma-1), param.dThermalCondSatCoeff/a,
                    param.feedback->dMultiPhaseMaxTime,
                    param.feedback->dMultiPhaseMinTemp,
@@ -1971,7 +1971,7 @@ void Main::advanceBigStep(int iStep) {
 	      bool bBuildTree = (iSub + 1 == driftSteps);
 	      treeProxy.drift(dDriftFac, param.bDoGas, param.bGasIsothermal,
 			      dKickFac, dTimeSub, nGrowMassDrift, bBuildTree,
-                              param.dMaxEnergy,
+                              param.dMaxEnergy, dTime,
 			      CkCallbackResumeThread());
               double tDrift = CkWallTimer() - startTime;
               timings[activeRung].tDrift += tDrift;
@@ -2156,10 +2156,12 @@ void Main::advanceBigStep(int iStep) {
 void Main::setupICs() {
   double startTime;
 
+  param.externalGravity.CheckParams(prm, param);
+
   treeProxy.setPeriodic(param.nReplicas, param.vPeriod, param.bEwald,
 			param.dEwCut, param.dEwhCut, param.bPeriodic,
                         param.csm->bComove,
-                        0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0);
+                        0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0, param.externalGravity.dOrbFreq);
 
   /******** Particles Loading ********/
   CkPrintf("Loading particles ...");
@@ -2289,8 +2291,6 @@ void Main::setupICs() {
           }
       } 
 
-  param.externalGravity.CheckParams(prm, param);
-
   string achLogFileName = string(param.achOutName) + ".log";
   ofstream ofsLog;
   if(bIsRestarting)
@@ -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,
+  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,
+        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,
+	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,
+	  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,
+        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,
+	    treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0,
                             CkCallbackResumeThread());
             domainDecomp(0);
             buildTree(0);
@@ -3675,7 +3675,7 @@ void Main::writeOutput(int iStep)
                                 // processors for continuing the simulation.
         // The following call is to get the particles in key order
         // before the sort.
-        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime,
                         CkCallbackResumeThread());
         domainDecomp(0);
     }

ParallelGravity.h

@@ -1192,6 +1192,7 @@ class TreePiece : public CBase_TreePiece {
 	/// Background density of the Universe
 	double dRhoFac;
 	Vector3D<cosmoType> fPeriod;
+	double dOrbFreq;
 	int nReplicas;
 	int bEwald;		/* Perform Ewald */
 	double fEwCut;
@@ -1570,7 +1571,7 @@ class TreePiece : public CBase_TreePiece {
 
 	void setPeriodic(int nReplicas, Vector3D<cosmoType> 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);
@@ -1690,12 +1691,12 @@ class TreePiece : public CBase_TreePiece {
   /*****************************/
 
   void kick(int iKickRung, double dDelta[MAXRUNG+1], int bClosing,
-	    int bNeedVPred, int bGasIsothermal, double dMaxEnergy, double duDelta[MAXRUNG+1],
+	    int bNeedVPred, int bGasIsothermal, double dOrbFreq, 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,
              double duDelta, int nGrowMass, bool buildTree, double dMaxEnergy,
-	     const CkCallback& cb);
+	         double dTime, const CkCallback& cb);
   void initAccel(int iKickRung, const CkCallback& cb);
 #ifdef COOLING_MOLECULARH
   void distribLymanWerner(const CkCallback& cb);
@@ -1995,14 +1996,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<cosmoType> 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<cosmoType> 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<cosmoType> offset(x*fPeriod.x, y*fPeriod.y, z*fPeriod.z);
-
+	    Vector3D<cosmoType> offset(x*fPeriod.x, y*fPeriod.y + SHEAR(x, totalTime, fPeriod[1], dOrbFreq), z*fPeriod.z);
 	    return offset;
 	    }
 

Sph.cpp

@@ -30,7 +30,7 @@ Main::initSph()
         // Starting is true
 	DenDvDxSmoothParams pDen(TYPE_GAS, 0, param.csm, dTime, 0,
 				 param.bConstantDiffusion, 1, bHaveAlpha,
-                                 param.dConstAlphaMax);
+                                 param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod);
 	double startTime = CkWallTimer();
 	double dfBall2OverSoft2 = 4.0*param.dhMinOverSoft*param.dhMinOverSoft;
 	treeProxy.startSmooth(&pDen, 1, param.nSmooth, dfBall2OverSoft2,
@@ -712,7 +712,7 @@ Main::doSph(int activeRung, int bNeedDensity)
 	// This also marks neighbors of actives
 	DenDvDxSmoothParams pDen(TYPE_GAS, activeRung, param.csm, dTime, 1,
 				 param.bConstantDiffusion, 0, 0,
-                                 param.dConstAlphaMax);
+                                 param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod);
 	double startTime = CkWallTimer();
 	treeProxy.startSmooth(&pDen, 1, param.nSmooth, dfBall2OverSoft2,
 			      CkCallbackResumeThread());
@@ -732,7 +732,7 @@ Main::doSph(int activeRung, int bNeedDensity)
 	// additional marking
 	DenDvDxNeighborSmParams pDenN(TYPE_GAS, activeRung, param.csm, dTime,
 				      param.bConstantDiffusion,
-                                      param.dConstAlphaMax);
+                                      param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod);
 	startTime = CkWallTimer();
 	treeProxy.startSmooth(&pDenN, 1, param.nSmooth, dfBall2OverSoft2,
 			      CkCallbackResumeThread());
@@ -745,7 +745,7 @@ Main::doSph(int activeRung, int bNeedDensity)
 	// actives, and those who have actives as neighbors.
 	DenDvDxSmoothParams pDen(TYPE_GAS, activeRung, param.csm, dTime, 0,
 				 param.bConstantDiffusion, 0, 0,
-                                 param.dConstAlphaMax);
+                                 param.dConstAlphaMax, param.externalGravity.dOrbFreq, param.fPeriod);
 	double startTime = CkWallTimer();
 	treeProxy.startSmooth(&pDen, 1, param.nSmooth, dfBall2OverSoft2,
 			      CkCallbackResumeThread());
@@ -779,7 +779,7 @@ Main::doSph(int activeRung, int bNeedDensity)
     PressureSmoothParams pPressure(TYPE_GAS, activeRung, param.csm, dTime,
                                    param.dConstAlpha, param.dConstBeta,
                                    param.dThermalDiffusionCoeff, param.dMetalDiffusionCoeff,
-                                   param.dEtaCourant, param.dEtaDiffusion);
+                                   param.dEtaCourant, param.dEtaDiffusion, param.externalGravity.dOrbFreq, param.fPeriod);
     double startTime = CkWallTimer();
     treeProxy.startReSmooth(&pPressure, CkCallbackResumeThread());
     ckout << " took " << (CkWallTimer() - startTime) << " seconds."
@@ -1119,6 +1119,14 @@ void DenDvDxSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth,
 		dvx = (-p->vPred().x + q->vPred().x)*vFac;
 		dvy = (-p->vPred().y + q->vPred().y)*vFac;
 		dvz = (-p->vPred().z + q->vPred().z)*vFac;
+#ifdef SLIDING_PATCH
+        if (dx < 0.0 && (p->position[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;
@@ -1575,6 +1583,14 @@ void PressureSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth,
         qParams.rNorm = rs1 * q->mass;
         params.dx = nnList[i].dx;
         dv = p->vPred() - q->vPred();
+#ifdef SLIDING_PATCH
+        if (params.dx.x < 0.0 && (p->position[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;

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<cosmoType> 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<cosmoType> _fPeriod) {
 	iType = _iType;
 	activeRung = am;
 	bActiveOnly = _bActiveOnly;
@@ -55,6 +58,8 @@ class DenDvDxSmoothParams : public SmoothParams
         bHaveAlpha = _bHaveAlpha;
 	dAlphaMax = _dAlphaMax;
         dTime = _dTime;
+        dOrbFreq = _dOrbFreq;
+        fPeriod = _fPeriod;
 	if(csm->bComove) {
 	    H = csmTime2Hub(csm,dTime);
 	    a = csmTime2Exp(csm,dTime);
@@ -76,6 +81,8 @@ class DenDvDxSmoothParams : public SmoothParams
 	p|bStarting;
 	p|bHaveAlpha;
 	p|dAlphaMax;
+    p|dOrbFreq;
+    p|fPeriod;
 	}
     };
 
@@ -107,9 +114,9 @@ class DenDvDxNeighborSmParams : public DenDvDxSmoothParams
     /// This calls the DenDvDx constructor, and we assume bActiveOnly,
     /// bStarting, and bHaveAlpha are not set.
     DenDvDxNeighborSmParams(int _iType, int am, CSM csm, double dTime,
-			    int bConstDiffusion, double dAlphaMax)
+			    int bConstDiffusion, double dAlphaMax, double dOrbFreq, Vector3D<cosmoType> _fPeriod)
       : DenDvDxSmoothParams(_iType, am, csm, dTime, 0, bConstDiffusion,
-                            0, 0, dAlphaMax) {}
+                            0, 0, dAlphaMax, dOrbFreq, fPeriod) {}
     PUPable_decl(DenDvDxNeighborSmParams);
     DenDvDxNeighborSmParams(CkMigrateMessage *m) : DenDvDxSmoothParams(m) {}
     virtual void pup(PUP::er &p) {
@@ -156,6 +163,8 @@ class PressureSmoothParams : public SmoothParams
     double dMetalDiffusionCoeff;
     double dtFacCourant; // Courant timestep factor
     double dtFacDiffusion; // Diffusion timestep factor
+    double dOrbFreq; // Orbital Frequency
+    Vector3D<cosmoType> 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<cosmoType> _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;
 	}
     };