Fix

doc/manual/source/user_guide/InSituPythonAnalysis.rst

@@ -36,17 +36,17 @@ Please follow the instructions in ``libyt`` documents:
 
   * `yt`_: An open-source, permissively-licensed python package for analyzing and visualizing volumetric data.
 
-  * `yt_libyt`_: ``libyt``'s yt frontend.
+  * `yt_libyt`_ (>=0.0.9): ``libyt``'s yt frontend.
 
   * `jupyter_libyt`_: A Jupyter Client plugin for connecting to libyt Jupyter kernel. This is only required in **jupyter kernel mode**.
 
-.. _libyt: https://libyt.readthedocs.io/en/latest/how-to-install.html#c-library-libyt
+.. _libyt: https://libyt.readthedocs.io/en/latest/
 
 .. _yt: https://yt-project.org
 
-.. _yt_libyt: https://libyt.readthedocs.io/en/latest/how-to-install.html#yt-libyt
+.. _yt_libyt: https://github.com/data-exp-lab/yt_libyt
 
-.. _jupyter_libyt: https://libyt.readthedocs.io/en/latest/how-to-install.html#jupyter-libyt
+.. _jupyter_libyt: https://github.com/yt-project/jupyter_libyt
 
 .. _mpi4py: https://mpi4py.readthedocs.io/en/stable/install.html#installation
 
@@ -295,7 +295,7 @@ Please add ``OMPI_MCA_osc=sm,pt2pt`` before ``mpirun``, for example:
 
     OMPI_MCA_osc=sm,pt2pt mpirun -np 4 ./enzo.exe -d CollapseTestNonCosmological.enzo
 
-This is something ``libyt`` will update and improve in the future.
+It is for one-sided MPI communication settings.
 
 .. _Doing In Situ Analysis:
 

src/enzo/ActiveParticle.h

@@ -462,6 +462,68 @@ namespace ActiveParticleHelpers {
 
   }
 
+  template <class APClass> std::vector<hid_t> GetParticleAttributesHDF5DataType() {
+      std::vector<hid_t> attribute_datatype;
+
+      AttributeVector &handlers = APClass::AttributeHandlers;
+      for (AttributeVector::iterator it = handlers.begin(); it != handlers.end(); ++it) {
+          attribute_datatype.push_back((*it)->hdf5type);
+      }
+
+      return attribute_datatype;
+  }
+
+  template <class APClass> std::vector<void*> GetParticleAttributes(ActiveParticleList<ActiveParticleType>& InList,
+          int ParticleTypeID, int TotalParticles, int Count, const std::string& particle_name) {
+
+      // This function allocates and returns the new buffer for particle attributes,
+      // and it is the caller's responsibility to free the buffer when it doesn't need it at all.
+      // It also ignores multi-dimensional arrays.
+      // Currently, this function is solely used in libyt.
+
+      std::vector<void*> attribute_values;
+
+      AttributeVector &handlers = APClass::AttributeHandlers;
+      for (AttributeVector::iterator it = handlers.begin(); it != handlers.end(); ++it) {
+          const char *attr_name = (*it)->name.c_str();
+
+          // ignores multi-dimensional arrays
+          if ((strcmp(attr_name, "AccretionRate") == 0 || strcmp(attr_name, "AccretionRateTime") == 0 || strcmp(attr_name, "Accreted_angmom") == 0)
+              && strcmp(particle_name.c_str(), "SmartStar") == 0) {
+              attribute_values.push_back(nullptr);
+          }
+          // allocates memory and get a copy of the attribute's value
+          else {
+              int size = Count * (*it)->element_size;
+              char *buffer = new char [size];
+              char *_buffer = buffer; // This is because GetAttribute _shifts_ the pointer after returning
+
+              for (int i = 0; i < TotalParticles; i++) {
+                  if (InList[i]->GetEnabledParticleID() != ParticleTypeID) {
+                      continue;
+                  }
+                  (*it)->GetAttribute(&buffer, InList[i]);
+              }
+
+              attribute_values.push_back((void*) _buffer);
+          }
+      }
+
+      return attribute_values;
+  }
+
+  template <class APClass> std::vector<std::string> GetParticleAttributeNames() {
+
+      std::vector<std::string> attributes;
+
+      AttributeVector &handlers = APClass::AttributeHandlers;
+      for (AttributeVector::iterator it = handlers.begin(); it != handlers.end(); ++it) {
+          attributes.push_back((*it)->name);
+      }
+
+      return attributes;
+  }
+
   template <class APClass> void WriteParticles(
           ActiveParticleList<ActiveParticleType>& InList, int ParticleTypeID,
           int TotalParticles,
@@ -687,6 +749,10 @@ class ActiveParticleType_info
    void (*unpack_buffer)(char *buffer, int offset, ActiveParticleList<ActiveParticleType> &Outlist,
                        int OutCount),
    int (*element_size)(void),
+   std::vector<hid_t> (*get_particle_attributes_hdf5_datatype)(void),
+   std::vector<void*> (*get_particle_attributes)(ActiveParticleList<ActiveParticleType> &particles,
+                                                 int type_id, int total_particles, int count, const std::string& particle_name),
+   std::vector<std::string> (*get_particle_attribute_names)(void),
    void (*write_particles)(ActiveParticleList<ActiveParticleType> &particles,
                          int type_id, int total_particles,
                          const std::string name, hid_t node),
@@ -709,6 +775,9 @@ class ActiveParticleType_info
     this->AllocateBuffer = allocate_buffer;
     this->UnpackBuffer = unpack_buffer;
     this->ReturnElementSize = element_size;
+    this->GetParticleAttributesHDF5DataType = get_particle_attributes_hdf5_datatype;
+    this->GetParticleAttributes = get_particle_attributes;
+    this->GetParticleAttributeNames = get_particle_attribute_names;
     this->WriteParticles = write_particles;
     this->ReadParticles = read_particles;
     this->ResetAcceleration = reset_acceleration;
@@ -754,6 +823,11 @@ class ActiveParticleType_info
                        int OutCount);
   int (*FillBuffer)(ActiveParticleList<ActiveParticleType> &InList, int InCount, char *buffer);
   int (*ReturnElementSize)(void);
