Shear

GravityParticle.h

@@ -61,6 +61,9 @@ class extraSPHData
     double _fTimeCoolIsOffUntil;/* time cooling is turned back on */
     Vector3D<double> _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<double> _curlv;	/* Curl of the velocity */
@@ -122,6 +125,9 @@ class extraSPHData
     inline double& fESNrate() {return _fESNrate;}
     inline double& fTimeCoolIsOffUntil() {return _fTimeCoolIsOffUntil;}
     inline Vector3D<double>& vPred() {return _vPred;}
+#ifdef SLIDING_PATCH
+    inline double& dPyPred() {return _dPyPred;}
+#endif
     inline double& uPred() {return _uPred;}
     inline double& divv() {return _divv;}
     inline Vector3D<double>& curlv() {return _curlv;}
@@ -182,6 +188,9 @@ class extraSPHData
 	p | _fESNrate;
 	p | _fTimeCoolIsOffUntil;
 	p | _vPred;
+#ifdef SLIDING_PATCH
+	p | _dPyPred;
+#endif
 	p | _uPred;
 	p | _divv;
 	p | _curlv;
@@ -311,6 +320,9 @@ class GravityParticle : public ExternalGravityParticle {
 public:
 	SFC::Key key;
 	Vector3D<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 +370,9 @@ class GravityParticle : public ExternalGravityParticle {
           ExternalGravityParticle::pup(p);
           p | key;
           p | velocity;
+#ifdef SLIDING_PATCH
+          p | dPy;
+#endif
           p | treeAcceleration;
           p | dtGrav;
           p | fDensity;
@@ -407,6 +422,9 @@ class GravityParticle : public ExternalGravityParticle {
 	inline double& fESNrate() {IMAGAS; return (((extraSPHData*)extraData)->fESNrate());}
 	inline double& fTimeCoolIsOffUntil() {IMAGAS; return (((extraSPHData*)extraData)->fTimeCoolIsOffUntil());}
 	inline Vector3D<double>& 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<double>& curlv() { IMAGAS; return (((extraSPHData*)extraData)->curlv());}

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);
@@ -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

@@ -1971,7 +1971,7 @@ void Main::advanceBigStep(int iStep) {
 	      bool bBuildTree = (iSub + 1 == driftSteps);
 	      treeProxy.drift(dDriftFac, param.bDoGas, param.bGasIsothermal,
 			      dKickFac, dTimeSub, nGrowMassDrift, bBuildTree,
-                              param.dMaxEnergy,
+                              param.dMaxEnergy, dTime,
 			      CkCallbackResumeThread());
               double tDrift = CkWallTimer() - startTime;
               timings[activeRung].tDrift += tDrift;
@@ -2156,10 +2156,12 @@ void Main::advanceBigStep(int iStep) {
 void Main::setupICs() {
   double startTime;
 
+  param.externalGravity.CheckParams(prm, param);
+
   treeProxy.setPeriodic(param.nReplicas, param.vPeriod, param.bEwald,
 			param.dEwCut, param.dEwhCut, param.bPeriodic,
                         param.csm->bComove,
-                        0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0);
+                        0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0, param.externalGravity.dOrbFreq);
 
   /******** Particles Loading ********/
   CkPrintf("Loading particles ...");
@@ -2289,8 +2291,6 @@ void Main::setupICs() {
           }
       } 
 
-  param.externalGravity.CheckParams(prm, param);
-
   string achLogFileName = string(param.achOutName) + ".log";
   ofstream ofsLog;
   if(bIsRestarting)
@@ -2423,6 +2423,12 @@ void Main::setupICs() {
 #endif
 #ifdef DAMPING
   ofsLog << " DAMPING";
+#endif
+#ifdef SLIDING_PATCH
+  ofsLog << " SLIDING_PATCH";
+#endif
+#ifdef NO_HILL
+  ofsLog << " NO_HILL";
 #endif
   ofsLog << endl;
   ofsLog << "# Key sizes: " << sizeof(KeyType) << " bytes particle "
@@ -2487,7 +2493,7 @@ void Main::setupICs() {
 
 // for periodic, puts all particles within the boundary
 // Also assigns keys and sorts.
-  treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+  treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime,
                   CkCallbackResumeThread());
 
   initialForces();
@@ -2625,7 +2631,7 @@ Main::restart(CkCheckpointStatusMsg *msg)
         } else {
             CkPrintf("Not Using CkLoop %d\n", param.bUseCkLoopPar);
         }
-        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime,
                         CkCallbackResumeThread());
 	if(param.bGasCooling || param.bStarForm) 
 	    initCooling();
@@ -2840,7 +2846,7 @@ Main::doSimulation()
 	    }
 	// The following drift is called because it deletes the tree
 	// so it won't be saved on disk.
-	treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, false, param.dMaxEnergy,
+	treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, false, param.dMaxEnergy, dTime,
                         CkCallbackResumeThread());
 	treeProxy[0].flushStarLog(CkCallbackResumeThread());
 	param.iStartStep = iStep; // update so that restart continues on
