Commits

Nathan Goldbaum committed ca7eb58

Active particle velocities are now updated correctly.

  • Participants
  • Parent commits b17e828

Comments (0)

Files changed (4)

File src/enzo/grid/Make.module.objects

 	Grid_FindMinimumParticleMass.o \
 	Grid_FindNewStarParticles.o \
 	Grid_GetActiveParticlePosition.o \
-	Grid_GetActiveParticleVelocity.o \
 	Grid_InterpolateParticlePositions.o \
 	Grid_MirrorStarParticles.o \
 	Grid_MoveAllParticles.o \

File src/enzo/grid/particles/Grid_GetActiveParticleVelocity.C

-/***********************************************************************
-/
-/  GRID CLASS (RETURN ACTIVE PARTICLE VELOCITIES AS A 2D POINTER ARRAY)
-/
-/  written by: Nathan Goldbaum
-/  date:       November, 2012
-/  modified1:  
-/
-/
-************************************************************************/
-#include "preincludes.h"
-
-#include "ErrorExceptions.h"
-#include "macros_and_parameters.h"
-#include "typedefs.h"
-#include "global_data.h"
-#include "Fluxes.h"
-#include "GridList.h"
-#include "ExternalBoundary.h"
-#include "Grid.h"
-#include "Hierarchy.h"
-#include "TopGridData.h"
-#include "fortran.def"
-#include "CosmologyParameters.h"
-
-#include "ActiveParticle.h"
-
-void grid::GetActiveParticleVelocity(float *ActiveParticleVelocity[]) 
-{
-  int i, dim;
-
-  for (i = 0; i < NumberOfActiveParticles; i++) {
-    FLOAT* vel = ActiveParticles[i]->ReturnVelocity();
-    for (dim = 0; dim < GridRank; dim++)
-      ActiveParticleVelocity[dim][i] = vel[dim];
-  }
-
-}

File src/enzo/grid/particles/Grid_UpdateParticleVelocity.C

 #endif
   int i, dim, dim1;
  
+  FLOAT coef, coef1, coef2;
+
   /* If using comoving coordinates, divide by a(t) first. */
- 
-  if (ComovingCoordinates)
+    
+  if (ComovingCoordinates) {
     if (CosmologyComputeExpansionFactor(Time + TimeStep, &a, &dadt)
 	== FAIL) {
             ENZO_FAIL("Error in CsomologyComputeExpansionFactors.");
     }
+
+    coef = 0.5*dadt/a*TimeStep;
+    coef1 = 1.0 - coef;
+    coef2 = 1.0 / (1.0 + coef);
+  }  
+
  
   /* Loop over dimensions. */
  
             ENZO_FAIL("No ParticleAccleration present.");
     }
  
-    if (NumberOfParticles > 0 && ParticleAcceleration[dim] == NULL) {
-            ENZO_FAIL("No ParticleAccleration present.");
-    }
- 
-
     /* Update velocities.  */
  
     if (ComovingCoordinates) {
       
-      FLOAT coef = 0.5*dadt/a*TimeStep;
-      FLOAT coef1 = 1.0 - coef;
-      FLOAT coef2 = 1.0 / (1.0 + coef);
-      
       /* If using comoving coordinates, subtract the (time-centered)
 	 drag-like term and add the acceleration. The acceleration has
 	 already been divided by a(t). */
 #endif /* VELOCITY_METHOD3 */
 	
       }
+    }
+    else {
+      
+      /* Otherwise, just add the acceleration. */
+   
+      for (i = 0; i < NumberOfParticles; i++)
+	ParticleVelocity[dim][i] += ParticleAcceleration[dim][i] * TimeStep;
+    }
+  }
 