+  std::vector<hid_t> (*GetParticleAttributesHDF5DataType)(void);
+  std::vector<void*> (*GetParticleAttributes)(ActiveParticleList<ActiveParticleType> &InList,
+                                              int ParticleTypeID, int TotalParticles, int Count,
+                                              const std::string& particle_name);
+  std::vector<std::string> (*GetParticleAttributeNames)(void);
   void (*WriteParticles)(ActiveParticleList<ActiveParticleType> &InList,
                        int ParticleTypeID, int TotalParticles,
                        const std::string name, hid_t node);
@@ -793,6 +867,9 @@ ActiveParticleType_info *register_ptype(std::string name)
      (&ActiveParticleHelpers::FillBuffer<APClass>),
      (&ActiveParticleHelpers::Unpack<APClass>),
      (&ActiveParticleHelpers::CalculateElementSize<APClass>),
+     (&ActiveParticleHelpers::GetParticleAttributesHDF5DataType<APClass>),
+     (&ActiveParticleHelpers::GetParticleAttributes<APClass>),
+     (&ActiveParticleHelpers::GetParticleAttributeNames<APClass>),
      (&ActiveParticleHelpers::WriteParticles<APClass>),
      (&ActiveParticleHelpers::ReadParticles<APClass>),
      (&APClass::ResetAcceleration),

src/enzo/CallInSitulibyt.C

@@ -18,6 +18,7 @@
 
 #include <stdlib.h>
 #include <stdio.h>
+#include <math.h>
 #include "ErrorExceptions.h"
 #include "macros_and_parameters.h"
 #include "typedefs.h"
@@ -43,6 +44,20 @@ void ExposeGridHierarchy(int NumberOfGrids);
 void ExportParameterFile(TopGridData *MetaData, FLOAT CurrentTime, FLOAT OldTime, float dtFixed);
 void CommunicationBarrier();
 
