Particle masses bugfix

src/Cello/data_Particle.hpp

@@ -125,6 +125,11 @@ class Particle {
   int num_attributes(int it) const
   { return particle_descr_->num_attributes(it); }
 
+  /// Check if attribute exists
+
+  bool has_attribute (int it, std::string attribute) const
+  { return particle_descr_->has_attribute(it,attribute); }
+
   /// Return the index for the given attribute
 
   int attribute_index (int it, std::string attribute) const
@@ -197,6 +202,10 @@ class Particle {
   int constant_offset(int it, int ic) const
   { return particle_descr_->constant_offset(it,ic); }
 
+  /// Check if given constant exists
+
+  bool has_constant(int it, std::string constant) const
+  { return particle_descr_->has_constant(it,constant); }
 
   //--------------------------------------------------
   // BYTES

src/Cello/data_ParticleDescr.cpp

@@ -253,6 +253,24 @@ int ParticleDescr::attribute_index (int it, std::string attribute_name) const
 
 //----------------------------------------------------------------------
 
+bool ParticleDescr::has_attribute(int it, std::string attribute_name) const
+{
+  ASSERT1("ParticleDescr::is_attribute",
+          "Trying to access unknown particle type '%s'",
+          type_name(it).c_str(),
+          (0 <= it && it < num_types()));
+
+  auto iter=attribute_index_[it].find(attribute_name);
+
+  int index = (iter != attribute_index_[it].end()) ? iter->second : -1;
+
+  static int count[CONFIG_NODE_SIZE] = {0};
+
+  return !((index == -1) && (count[cello::index_static()]++ < 10));
+}
+
+//----------------------------------------------------------------------
+
 std::string ParticleDescr::attribute_name (int it, int ia) const
 {
   ASSERT2("ParticleDescr::attribute_name",
@@ -390,6 +408,23 @@ int ParticleDescr::constant_offset (int it, int ic) const
   return constant_offset_[it][ic]; 
 }
 
+//----------------------------------------------------------------------
+
+bool ParticleDescr::has_constant(int it, std::string constant_name) const
+{
+  ASSERT1("ParticleDescr::constant_index",
+          "Trying to access unknown particle type '%s'",
+          type_name(it).c_str(),
+          (0 <= it && it < num_types()));
+
+  auto iter=constant_index_[it].find(constant_name);
+
+  bool check = (iter != constant_index_[it].end()) ? true : false;
+
+  return check;
+}
+
+
 //----------------------------------------------------------------------
 // Bytes
 //----------------------------------------------------------------------

src/Cello/data_ParticleDescr.hpp

@@ -82,6 +82,10 @@ class ParticleDescr {
    /// Return a pointer to the given constant for the given type
   char * constant_value (int it, int ic);
 
+  /// Check if given constant exists for the given particle type
+
+  bool has_constant (int it, std::string constant) const;
+
   //--------------------------------------------------
   // ATTRIBUTES
   //--------------------------------------------------
@@ -94,6 +98,10 @@ class ParticleDescr {
 
   int num_attributes(int it) const;
 
+  /// Check if attribute exists
+
+  bool has_attribute (int it, std::string attribute) const;
+
   /// Return the index for the given attribute
 
   int attribute_index (int it, std::string attribute) const;

src/Enzo/enzo_EnzoMethodPmDeposit.cpp

@@ -2,6 +2,7 @@
 
 /// @file     enzo_EnzoMethodPmDeposit.cpp
 /// @author   James Bordner (jobordner@ucsd.edu)
+/// @author   Stefan Arridge (stefan.arridge@gmail.com)
 /// @date     Fri Apr  2 17:05:23 PDT 2010
 /// @brief    Implements the EnzoMethodPmDeposit class
 ///
@@ -104,22 +105,38 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
     block->lower(&xm,&ym,&zm);
     block->upper(&xp,&yp,&zp);
     block->cell_width(&hx,&hy,&hz);
+
 
-    // declare particle position arrays
-
-    //--------------------------------------------------
-    // Add "dark" particle density
-    //--------------------------------------------------
+    // Get cosmological scale factors, if cosmology is turned on
+    enzo_float cosmo_a=1.0;
+    enzo_float cosmo_dadt=0.0;
+    EnzoPhysicsCosmology * cosmology = enzo::cosmology();
 
-    if (particle.type_exists ("dark")) {
-      const int it = particle.type_index ("dark");
+    if (cosmology) {
+      cosmology->compute_expansion_factor(&cosmo_a,&cosmo_dadt,
+					  block->time() + alpha_*block->dt());
+    }
 
-      const int ic_mass = particle.constant_index (it,"mass");
+    double inv_vol = 1.0;
+    if (rank >= 1) inv_vol /= hx;
+    if (rank >= 2) inv_vol /= hy;
+    if (rank >= 3) inv_vol /= hz;
 
-      enzo_float dens = *((enzo_float *)(particle.constant_value (it,ic_mass)));
+    const double dt = alpha_ * block->dt() / cosmo_a;
 
-      // check precisions match
-
+    // Get the number of particle types in the 'has_mass' group
+    ParticleDescr * particle_descr = cello::particle_descr();
+    Grouping * particle_groups = particle_descr->groups();
+    int num_mass = particle_groups->size("has_mass");
+
+    // Loop over all particles that have mass
+    for (int ipt = 0; ipt < num_mass; ipt++) {
+      const int it = particle.type_index(particle_groups->item("has_mass",ipt));
+
+      // Index for mass attribute / constant
+      int im = 0;
+
+      // check correct precision for position
       int ia = particle.attribute_index(it,"x");
       int ba = particle.attribute_bytes(it,ia); // "bytes (actual)"
       int be = sizeof(enzo_float);                // "bytes (expected)"
@@ -132,40 +149,51 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
                ((be == 4) ? "single" : ((be == 8) ? "double" : "quadruple")),
                (ba == be));
 
-      // Accumulate particle density using CIC
-
-      enzo_float cosmo_a=1.0;
-      enzo_float cosmo_dadt=0.0;
-      EnzoPhysicsCosmology * cosmology = enzo::cosmology();
-
-      if (cosmology) {
-
-        double time = block->time();
-        double dt   = block->dt();
-        cosmology->compute_expansion_factor (&cosmo_a,&cosmo_dadt,time+alpha_*dt);
-
-      }
-      if (rank >= 1) hx *= cosmo_a;
-      if (rank >= 2) hy *= cosmo_a;
-      if (rank >= 3) hz *= cosmo_a;
-
-      const double dt = alpha_ * block->dt() / cosmo_a;
-
-      // Scale mass by volume if particle value is mass instead of density
 
-      int level = block->level();
-
-      // Required for Cosmology ("mass" is mass)
-      // Not for Collapse ("mass" is density)
-      if (cosmology) {
-        dens *= std::pow(2.0,rank*level);
-      }
+      // Loop over batches
+      for (int ib=0; ib<particle.num_batches(it); ib++) {
+	// Pointer to array of particle masses
+	enzo_float * pmass = NULL;
+        const int np = particle.num_particles(it,ib);
 
-      // Accumulated single velocity array for Baryon deposit
+	// Particle masses can be either a constant value for all particles,
+	// or an attribrute array
+	if (particle.has_constant (it,"mass")){
 
-      for (int ib=0; ib<particle.num_batches(it); ib++) {
+	  // Check if particle also has a mass attribute
+	  ASSERT1("EnzoMethodPmDeposit::compute",
+		  "Particle type %s has both a constant and an attribute called"
+		  "'mass'. Exiting.", particle.type_name(it).c_str(),
+		  particle.has_attribute(it,"mass"));
+
+	  // In this case we fill an array of length np with the constant
+	  // value
+	  pmass = new enzo_float[np];
+	  im = particle.constant_index(it,"mass");
+	  for (int ip = 0; ip<np; ip++)
+	    pmass[ip] = *((enzo_float *)(particle.constant_value (it,im)));
+
+        } else if (particle.has_attribute(it,"mass")) {
 
-        const int np = particle.num_particles(it,ib);
+	  // In this case we set to pointer to the attribute array
+	  im = particle.attribute_index(it,"mass");
+	  pmass = (enzo_float *) particle.attribute_array( it, im, ib);
+        }
+	else {
+	  ERROR1("EnzoMethodPmDeposit::compute",
+		"Particle type %s has neither a constant nor an attribute" 
+                "called 'mass', yet is in the 'has_mass' group",
+		 particle.type_name(it).c_str());
+	}
+
+	// If mass is a constant, we simply loop through the pmass
+	// array, and so dm = 1. If its an attribute, we need to get
+	// the stride length
+	const int dm = particle.has_attribute(it,"mass") ?
+	               particle.stride(it,im) : 1;
+
+
+	// Allocate particle masses to the grid with CIC scheme
 
         if (rank == 1) {
 
@@ -180,6 +208,9 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
 
           for (int ip=0; ip<np; ip++) {
 
+	    // Stefan: Unsure about this line, doesn't match what's in the
+	    // Enzo paper, which says that particles are drifted by half a timestep
+	    // before being deposited to the grid
             double x = xa[ip*dp] + vxa[ip*dv]*dt;
 
             double tx = nx*(x - xm) / (xp - xm) - 0.5;
@@ -191,8 +222,8 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
             double x0 = 1.0 - (tx - floor(tx));
             double x1 = 1.0 - x0;
 
-            de_p[ix0] += dens*x0;
-            de_p[ix1] += dens*x1;
+            de_p[ix0] += pmass[ip*dm]*inv_vol*x0;
+            de_p[ix1] += pmass[ip*dm]*inv_vol*x1;
 
           }
 
@@ -215,7 +246,7 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
 
           const int dp =  particle.stride(it,ia_x);
           const int dv =  particle.stride(it,ia_vx);
-
+ 
           for (int ip=0; ip<np; ip++) {
 
             double x = xa[ip*dp] + vxa[ip*dv]*dt;
@@ -236,30 +267,31 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
             double x1 = 1.0 - x0;
             double y1 = 1.0 - y0;
 
-            if ( dens < 0.0) {
-              CkPrintf ("%s:%d ERROR: dens = %f\n", __FILE__,__LINE__,dens);
+            if ( pmass[ip*dm]*inv_vol < 0.0) {
+              CkPrintf ("%s:%d ERROR: pmass[%d]*inv_vol = %f\n",
+			__FILE__,__LINE__,ip*dm,pmass[ip*dm]*inv_vol);
             }
 
-            de_p[ix0+mx*iy0] += dens*x0*y0;
-            de_p[ix1+mx*iy0] += dens*x1*y0;
-            de_p[ix0+mx*iy1] += dens*x0*y1;
-            de_p[ix1+mx*iy1] += dens*x1*y1;
+            de_p[ix0+mx*iy0] += pmass[ip*dm]*inv_vol*x0*y0;
+            de_p[ix1+mx*iy0] += pmass[ip*dm]*inv_vol*x1*y0;
+            de_p[ix0+mx*iy1] += pmass[ip*dm]*inv_vol*x0*y1;
+            de_p[ix1+mx*iy1] += pmass[ip*dm]*inv_vol*x1*y1;
 
             if ( de_p[ix0+mx*iy0] < 0.0) {
               CkPrintf ("%s:%d ERROR: de_p %d %d = %f\n",
-                        __FILE__,__LINE__,ix0,iy0,dens);
+                        __FILE__,__LINE__,ix0,iy0,pmass[ip*dm]*inv_vol);
             }
             if ( de_p[ix1+mx*iy0] < 0.0) {
               CkPrintf ("%s:%d ERROR: de_p %d %d = %f\n",
-                        __FILE__,__LINE__,ix1,iy0,dens);
+                        __FILE__,__LINE__,ix1,iy0,pmass[ip*dm]*inv_vol);
             }
             if ( de_p[ix0+mx*iy1] < 0.0) {
               CkPrintf ("%s:%d ERROR: de_p %d %d = %f\n",
-                        __FILE__,__LINE__,ix0,iy1,dens);
+                        __FILE__,__LINE__,ix0,iy1,pmass[ip*dm]*inv_vol);
             }
             if ( de_p[ix1+mx*iy1] < 0.0) {
               CkPrintf ("%s:%d ERROR: de_p %d %d = %f\n",
-                        __FILE__,__LINE__,ix1,iy1,dens);
+                        __FILE__,__LINE__,ix1,iy1,pmass[ip*dm]*inv_vol);
             }
           }
 
@@ -318,19 +350,23 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
             double y1 = 1.0 - y0;
             double z1 = 1.0 - z0;
 
-            de_p[ix0+mx*(iy0+my*iz0)] += dens*x0*y0*z0;
-            de_p[ix1+mx*(iy0+my*iz0)] += dens*x1*y0*z0;
-            de_p[ix0+mx*(iy1+my*iz0)] += dens*x0*y1*z0;
-            de_p[ix1+mx*(iy1+my*iz0)] += dens*x1*y1*z0;
-            de_p[ix0+mx*(iy0+my*iz1)] += dens*x0*y0*z1;
-            de_p[ix1+mx*(iy0+my*iz1)] += dens*x1*y0*z1;
-            de_p[ix0+mx*(iy1+my*iz1)] += dens*x0*y1*z1;
-            de_p[ix1+mx*(iy1+my*iz1)] += dens*x1*y1*z1;
-
-          }
-        }
-      }
-    }
+            de_p[ix0+mx*(iy0+my*iz0)] += pmass[ip*dm]*inv_vol*x0*y0*z0;
+            de_p[ix1+mx*(iy0+my*iz0)] += pmass[ip*dm]*inv_vol*x1*y0*z0;
+            de_p[ix0+mx*(iy1+my*iz0)] += pmass[ip*dm]*inv_vol*x0*y1*z0;
+            de_p[ix1+mx*(iy1+my*iz0)] += pmass[ip*dm]*inv_vol*x1*y1*z0;
+            de_p[ix0+mx*(iy0+my*iz1)] += pmass[ip*dm]*inv_vol*x0*y0*z1;
+            de_p[ix1+mx*(iy0+my*iz1)] += pmass[ip*dm]*inv_vol*x1*y0*z1;
+            de_p[ix0+mx*(iy1+my*iz1)] += pmass[ip*dm]*inv_vol*x0*y1*z1;
+            de_p[ix1+mx*(iy1+my*iz1)] += pmass[ip*dm]*inv_vol*x1*y1*z1;
+
+          } // Loop over particle in batch
+        } // if rank == 3
+
+	// If constant mass, delete the pmass array
+	if (particle.has_constant(it,"mass")) delete [] pmass;
+
+      } // Loop over batches
+    } // Loop over particle types in 'has_mass' group
 
     //--------------------------------------------------
     // Add gas density