-      if (NumberOfActiveParticles > 0) {
-	  
-	float **ActiveParticleVelocity = new float*[GridRank];
-	for (dim1 = 0; dim1 < GridRank; dim1++)
-	  ActiveParticleVelocity[dim1] = new float[NumberOfActiveParticles];
-	this->GetActiveParticleVelocity(ActiveParticleVelocity);
+
+  if (ProblemType == 29)
+    for (i = 0; i < NumberOfParticles; i++)
+      printf("id=%"PISYM"  %"PSYM" %"PSYM" %"PSYM"\n", ParticleNumber[i],
+	     ParticlePosition[0][i], ParticlePosition[1][i], ParticlePosition[2][i]);
+
+  
+  if (NumberOfActiveParticles > 0) {
+    if (ComovingCoordinates) {
+      if (CosmologyComputeExpansionFactor(Time + TimeStep, &a, &dadt)
+	  == FAIL) {
+	ENZO_FAIL("Error in CsomologyComputeExpansionFactors.");
+      }
 	
-	for (i = 0; i < NumberOfActiveParticles; i++) {
+      for (i = 0; i < NumberOfActiveParticles; i++) {
 	
+	float* apvel = ActiveParticles[i]->ReturnVelocity();
+
+	for (dim = 0; dim < GridRank; dim++) {
+
 #ifdef VELOCITY_METHOD1
 	
-          /* i) partially time-centered. */
+	  /* i) partially time-centered. */
 	
-	  VelocityMidStep = ActiveParticleVelocity[dim][i] +
+	  VelocityMidStep = apvel[dim] +
 	    ActiveParticleAcceleration[dim][i]*0.5*TimeStep;
 	
-	  ActiveParticleVelocity[dim][i] +=
+	  apvel[dim] += 
 	    (-VelocityMidStep*dadt/a + ActiveParticleAcceleration[dim][i]) * TimeStep;
 	
 #endif /* VELOCITY_METHOD1 */
 	  
 	  /* ii) partially backward. */
 	  
-	  VelocityMidStep = ActiveParticleVelocity[dim][i] ;
+	  VelocityMidStep = apvel[dim];
 	  
-	  ActiveParticleVelocity[dim][i] +=
+	  apvel[dim] +=
 	    (-VelocityMidStep*dadt/a + ActiveParticleAcceleration[dim][i]) * TimeStep;
 	  
 #endif /* VELOCITY_METHOD2 */
 	  
 	  /* iii) Semi-implicit way */
 	  
-	  ActiveParticleVelocity[dim][i] = (coef1*ActiveParticleVelocity[dim][i] +
-	                         ActiveParticleAcceleration[dim][i]*TimeStep)*coef2;
+	  apvel[dim] = (coef1*apvel[dim] + ActiveParticleAcceleration[dim][i]*TimeStep)*coef2;
 	  
  
 #endif /* VELOCITY_METHOD3 */
 	  
         }
-	for (dim1 = 0; dim1 < GridRank; dim1++)
-	  delete [] ActiveParticleVelocity[dim1];
-	delete [] ActiveParticleVelocity;
+	ActiveParticles[i]->SetVelocity(apvel);
       }
     }
-    else
- 
+    else {
+      
       /* Otherwise, just add the acceleration. */
- 
-    for (i = 0; i < NumberOfParticles; i++)
-      ParticleVelocity[dim][i] += ParticleAcceleration[dim][i] * TimeStep;
-    for (i = 0; i < NumberOfActiveParticles; i++) {
-      float **ActiveParticleVelocity = new float*[GridRank];
-      for (dim1 = 0; dim1 < GridRank; dim1++)
-	ActiveParticleVelocity[dim1] = new float[NumberOfActiveParticles];
 
-      this->GetActiveParticleVelocity(ActiveParticleVelocity);
-
-      ActiveParticleVelocity[dim][i] += ActiveParticleAcceleration[dim][i] * TimeStep;
-
-      for (int dim = 0; dim1 < GridRank; dim1++)
-	delete [] ActiveParticleVelocity[dim1];
-      delete [] ActiveParticleVelocity;
+      for (i = 0; i < NumberOfActiveParticles; i++) {
+	float* apvel = ActiveParticles[i]->ReturnVelocity();
+      
+	for (dim = 0; dim < GridRank; dim++)
+	  apvel[dim] += ActiveParticleAcceleration[dim][i] * TimeStep;
+    
+	ActiveParticles[i]->SetVelocity(apvel);
+      }
     }
   }
 
-
-  if (ProblemType == 29)
-    for (i = 0; i < NumberOfParticles; i++)
-      printf("id=%"PISYM"  %"PSYM" %"PSYM" %"PSYM"\n", ParticleNumber[i],
-	     ParticlePosition[0][i], ParticlePosition[1][i], ParticlePosition[2][i]);
-
- 
   return SUCCESS;
 }

File src/enzo/headers/Grid.h

   /* Create flat arrays of active particle data */
 
   void GetActiveParticlePosition(FLOAT *ActiveParticlePosition[]);
-  void GetActiveParticleVelocity(float *ActiveParticleVelocity[]);
 
   /* Returns averaged velocity from the 6 neighbor cells and itself */