Shearing patch

GravityParticle.h

@@ -60,7 +60,10 @@ class extraSPHData
     double _fESNrate;		/* SN energy rate  */
     double _fTimeCoolIsOffUntil;/* time cooling is turned back on */
     Vector3D<double> _vPred;	/* Predicted velocities for velocity
-				   dependent forces */
+                                   dependent forces */
+#ifdef SLIDING_PATCH
+    double _dPyPred;            ///< Predicted canonical momentum
+#endif
     double _uPred;		/* Predicted internal energy */
     double _divv;		/* Diverence of the velocity */
     Vector3D<double> _curlv;	/* Curl of the velocity */
@@ -124,6 +127,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;}
@@ -186,6 +192,9 @@ class extraSPHData
 	p | _fESNrate;
 	p | _fTimeCoolIsOffUntil;
 	p | _vPred;
+#ifdef SLIDING_PATCH
+        p | _dPyPred;
+#endif
 	p | _uPred;
 	p | _divv;
 	p | _curlv;
@@ -317,6 +326,10 @@ class GravityParticle : public ExternalGravityParticle {
 public:
 	SFC::Key key;
 	Vector3D<double> velocity;
+#ifdef SLIDING_PATCH
+        double dPy;         ///< Canonical momentum used to update y-velocity
+#endif
+    // inline Vector3D<double>& vPred() { return _vPred; }
 	Vector3D<cosmoType> treeAcceleration;
 	cosmoType potential;
         cosmoType dtGrav;       ///< timestep from gravity; N.B., this
@@ -364,6 +377,9 @@ class GravityParticle : public ExternalGravityParticle {
           ExternalGravityParticle::pup(p);
           p | key;
           p | velocity;
+#ifdef SLIDING_PATCH
+          p | dPy;
+#endif
           p | treeAcceleration;
           p | dtGrav;
           p | fDensity;
@@ -413,6 +429,9 @@ class GravityParticle : public ExternalGravityParticle {
 	inline double& fESNrate() {IMAGAS; return (((extraSPHData*)extraData)->fESNrate());}
 	inline double& fTimeCoolIsOffUntil() {IMAGAS; return (((extraSPHData*)extraData)->fTimeCoolIsOffUntil());}
 	inline Vector3D<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

@@ -126,6 +126,7 @@ defines := $(strip @HEXADECAPOLE@ @FLAG_GDFORCE@ @FLAG_ARCH@  \
           @FLAG_JEANSSOFTONLY@ \
           @FLAG_TREE_BUILD@ $(debug_defines) @FLAG_INTERLIST@ \
           @FLAG_NSMOOTHINNER@ @FLAG_SPLITGAS@ @FLAG_SIDMINTERACT@ \
+          @FLAG_SLIDING_PATCH@ \
           @FLAG_SUPERBUBBLE@ $(cuda_defines) -DREDUCTION_HELPER)
 
 modules := $(strip -language charm++ -balancer @DEFAULT_LB@   \

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);
@@ -440,9 +441,9 @@ mainmodule ParallelGravity {
     entry void liveVizDumpFrameInit(liveVizRequestMsg * msg);
     entry void setProjections(int bOn);
 #ifdef PUSH_GRAVITY
-    entry void buildTree(int bucketSize, const CkCallback& cb, bool merge);
+    entry void buildTree(int bucketSize, double dTime, const CkCallback& cb, bool merge);
 #else
-    entry void buildTree(int bucketSize, const CkCallback& cb);
+    entry void buildTree(int bucketSize, double dTime, const CkCallback& cb);
 #endif
 
     entry void startOctTreeBuild(CkReductionMsg* m);

ParallelGravity.cpp

@@ -1695,9 +1695,9 @@ void Main::buildTree(int iPhase)
     CkPrintf("Building trees ... ");
     double startTime = CkWallTimer();
 #ifdef PUSH_GRAVITY
-    treeProxy.buildTree(bucketSize, CkCallbackResumeThread(),!bDoPush);
+    treeProxy.buildTree(bucketSize, dTime, CkCallbackResumeThread(),!bDoPush);
 #else
-    treeProxy.buildTree(bucketSize, CkCallbackResumeThread());
+    treeProxy.buildTree(bucketSize, dTime, CkCallbackResumeThread());
 #endif
     double tTB =  CkWallTimer()-startTime;
     timings[iPhase].tTBuild += tTB;
@@ -1970,7 +1970,7 @@ void Main::advanceBigStep(int iStep) {
 	      bool bBuildTree = (iSub + 1 == driftSteps);
 	      treeProxy.drift(dDriftFac, param.bDoGas, param.bGasIsothermal,
 			      dKickFac, dTimeSub, nGrowMassDrift, bBuildTree,
-                              param.dMaxEnergy,
+                              param.dMaxEnergy, dTime,
 			      CkCallbackResumeThread());
               double tDrift = CkWallTimer() - startTime;
               timings[activeRung].tDrift += tDrift;
@@ -2160,10 +2160,13 @@ void Main::advanceBigStep(int iStep) {
 void Main::setupICs() {
   double startTime;
 
+  param.externalGravity.CheckParams(prm, param);
+
   treeProxy.setPeriodic(param.nReplicas, param.vPeriod, param.bEwald,
 			param.dEwCut, param.dEwhCut, param.bPeriodic,
                         param.csm->bComove,
-                        0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0);
+                        0.5*param.csm->dHubble0*param.csm->dHubble0*param.csm->dOmega0,
+                        param.externalGravity.dOrbFreq);
 
   /******** Particles Loading ********/
   CkPrintf("Loading particles ...");
@@ -2300,8 +2303,6 @@ void Main::setupICs() {
       if(param.iStartStep) restartNSIDM();
   }
 
-  param.externalGravity.CheckParams(prm, param);
-
   string achLogFileName = string(param.achOutName) + ".log";
   ofstream ofsLog;
   if(bIsRestarting)
@@ -2433,6 +2434,12 @@ void Main::setupICs() {
 #endif
 #ifdef DAMPING
   ofsLog << " DAMPING";
+#endif
+#ifdef SLIDING_PATCH
+  ofsLog << " SLIDING_PATCH";
+#endif
+#ifdef NO_HILL
+  ofsLog << " NO_HILL";
 #endif
   ofsLog << endl;
   ofsLog << "# Key sizes: " << sizeof(KeyType) << " bytes particle "
@@ -2496,7 +2503,7 @@ void Main::setupICs() {
 
 // for periodic, puts all particles within the boundary
 // Also assigns keys and sorts.
-  treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+  treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime,
                   CkCallbackResumeThread());
 
   initialForces();
@@ -2634,7 +2641,7 @@ Main::restart(CkCheckpointStatusMsg *msg)
         } else {
             CkPrintf("Not Using CkLoop %d\n", param.bUseCkLoopPar);
         }
-        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime,
                         CkCallbackResumeThread());
 	if(param.bGasCooling || param.bStarForm) 
 	    initCooling();
@@ -2849,7 +2856,7 @@ Main::doSimulation()
 	    }
 	// The following drift is called because it deletes the tree
 	// so it won't be saved on disk.
-	treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, false, param.dMaxEnergy,
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, false, param.dMaxEnergy, dTime,
                         CkCallbackResumeThread());
 	treeProxy[0].flushStarLog(CkCallbackResumeThread());
 	param.iStartStep = iStep; // update so that restart continues on
@@ -3039,7 +3046,7 @@ Main::doSimulation()
       if(param.bDoGas && param.bDoDensity) {
 	  // The following call is to get the particles in key order
 	  // before the sort.
-	  treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+          treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0,
                           CkCallbackResumeThread());
           domainDecomp(0);
           buildTree(0);
@@ -3613,8 +3620,8 @@ void Main::writeOutput(int iStep)
     if(param.nSteps != 0 && param.bDoDensity) {
         // The following call is to get the particles in key order
         // before the sort.
-        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
-                        CkCallbackResumeThread());
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0,
+            CkCallbackResumeThread());
         domainDecomp(0);
         buildTree(0);
 
@@ -3658,7 +3665,7 @@ void Main::writeOutput(int iStep)
 	    startTime = CkWallTimer();
 	    // The following call is to get the particles in key order
 	    // before the sort.
-	    treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+	    treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, 0,
                             CkCallbackResumeThread());
             domainDecomp(0);
             buildTree(0);
@@ -3681,11 +3688,10 @@ void Main::writeOutput(int iStep)
                                 // processors for continuing the simulation.
         // The following call is to get the particles in key order
         // before the sort.
-        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy,
+        treeProxy.drift(0.0, 0, 0, 0.0, 0.0, 0, true, param.dMaxEnergy, dTime,
                         CkCallbackResumeThread());
         domainDecomp(0);
     }