@@ -3033,7 +3039,7 @@ Main::doSimulation()
       if(param.bDoGas && param.bDoDensity) {
 	  // The following call is to get the particles in key order
 	  // before the sort.
-	  treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+	  treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0,
                           CkCallbackResumeThread());
           domainDecomp(0);
           buildTree(0);
@@ -3607,7 +3613,7 @@ void Main::writeOutput(int iStep)
     if(param.nSteps != 0 && param.bDoDensity) {
         // The following call is to get the particles in key order
         // before the sort.
-        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0,
                         CkCallbackResumeThread());
         domainDecomp(0);
         buildTree(0);
@@ -3652,7 +3658,7 @@ void Main::writeOutput(int iStep)
 	    startTime = CkWallTimer();
 	    // The following call is to get the particles in key order
 	    // before the sort.
-	    treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+	    treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0,
                             CkCallbackResumeThread());
             domainDecomp(0);
             buildTree(0);
@@ -3675,7 +3681,7 @@ void Main::writeOutput(int iStep)
                                 // processors for continuing the simulation.
         // The following call is to get the particles in key order
         // before the sort.
-        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime,
                         CkCallbackResumeThread());
         domainDecomp(0);
     }

ParallelGravity.h

@@ -1192,6 +1192,7 @@ class TreePiece : public CBase_TreePiece {
 	/// Background density of the Universe
 	double dRhoFac;
 	Vector3D<cosmoType> fPeriod;
+	double dOrbFreq; ///< Orbital frequency of patch
 	int nReplicas;
 	int bEwald;		/* Perform Ewald */
 	double fEwCut;
@@ -1570,7 +1571,7 @@ class TreePiece : public CBase_TreePiece {
 
 	void setPeriodic(int nReplicas, Vector3D<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);
@@ -1695,7 +1696,7 @@ class TreePiece : public CBase_TreePiece {
         double dMultiPhaseMaxTime, double dMultiPhaseMinTemp, double dEvapCoeff, const CkCallback& cb);
   void drift(double dDelta, int bNeedVPred, int bGasIsothermal, double dvDelta,
              double duDelta, int nGrowMass, bool buildTree, double dMaxEnergy,
-	     const CkCallback& cb);
+	         double dTime, const CkCallback& cb);
   void initAccel(int iKickRung, const CkCallback& cb);
 #ifdef COOLING_MOLECULARH
   void distribLymanWerner(const CkCallback& cb);
@@ -1995,14 +1996,26 @@ class TreePiece : public CBase_TreePiece {
         // need this in TreeWalk
         GenericTreeNode *getRoot() {return root;}
         // need this in Compute
+	inline double SHEAR(int ix,                      ///< Interior or exterior?
+			double t,                    ///< Current simulation time
+			Vector3D<cosmoType> 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<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.y, 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.x - q->position.x > 0.0)) {
+            dvy -= 1.5 * dOrbFreq * fPeriod.x;
+        }
+        else if (dx > 0.0 && (p->position.x - q->position.x < 0.0)) {
+            dvy += 1.5 * dOrbFreq * fPeriod.x;
+        }
+#endif
 		dvxdx += dvx*dx*rs1;
 		dvxdy += dvx*dy*rs1;
 		dvxdz += dvx*dz*rs1;
@@ -1575,6 +1583,14 @@ void PressureSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth,
         qParams.rNorm = rs1 * q->mass;
         params.dx = nnList[i].dx;
         dv = p->vPred() - q->vPred();
+#ifdef SLIDING_PATCH
+        if (params.dx.x < 0.0 && (p->position.x - q->position.x > 0.0)) {
+            dv[1] += 1.5 * dOrbFreq * fPeriod.x;
+        }
+        else if (params.dx.x > 0.0 && (p->position.x - q->position.x < 0.0)) {
+            dv[1] -= 1.5 * dOrbFreq * fPeriod.x;
+        }
+#endif
         params.dvdotdr = vFac*dot(dv, params.dx) + fDist2*H;
 #ifdef RTFORCE
         pParams.PoverRho2 = p->PoverRho2()*p->fDensity/q->fDensity;