+static yt_dtype MapHDF5TypeToYTType(hid_t hdf5type) {
+    if (hdf5type == HDF5_INT) {
+        return EYT_INT;
+    } else if (hdf5type == HDF5_REAL) {
+        return EYT_BFLOAT;
+    } else if (hdf5type == HDF5_R8) {
+        return YT_DOUBLE;
+    } else if (hdf5type == HDF5_PREC) {
+        return EYT_PFLOAT;
+    } else {
+        return YT_DTYPE_UNKNOWN;
+    }
+}
+
 #endif
 
 int CallInSitulibyt(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
@@ -149,10 +164,11 @@ int CallInSitulibyt(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
     params->length_unit = LengthUnits;
     params->mass_unit = DensityUnits * LengthUnits * LengthUnits * LengthUnits; /* this right? */
     params->time_unit = TimeUnits;
-    params->magnetic_unit = 0.0; /* Not right */
-    params->periodicity[0] = 1; /* also not right */
-    params->periodicity[1] = 1;
-    params->periodicity[2] = 1;
+    params->velocity_unit = VelocityUnits;
+    params->magnetic_unit = sqrt(4.0 * 3.141592653589793238462643383279502884L * DensityUnits) * VelocityUnits;
+    params->periodicity[0] = MetaData->LeftFaceBoundaryCondition[0];
+    params->periodicity[1] = MetaData->LeftFaceBoundaryCondition[1];
+    params->periodicity[2] = MetaData->LeftFaceBoundaryCondition[2];
     params->dimensionality = MetaData->TopGridRank;
     params->domain_dimensions[0] = MetaData->TopGridDims[0];
     params->domain_dimensions[1] = MetaData->TopGridDims[1];
@@ -170,23 +186,44 @@ int CallInSitulibyt(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
         }
     }
 
-    params->num_par_types = 1;
-    yt_par_type par_type_list[1];
-    par_type_list[0].par_type = "io";
+    // Add fields for Temperature/Cooling_Time derived field from enzo
+    params->num_fields += 2;
+
+    // Add active particle ptypes
+    params->num_par_types = 1 + EnabledActiveParticlesCount; // DarkMatter and Other ActiveParticle (ex: SmartStar)
+    yt_par_type *par_type_list = new yt_par_type [params->num_par_types];
+
+    par_type_list[0].par_type = "DarkMatter";
     par_type_list[0].num_attr = 3 + 3 + 1 + 1 + 1 + NumberOfParticleAttributes;
+
+    // the attributes name should be alive within the entire libyt in situ analysis,
+    // because libyt does not make a copy of the names.
+    std::vector<std::vector<std::string>> active_particles_attributes;
+    std::vector<std::vector<hid_t>> active_particles_attributes_hdf5type;
+
+    for (int i = 0; i < EnabledActiveParticlesCount; i++) {
+        ActiveParticleType_info *ActiveParticleTypeToEvaluate = EnabledActiveParticles[i];
+        active_particles_attributes.emplace_back(ActiveParticleTypeToEvaluate->GetParticleAttributeNames());
+        active_particles_attributes_hdf5type.emplace_back(ActiveParticleTypeToEvaluate->GetParticleAttributesHDF5DataType());
+
+        par_type_list[1 + i].par_type = ActiveParticleTypeToEvaluate->particle_name.c_str();
+        par_type_list[1 + i].num_attr = active_particles_attributes[i].size();
+    }
+
     params->par_type_list = par_type_list;
 
     if (yt_set_Parameters(params) != YT_SUCCESS){
         fprintf(stderr, "Error in libyt API yt_set_Parameters\n");
         return FAIL;
     }
 
+    delete [] par_type_list;
+
     yt_particle *particle_list;
     yt_get_ParticlesPtr(&particle_list);
 
-    // Particle type "io", each particle has position in the center of the grid it belongs to with value grid level.
-    // par_type and num_attr will be assigned by libyt with the same value we passed in par_type_list at yt_set_Parameters.
-    particle_list[0].par_type = "io";
+    // TODO: make sure enzo's particle is always DarkMatter
+    particle_list[0].par_type = "DarkMatter";
 
     // We have the attributes: 3 positions, 3 velocities, "mass", ID and Type
     // and extras.
@@ -226,6 +263,19 @@ int CallInSitulibyt(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
     particle_list[0].coor_y = attr_name[1];
     particle_list[0].coor_z = attr_name[2];
 
+    // we need to map hdf5type to yt datatype,
+    // since we want to avoid storing yt data type in ParticleAttributeHandler class and causing
+    for (int i = 0; i < EnabledActiveParticlesCount; i++) {
+        for (int v = 0; v < active_particles_attributes[i].size(); v++) {
+            particle_list[1 + i].attr_list[v].attr_name = active_particles_attributes[i][v].c_str();
+            particle_list[1 + i].attr_list[v].attr_dtype = MapHDF5TypeToYTType(active_particles_attributes_hdf5type[i][v]);
+        }
+
+        particle_list[1 + i].coor_x = "particle_position_x";
+        particle_list[1 + i].coor_y = "particle_position_y";
+        particle_list[1 + i].coor_z = "particle_position_z";
+    }
+
     /* Set code-specific parameter */
     char tempname[256];
     #include "InitializeLibytInterface_finderfunctions.inc"
@@ -277,14 +327,30 @@ int CallInSitulibyt(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
     /* Now we add on the following fields, as per grid::WriteGrid:
      *
      *  - Temperature
-     *  - Dust_Temperature
+     *  - Dust_Temperature (ignore for now)
      *  - Cooling_Time
      *
      *  Each of these is predicated on the global parameter associated with
      *  them.
      *
      * */
 
+    field_list[libyt_field_i].field_name = "Temperature";
+    field_list[libyt_field_i].field_type = "cell-centered";
+    field_list[libyt_field_i].field_dtype = EYT_BFLOAT;
+    for (j = 0; j < 6; j++) {
+        field_list[libyt_field_i].field_ghost_cell[j] = NumberOfGhostZones;
+    }
+    libyt_field_i = libyt_field_i + 1;
+
+    field_list[libyt_field_i].field_name = "Cooling_Time";
+    field_list[libyt_field_i].field_type = "cell-centered";
+    field_list[libyt_field_i].field_unit = "code_time";
+    field_list[libyt_field_i].field_dtype = EYT_BFLOAT;
+    for (j = 0; j < 6; j++) {
+        field_list[libyt_field_i].field_ghost_cell[j] = NumberOfGhostZones;
+    }
+
     /* We now have to do everything we do in CallPython.C, which amounts to
      *
      *  - ExposeGridHierarchy (not necessary anymore)
@@ -371,6 +437,11 @@ int CallInSitulibyt(LevelHierarchyEntry *LevelArray[], TopGridData *MetaData,
         return FAIL;
     }
 
+    for (std::size_t i = 0; i < libyt_generated_data.size(); i++) {
+        delete [] libyt_generated_data[i];
+    }
+    libyt_generated_data.clear();
+
     CommunicationBarrier();
     return SUCCESS;
 #endif

src/enzo/Grid_ConvertToLibyt.C

@@ -17,6 +17,7 @@
 #include <map>
 #include <stdlib.h>
 #include <stdio.h>
+#include <math.h>
 #include <unistd.h>
 #include "libyt.h"
 #include "ErrorExceptions.h"
@@ -29,6 +30,8 @@
 #include "Grid.h"
 #include "CosmologyParameters.h"
 
+int GetUnits(float *DensityUnits, float *LengthUnits, float *TemperatureUnits, float *TimeUnits, float *VelocityUnits, FLOAT Time);
+
 /* Note that previously we allowed WriteTime to be supplied and we do not now */
 void grid::ConvertToLibyt(int LocalGridID, int GlobalGridID, int ParentID, int level, yt_grid &GridInfo)
 {
@@ -78,10 +81,30 @@ void grid::ConvertToLibyt(int LocalGridID, int GlobalGridID, int ParentID, int l
          *
          * */
         int libyt_field = libyt_field_lookup[field];
-        if (libyt_field == -1) continue;
+        if (libyt_field == -1) continue; // TODO: will this ever be -1?
         GridInfo.field_data[libyt_field].data_ptr = BaryonField[field];
     }
 
+    long grid_size = GridDimension[0] * GridDimension[1] * GridDimension[2];
+    float *temp_field = new float[grid_size];
+    float *cooling_time_field = new float[grid_size];
+
+    float TemperatureUnits = 1, DensityUnits = 1, LengthUnits = 1, VelocityUnits = 1, TimeUnits = 1, aUnits = 1;
+    GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits, &VelocityUnits, this->Time);
+
+    if (this->ComputeTemperatureField(temp_field) == SUCCESS) {
+        libyt_generated_data.push_back(temp_field);
+        int temp_field_i = ((yt_param_yt*)param_yt)->num_fields - 2;
+        GridInfo.field_data[temp_field_i].data_ptr = temp_field;
+    }
+
+    if (this->ComputeCoolingTime(cooling_time_field) == SUCCESS) {
+        for (long i = 0; i < grid_size; i++) { cooling_time_field[i] = fabs(cooling_time_field[i]) * TimeUnits; }
+        libyt_generated_data.push_back(cooling_time_field);
+        int cooling_time_field_i = ((yt_param_yt*)param_yt)->num_fields - 1;
+        GridInfo.field_data[cooling_time_field_i].data_ptr = cooling_time_field;
+    }
+
     /* par_count_list can take multiple particle types */
     if (this->ReturnNumberOfParticles() > 0) {
         GridInfo.par_count_list[0] = this->ReturnNumberOfParticles();
@@ -100,5 +123,35 @@ void grid::ConvertToLibyt(int LocalGridID, int GlobalGridID, int ParentID, int l
     else {
         GridInfo.par_count_list[0] = 0;
     }
+
+    /* fill in active particle count */
+    std::vector<int> active_particle_count(EnabledActiveParticlesCount, 0);
+    for (int i = 0; i < NumberOfActiveParticles; i++) {
+        active_particle_count[(this->ActiveParticles)[i]->GetEnabledParticleID()]++;
+    }
+    for (int i = 0; i < EnabledActiveParticlesCount; i++) {
+        GridInfo.par_count_list[1 + i] = active_particle_count[i];
+    }
+
+    /* fill in buffer only if there is active particle, which is sum of active particle count > 0
+     * ignore AccretionRate/AccretionRateTime/Accreted_angmom for now since they are not one dimensional */
+    if (NumberOfActiveParticles > 0) {
+        for (int i = 0; i < EnabledActiveParticlesCount; i++) {
+            ActiveParticleType_info *ActiveParticleTypeToEvaluate = EnabledActiveParticles[i];
+
+            // the returned buffer ignores multi-dimensional array for now, it set them as nullptr.
+            std::vector<void*> par_attr_buffer = ActiveParticleTypeToEvaluate->GetParticleAttributes(
+                    this->ActiveParticles, i, NumberOfActiveParticles, active_particle_count[i],
+                    ActiveParticleTypeToEvaluate->particle_name
+                    );
+
+            for (int a = 0; a < par_attr_buffer.size(); a++) {
+                GridInfo.particle_data[1 + i][a].data_ptr = par_attr_buffer[a];
+            }
+
+            // free these pre-allocated buffer after libyt in situ analysis is done.
+            libyt_generated_data.insert(libyt_generated_data.end(), par_attr_buffer.begin(), par_attr_buffer.end());
+        }
+    }
 }
 #endif

src/enzo/InitializeLibytInterface_finderfunctions.inc

@@ -33,7 +33,7 @@ YT_SET_USERPARAM_INTEGER("ProblemType", 1, &ProblemType);
 YT_SET_USERPARAM_INTEGER("HydroMethod", 1, &HydroMethod);
 YT_SET_USERPARAM_NONINTEGER("huge_number", 1, &huge_number);
 YT_SET_USERPARAM_NONINTEGER("tiny_number", 1, &tiny_number);
-YT_SET_USERPARAM_NONINTEGER("Gamma", 1, &Gamma);
+YT_SET_USERPARAM_NONINTEGER("gamma", 1, &Gamma);
 YT_SET_USERPARAM_INTEGER("PressureFree", 1, &PressureFree);
 YT_SET_USERPARAM_INTEGER("QuantumPressure", 1, &QuantumPressure);
 YT_SET_USERPARAM_NONINTEGER("FDMMass", 1, &FDMMass);