-
     if(param.iBinaryOut == 6)
         writeNCXML(achFile);
     }

ParallelGravity.h

@@ -1203,6 +1203,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;
@@ -1330,10 +1331,9 @@ class TreePiece : public CBase_TreePiece {
   // called when a chunk has been used completely (chunkRemaining[chunk] == 0)
   void finishedChunk(int chunk);
 
-  ///Array of comparison function pointers
+  /// Array of comparison function pointers for ORB tree build
   bool (*compFuncPtr[3])(GravityParticle,GravityParticle);
-  double tmpTime;
-  double totalTime;
+  double dTimeTree;             ///< Simulation time tree is built
  public:
 
 #ifdef SPCUDA
@@ -1468,8 +1468,7 @@ class TreePiece : public CBase_TreePiece {
 #endif
 #endif
 
-	  tmpTime=0.0;
-	  totalTime=0.0;
+	  dTimeTree = 0.0;
 	  // temporarely set to -1, it will updated after the tree is built
 	  numChunks=-1;
 	  prefetchRoots = NULL;
@@ -1581,7 +1580,7 @@ class TreePiece : public CBase_TreePiece {
 
 	void setPeriodic(int nReplicas, Vector3D<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);
@@ -1706,7 +1705,7 @@ class TreePiece : public CBase_TreePiece {
         double dMultiPhaseMaxTime, double dMultiPhaseMinTemp, double dEvapCoeff, const CkCallback& cb);
   void drift(double dDelta, int bNeedVPred, int bGasIsothermal, double dvDelta,
              double duDelta, int nGrowMass, bool buildTree, double dMaxEnergy,
-	     const CkCallback& cb);
+             double dTime, const CkCallback& cb);
   void initAccel(int iKickRung, const CkCallback& cb);
 #ifdef COOLING_MOLECULARH
   void distribLymanWerner(const CkCallback& cb);
@@ -1825,10 +1824,13 @@ class TreePiece : public CBase_TreePiece {
 	void setProjections(int bOn);
 
 	/// \brief Charm entry point to build the tree (called by Main).
+        /// \param bucketSize Maximum particles per bucket
+        /// \param dTime Simulation time tree is built
+        /// \param cb Callback after tree is built
 #ifdef PUSH_GRAVITY
-	void buildTree(int bucketSize, const CkCallback& cb, bool merge);
+        void buildTree(int bucketSize, double dTime, const CkCallback& cb, bool merge);
 #else
-	void buildTree(int bucketSize, const CkCallback& cb);
+        void buildTree(int bucketSize, double dTime, const CkCallback& cb);
 #endif
 
 	/// \brief Real tree build, independent of other TreePieces.
@@ -2005,15 +2007,32 @@ class TreePiece : public CBase_TreePiece {
         // jetley
         // need this in TreeWalk
         GenericTreeNode *getRoot() {return root;}
-        // need this in Compute
+        /// @brief adjust y (azimuthal) position if we are crossing an
+        /// x (radial) boundary
+        /// @return change in y position
+        inline double SHEAR(int ix,                  ///< Interior or exterior?
+                            double t,                ///< Current simulation time
+                            Vector3D<cosmoType> fPeriod, ///< patch size
+                            double dOrbFreq)         ///< Orbital frequency
+        {
+#ifdef SLIDING_PATCH
+            return
+                (ix < 0) ? fmod(0.5 * fPeriod.y - 1.5 * ix * dOrbFreq * fPeriod.x * t, fPeriod.y) - 0.5 * fPeriod.y :
+                (ix > 0) ? 0.5 * fPeriod.y - fmod(0.5 * fPeriod.y + 1.5 * ix * dOrbFreq * fPeriod.x * t, fPeriod.y) : 0.0;
+#else
+            return 0.0;
+#endif
+        }
+
 	inline Vector3D<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, dTimeTree, fPeriod.y, dOrbFreq),
+                                       z*fPeriod.z);
 	    return offset;
 	    }