Commits

Matthew Turk committed c1f7ae3 Merge

Automated merge with file:///Users/matthewturk/enzo/week-of-code

  • Participants
  • Parent commits cf28d9e, 1c4d8c1
  • Branches trunk

Comments (0)

Files changed (12)

File src/enzo/EvolveHierarchy.C

 				TopGridData* MetaData = NULL);
 double ReturnWallTime(void);
 int Enzo_Dims_create(int nnodes, int ndims, int *dims);
+#ifdef EMBEDDED_PYTHON
+int CallPython();
+#endif
 
 #ifdef MEM_TRACE
 Eint64 mused(void);

File src/enzo/EvolveLevel.C

 /                Optional StaticSiblingList for root grid
 /  modified8:  April, 2009 by John Wise
 /                Added star particle class and radiative transfer
+/  modified9:  June, 2009 by Matthew Turk
+/                Added Python
 /
 /  PURPOSE:
 /    This routine is the main grid evolution function.  It assumes that the
 static int StaticLevelZero = 0;
 #endif
 
+#ifdef EMBEDDED_PYTHON
+int CallPython(TopGridData &MetaData, FLOAT CurrentTime);
+int ExposeDataHierarchy(TopGridData &MetaData, HierarchyEntry *Grid, 
+		       int &GridID, FLOAT WriteTime, int reset, int ParentID, int level);
+void ExposeGridHierarchy(int NumberOfGrids);
+#endif
 #define TIME_MESSAGING 
 
 
 #endif
 
     }
+
+#ifdef EMBEDDED_PYTHON
+#ifdef USE_MPI
+    MPI_Barrier(MPI_COMM_WORLD);
+#endif
+    if ( LevelArray[level+1] == NULL) {
+    //if ( level == 0) {
+#ifdef USE_MPI
+    MPI_Barrier(MPI_COMM_WORLD);
+#endif
+        PyDict_Clear(grid_dictionary);
+        PyDict_Clear(old_grid_dictionary);
+
+        LevelHierarchyEntry *Temp2 = LevelArray[0];
+        /* Count the grids */
+        int num_grids = 0;
+        /* I think there is a better idiom for this somewhere
+           but I couldn't find it, and I think this works     */
+        int start_index = 1;
+        FLOAT WriteTime;
+        for (int lc = 0; LevelArray[lc] != NULL; lc++){
+            Temp2 = LevelArray[lc];
+            while (Temp2 != NULL) {
+                num_grids++; Temp2 = Temp2->NextGridThisLevel;
+            }
+        }
+        ExposeGridHierarchy(num_grids);
+/*
+        if (ExposeGridHierarchy(num_grids) == FAIL) {
+            fprintf(stderr, "Error in ExposeGridHierarchy\n");
+            return FAIL;
+        }
+*/
+        Temp2 = LevelArray[0];
+        if(Temp2->NextGridThisLevel == NULL){fprintf(stderr,"NULLness\n");}
+        while (Temp2->NextGridThisLevel != NULL)
+            Temp2 = Temp2->NextGridThisLevel; /* ugh: find last in linked list */
+        WriteTime = LevelArray[level]->GridData->ReturnTime();
+        fprintf(stderr, "Times: %"GOUTSYM" %"GOUTSYM"\n", WriteTime, MetaData->Time);
+        if (ExposeDataHierarchy(*MetaData, Temp2->GridHierarchyEntry, start_index,
+                    WriteTime, 1, 0, 0) == FAIL) {
+            fprintf(stderr, "Error in ExposeDataHierarchy\n");
+            return FAIL;
+        }
+        CallPython(*MetaData, WriteTime);
+    }
+#endif
  
     /* Check for new output from new file in file system, but only
        when at bottom of hierarchy */

File src/enzo/ExposeDataHierarchy.C

+/***********************************************************************
+/
+/  CONVERT DATA HIERARCHY TO NUMPY FORMAT
+/
+/  written by: Matthew Turk
+/  date:       September, 2008
+/  modified1:
+/
+/  PURPOSE:
+/
+************************************************************************/
+
+// This function writes out the data hierarchy (TopGrid)
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.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"
+
+/* function prototypes */
+
+
+
+int ExposeDataHierarchy(TopGridData &MetaData, HierarchyEntry *Grid, 
+		       int &GridID, FLOAT WriteTime, int reset, int ParentID, int level)
+{
+
+#ifdef EMBEDDED_PYTHON
+  int OriginalID, NextGridThisLevelID, NextGridNextLevelID;
+  int flagged, noParent = 0;
+  static PyArrayObject *container[11];
+  
+  /* Borrowed reference */
+  if(reset == 1) {
+        int i = 0;
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridDimensions");
+        container[i++] =  (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridStartIndices");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridEndIndices");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridLeftEdge");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridRightEdge");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridLevels");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridTimes");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridOldTimes");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridProcs");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridNumberOfParticles");
+        container[i++] = (PyArrayObject *)
+            PyDict_GetItemString(hierarchy_information, "GridParentIDs");
+      }
+
+  /* Write out header info for this grid */
+
+  OriginalID = GridID;
+
+  /* We should have a call for interpolation here */
+
+  Grid->GridData->ConvertToNumpy(GridID, container, ParentID, level, WriteTime);
+
+  /* Write out pointer information for the next grid this level */
+
+  NextGridThisLevelID = GridID + 1;
+  if (Grid->NextGridThisLevel == NULL) NextGridThisLevelID = 0;
+
+  if (NextGridThisLevelID != 0) {
+    GridID++;
+    if (ExposeDataHierarchy(MetaData, Grid->NextGridThisLevel, GridID,
+			   WriteTime, 0, ParentID, level) == FAIL) {
+      fprintf(stderr, "Error in ExposeDataHierarchy(1).\n");
+      return FAIL;
+    }
+  }
+
+  /* Write out pointer information for the next grid next level */
+
+  NextGridNextLevelID = GridID + 1;
+  if (Grid->NextGridNextLevel == NULL) NextGridNextLevelID = 0;
+
+  if (NextGridNextLevelID != 0) {
+    GridID++;
+    if (ExposeDataHierarchy(MetaData, Grid->NextGridNextLevel, GridID,
+			   WriteTime, 0, OriginalID, level+1) == FAIL) {
+      fprintf(stderr, "Error in ExposeDataHierarchy(1).\n");
+      return FAIL;
+    }
+  }
+
+#endif
+  return SUCCESS;
+}

File src/enzo/ExposeGridHierarchy.C

+/***********************************************************************
+/
+/  GENERATE GRID HIERARCHY HANDLERS
+/
+/  written by: Matthew Turk
+/  date:       September, 2008
+/  modified1:
+/
+/  PURPOSE:
+/
+************************************************************************/
+
+// This function writes out the data hierarchy (TopGrid)
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.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"
+
+void ExposeGridHierarchy(int NumberOfGrids)
+{
+#undef int
+
+  /* This function just fills the dictionaries */
+
+#ifdef EMBEDDED_PYTHON
+
+  // We'll set up some common variables
+  npy_intp flat_dimensions[2];
+  PyArrayObject *temp_array;
+
+  PyDict_Clear(hierarchy_information);
+  flat_dimensions[0] = (npy_intp) NumberOfGrids;
+
+  /* Here are the arrays yt expects from the data file:
+    GridDimensions[:] = harray[:,0:3]
+    GridStartIndices[:] = harray[:,3:6]
+    GridEndIndices[:] = harray[:,6:9]
+    GridLeftEdge[:] = harray[:,9:12]
+    GridRightEdge[:] = harray[:,12:15]
+    GridLevels[:] = harray[:,15:16]
+    GridTimes[:] = harray[:,16:17]
+    GridNumberOfParticles[:] = harray[:,17:18] */
+
+  int counter=0;
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 0
+  flat_dimensions[1] = 3;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, NPY_INT);
+  PyDict_SetItemString(hierarchy_information, "GridDimensions", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 1
+  flat_dimensions[1] = 3;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, NPY_INT);
+  PyDict_SetItemString(hierarchy_information, "GridStartIndices", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 2
+  flat_dimensions[1] = 3;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, NPY_INT);
+  PyDict_SetItemString(hierarchy_information, "GridEndIndices", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 3
+  flat_dimensions[1] = 3;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, ENPY_FLOAT);
+  PyDict_SetItemString(hierarchy_information, "GridLeftEdge", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 4
+  flat_dimensions[1] = 3;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, ENPY_FLOAT);
+  PyDict_SetItemString(hierarchy_information, "GridRightEdge", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 5
+  flat_dimensions[1] = 1; /* a bit iffy */
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, NPY_INT);
+  PyDict_SetItemString(hierarchy_information, "GridLevels", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 6
+  flat_dimensions[1] = 1;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, ENPY_FLOAT);
+  PyDict_SetItemString(hierarchy_information, "GridTimes", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 7
+  flat_dimensions[1] = 1;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, ENPY_FLOAT);
+  PyDict_SetItemString(hierarchy_information, "GridOldTimes", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 8
+  flat_dimensions[1] = 1;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, NPY_INT);
+  PyDict_SetItemString(hierarchy_information, "GridProcs", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 9
+  flat_dimensions[1] = 1;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, NPY_INT);
+  PyDict_SetItemString(hierarchy_information, "GridNumberOfParticles", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+  //fprintf(stderr, "counter: %d\n", counter++); // 10
+  flat_dimensions[1] = 1;
+  temp_array = (PyArrayObject *) PyArray_SimpleNew(2, flat_dimensions, NPY_INT);
+  PyDict_SetItemString(hierarchy_information, "GridParentIDs", (PyObject *) temp_array);
+  Py_DECREF(temp_array);
+
+#endif
+  //return SUCCESS;
+}

File src/enzo/Grid.h

 #   include "Grid_AnalyzeClusters.h"
 #endif
 
+#ifdef EMBEDDED_PYTHON
+    void ConvertToNumpy(int GridID, PyArrayObject *container[],
+                        int ParentID, int level, FLOAT WriteTime);
+#endif
 //------------------------------------------------------------------------
 // Methods for star formation
 //------------------------------------------------------------------------

File src/enzo/Grid_ConvertToNumpy.C

+#ifdef EMBEDDED_PYTHON
+/***********************************************************************
+/
+/  GRID CLASS (TURN ALL OF OUR FIELDS INTO NUMPY ARRAYS)
+/
+/  written by: Matthew Turk
+/  date:       September, 2008
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    ENZO_SUCCESS or FAIL
+/
+************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.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 "CosmologyParameters.h"
+
+void grid::ConvertToNumpy(int GridID, PyArrayObject *container[], int ParentID, int level, FLOAT WriteTime)
+{
+
+    this->DebugCheck("Converting to NumPy arrays");
+
+    PyObject *grid_data, *old_grid_data, *field_name, *grid_id;
+
+    /* Declarations */
+
+  /* Only set up baryonfields if it's on our proc */
+  
+    if (ProcessorNumber == MyProcessorNumber) {
+
+        int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num;
+        int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
+            DINum, DIINum, HDINum;
+        int field;
+        FLOAT a = 1.0, dadt;
+
+        /* Find fields: density, total energy, velocity1-3. */
+
+        if (this->IdentifyPhysicalQuantities(DensNum, GENum, Vel1Num, Vel2Num, 
+                    Vel3Num, TENum) == FAIL) {
+            fprintf(stderr, "Error in IdentifyPhysicalQuantities.\n");
+            return ;
+        }
+
+        /* Find Multi-species fields. */
+
+        if (MultiSpecies)
+            if (IdentifySpeciesFields(DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, 
+                        HMNum, H2INum, H2IINum, DINum, DIINum, HDINum) == FAIL) {
+                fprintf(stderr, "Error in grid->IdentifySpeciesFields.\n");
+                return ;
+            }
+
+#undef int
+        grid_data = PyDict_New();
+        old_grid_data = PyDict_New();
+        PyArrayObject *dataset;
+        int nd = 3;
+        npy_intp dims[3];
+        dims[0]=GridDimension[2];dims[1]=GridDimension[1];dims[2]=GridDimension[0];
+
+        for (field = 0; field < NumberOfBaryonFields; field++) {
+            /* This gives back a new reference 
+               So we need to decref it after we add it to the dict */
+            dataset = (PyArrayObject *) PyArray_SimpleNewFromData(
+                    3, dims, ENPY_FLOAT, BaryonField[field]);
+            dataset->flags &= ~NPY_OWNDATA;
+            field_name = PyString_FromString(DataLabel[field]);
+            PyDict_SetItem(grid_data, field_name, (PyObject*) dataset);
+            Py_DECREF(dataset);
+
+            dataset = (PyArrayObject *) PyArray_SimpleNewFromData(
+                    3, dims, ENPY_FLOAT, OldBaryonField[field]);
+            field_name = PyString_FromString(DataLabel[field]);
+            PyDict_SetItem(grid_data, field_name, (PyObject*) dataset);
+            Py_DECREF(dataset);
+            PyDict_SetItem(old_grid_data, field_name, (PyObject*) dataset);
+        }
+
+        grid_id = PyLong_FromLong((long) GridID);
+        PyDict_SetItem(grid_dictionary, grid_id, grid_data);
+        PyDict_SetItem(old_grid_dictionary, grid_id, grid_data);
+        /* New reference from setting, so we decref */
+        Py_DECREF(grid_data);
+        Py_DECREF(old_grid_data);
+
+    }
+    int j = 0;
+    /* Fill our hierarchy information */
+    for (int i = 0; i < 3; i++) {
+        *(int *) PyArray_GETPTR2(container[j], GridID-1, i) =
+            (int) this->GridDimension[i];
+    }
+    j++;
+
+    for (int i = 0; i < 3; i++) {
+        *(int *) PyArray_GETPTR2(container[j], GridID-1, i) =
+            (int) this->GridStartIndex[i];
+    }
+    j++;
+
+    for (int i = 0; i < 3; i++) {
+        *(int *) PyArray_GETPTR2(container[j], GridID-1, i) =
+            (int) this->GridEndIndex[i];
+    }
+    j++;
+
+    for (int i = 0; i < 3; i++) {
+        *(FLOAT *) PyArray_GETPTR2(container[j], GridID-1, i) =
+            (FLOAT) this->GridLeftEdge[i];
+    }
+    j++;
+
+    for (int i = 0; i < 3; i++) {
+        *(FLOAT *) PyArray_GETPTR2(container[j], GridID-1, i) =
+            (FLOAT) this->GridRightEdge[i];
+    }
+    j++;
+
+    *(int *) PyArray_GETPTR2(container[j], GridID-1, 0) =
+        (int) level; j++;
+
+    *(FLOAT *) PyArray_GETPTR2(container[j], GridID-1, 0) =
+        (FLOAT) this->Time; j++;
+
+    *(FLOAT *) PyArray_GETPTR2(container[j], GridID-1, 0) =
+        (FLOAT) this->OldTime; j++;
+
+    *(int *) PyArray_GETPTR2(container[j], GridID-1, 0) =
+        (int) this->ProcessorNumber; j++;
+
+    *(int *) PyArray_GETPTR2(container[j], GridID-1, 0) =
+        (int) this->ReturnNumberOfParticles(); j++;
+
+    *(int *) PyArray_GETPTR2(container[j], GridID-1, 0) =
+        (int) ParentID; j++;
+
+}
+#endif

File src/enzo/InitializePythonInterface.C

+#ifdef EMBEDDED_PYTHON
+/***********************************************************************
+/
+/  INITIALIZE PYTHON INTERFACE AND START INTERPRETER
+/
+/  written by: Matthew Turk
+/  date:       September, 2008
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    SUCCESS or FAIL
+/
+************************************************************************/
+
+#ifdef EMBEDDED_PYTHON
+#define PY_ARRAY_UNIQUE_SYMBOL enzo_ARRAY_API
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#define ENZO_PYTHON_IMPORTED
+#endif
+
+#include <stdlib.h>
+#include <stdio.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 "CosmologyParameters.h"
+#include "TopGridData.h"
+
+int  CosmologyGetUnits(float *DensityUnits, float *LengthUnits,
+		       float *TemperatureUnits, float *TimeUnits,
+		       float *VelocityUnits, FLOAT Time);
+int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
+
+
+static PyMethodDef _EnzoModuleMethods[] = {
+  {NULL, NULL, 0, NULL}
+};
+
+int InitializePythonInterface(int argc, char *argv[])
+{
+#undef int
+
+  //Py_SetProgramName(argv[0]);
+  Py_SetProgramName("embed_enzo");
+
+  Py_Initialize();
+
+  PySys_SetArgv(argc, argv);
+  PyRun_SimpleString("import sys\nsys.path.insert(0,'.')\nsys._parallel = True\n");
+  PyObject *enzo_module, *enzo_module_dict; 
+  enzo_module = Py_InitModule("enzo", _EnzoModuleMethods);
+  enzo_module_dict = PyModule_GetDict(enzo_module);
+  if(enzo_module == NULL){fprintf(stderr, "Failed on Enzo Module!\n");return FAIL;}
+  if(enzo_module_dict == NULL){fprintf(stderr, "Failed on Dict!\n");return FAIL;}
+  PyDict_SetItemString(enzo_module_dict, "grid_data", grid_dictionary);
+  PyDict_SetItemString(enzo_module_dict, "old_grid_data", old_grid_dictionary);
+  PyDict_SetItemString(enzo_module_dict, "hierarchy_information", hierarchy_information);
+  PyDict_SetItemString(enzo_module_dict, "yt_parameter_file", yt_parameter_file);
+  PyDict_SetItemString(enzo_module_dict, "conversion_factors", conversion_factors);
+  import_array1(FAIL);
+  npy_intp flat_dimensions[1];
+  flat_dimensions[0] = (npy_intp) 5;
+  PyObject *temp = PyArray_ZEROS(1, flat_dimensions, NPY_INT, 0);
+  fprintf(stderr, "Completed initialization\n");
+  return SUCCESS;
+}
+
+
+#define TEMP_PYINT(A) Py_XDECREF(temp_int); temp_int = PyLong_FromLong((long) A);
+#define TEMP_PYFLOAT(A) Py_XDECREF(temp_float); temp_float = PyFloat_FromDouble((double) A);
+#define TEMP_PYSTRING(A) Py_XDECREF(temp_string); temp_string = PyString_FromString(A);
+
+int ExportParameterFile(TopGridData &MetaData, FLOAT CurrentTime)
+{
+  /* We need: */
+
+  float DensityUnits = 1, LengthUnits = 1, TemperatureUnits = 1, TimeUnits = 1,
+        VelocityUnits = 1;
+
+  if (ComovingCoordinates)
+    if (CosmologyGetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
+			  &TimeUnits, &VelocityUnits, CurrentTime) == FAIL) {
+      fprintf(stderr, "Error in CosmologyGetUnits.\n");
+      return FAIL;
+    }
+
+  FLOAT a, dadt, FinalRedshift, CurrentRedshift;
+  if (CosmologyComputeExpansionFactor(MetaData.StopTime, &a, &dadt) == FAIL) {
+    fprintf(stderr, "Error in CosmologyComputeExpansionFactor.\n");
+    return FAIL;
+  }
+  FinalRedshift = (1 + InitialRedshift)/a - 1;
+
+  /* Compute the current redshift (for information only). */
+
+  CosmologyComputeExpansionFactor(CurrentTime, &a, &dadt);
+  CurrentRedshift = (1 + InitialRedshift)/a - 1;
+
+
+  PyObject *temp_int = NULL;
+  PyObject *temp_float = NULL;
+  PyObject *temp_string = NULL;
+
+  if (ComovingCoordinates) {
+    TEMP_PYFLOAT(CurrentRedshift);
+    PyDict_SetItemString(yt_parameter_file, "CosmologyCurrentRedshift", temp_float);
+
+    TEMP_PYFLOAT(ComovingBoxSize);
+    PyDict_SetItemString(yt_parameter_file, "CosmologyComovingBoxSize", temp_float);
+
+    TEMP_PYFLOAT(OmegaMatterNow);
+    PyDict_SetItemString(yt_parameter_file, "CosmologyOmegaMatterNow",temp_float);
+
+    TEMP_PYFLOAT(HubbleConstantNow);
+    PyDict_SetItemString(yt_parameter_file, "CosmologyHubbleConstantNow",temp_float);
+
+    TEMP_PYFLOAT(InitialRedshift);
+    PyDict_SetItemString(yt_parameter_file, "CosmologyInitialRedshift",temp_float);
+  }
+
+  TEMP_PYFLOAT(DensityUnits);
+  PyDict_SetItemString(yt_parameter_file, "DensityUnits", temp_float);
+
+  TEMP_PYFLOAT(LengthUnits);
+  PyDict_SetItemString(yt_parameter_file, "LengthUnits", temp_float);
+
+  TEMP_PYFLOAT(TemperatureUnits);
+  PyDict_SetItemString(yt_parameter_file, "TemperatureUnits", temp_float);
+
+  TEMP_PYFLOAT(TimeUnits);
+  PyDict_SetItemString(yt_parameter_file, "TimeUnits", temp_float);
+
+  /*TEMP_PYSTRING(MetaData.MetaDataString);
+  PyDict_SetItemString(yt_parameter_file, "MetaDataString", temp_string);*/
+
+  TEMP_PYINT(HydroMethod);
+  PyDict_SetItemString(yt_parameter_file, "HydroMethod", temp_int);
+
+  TEMP_PYINT(DualEnergyFormalism);
+  PyDict_SetItemString(yt_parameter_file, "DualEnergyFormalism", temp_int);
+
+  TEMP_PYFLOAT(CurrentTime);
+  PyDict_SetItemString(yt_parameter_file, "InitialTime", temp_float);
+
+  TEMP_PYINT(ComovingCoordinates);
+  PyDict_SetItemString(yt_parameter_file, "ComovingCoordinates", temp_int);
+
+  TEMP_PYINT(MultiSpecies);
+  PyDict_SetItemString(yt_parameter_file, "MultiSpecies", temp_int);
+
+  TEMP_PYINT(MetaData.TopGridRank);
+  PyDict_SetItemString(yt_parameter_file, "TopGridRank", temp_int);
+
+  TEMP_PYINT(RefineBy);
+  PyDict_SetItemString(yt_parameter_file, "RefineBy", temp_int);
+
+  TEMP_PYINT(NumberOfPythonCalls);
+  PyDict_SetItemString(yt_parameter_file, "NumberOfPythonCalls", temp_int);
+
+  PyObject *tgd_tuple, *tgd0, *tgd1, *tgd2;
+  
+  /* Construct a tuple */
+  tgd0 = PyLong_FromLong((long) DomainLeftEdge[0]);
+  tgd1 = PyLong_FromLong((long) DomainLeftEdge[1]);
+  tgd2 = PyLong_FromLong((long) DomainLeftEdge[2]);
+  tgd_tuple = PyTuple_Pack(3, tgd0, tgd1, tgd2);
+  PyDict_SetItemString(yt_parameter_file, "DomainLeftEdge", tgd_tuple);
+  Py_XDECREF(tgd_tuple); Py_XDECREF(tgd0); Py_XDECREF(tgd1); Py_XDECREF(tgd2);
+
+  tgd0 = PyLong_FromLong((long) DomainRightEdge[0]);
+  tgd1 = PyLong_FromLong((long) DomainRightEdge[1]);
+  tgd2 = PyLong_FromLong((long) DomainRightEdge[2]);
+  tgd_tuple = PyTuple_Pack(3, tgd0, tgd1, tgd2);
+  PyDict_SetItemString(yt_parameter_file, "DomainRightEdge", tgd_tuple);
+  Py_XDECREF(tgd_tuple); Py_XDECREF(tgd0); Py_XDECREF(tgd1); Py_XDECREF(tgd2);
+
+  tgd0 = PyLong_FromLong((long) MetaData.TopGridDims[0]);
+  tgd1 = PyLong_FromLong((long) MetaData.TopGridDims[1]);
+  tgd2 = PyLong_FromLong((long) MetaData.TopGridDims[2]);
+  tgd_tuple = PyTuple_Pack(3, tgd0, tgd1, tgd2);
+  PyDict_SetItemString(yt_parameter_file, "TopGridDimensions", tgd_tuple);
+  Py_XDECREF(tgd_tuple); Py_XDECREF(tgd0); Py_XDECREF(tgd1); Py_XDECREF(tgd2);
+
+  /* Do the conversion factors */
+  for (int dim = 0; dim < MAX_NUMBER_OF_BARYON_FIELDS; dim++) {
+    if (DataLabel[dim]) {
+      if (strstr(DataLabel[dim], "Density") != NULL) {
+        TEMP_PYFLOAT(DensityUnits);
+	    PyDict_SetItemString(conversion_factors, DataLabel[dim], temp_float);
+        }
+      if (strstr(DataLabel[dim], "velocity") != NULL) {
+        TEMP_PYFLOAT(VelocityUnits);
+	    PyDict_SetItemString(conversion_factors, DataLabel[dim], temp_float);
+        }
+      if (strstr(DataLabel[dim], "_kph") != NULL)  {
+        TEMP_PYFLOAT(1./TimeUnits);
+	    PyDict_SetItemString(conversion_factors, DataLabel[dim], temp_float);
+        }
+    }
+  }
+
+
+  Py_XDECREF(temp_int); Py_XDECREF(temp_float); Py_XDECREF(temp_string);
+  return SUCCESS;
+}
+
+int CallPython(TopGridData &MetaData, FLOAT CurrentTime)
+{
+  ExportParameterFile(MetaData, CurrentTime);
+
+  PyRun_SimpleString("import user_script\nuser_script.main()\n");
+
+  NumberOfPythonCalls++;
+
+  PyDict_Clear(grid_dictionary);
+  PyDict_Clear(hierarchy_information);
+  PyDict_Clear(yt_parameter_file);
+  PyDict_Clear(conversion_factors);
+  return SUCCESS;
+}
+
+#endif

File src/enzo/Make.config.objects

 	EvolveHierarchy.o \
         expand_terms.o \
 	ExternalBoundary_AddField.o \
+	ExposeDataHierarchy.o \
+	ExposeGridHierarchy.o \
 	ExternalBoundary_AppendForcingToBaryonFields.o \
 	ExternalBoundary_constructor.o \
 	ExternalBoundary_DetachForcingFromBaryonFields.o \
 	Grid_constructor.o \
 	Grid_ConvertTotalEnergyToGasEnergy.o \
 	Grid_CoolingTestInitializeGrid.o \
+	Grid_ConvertToNumpy.o \
 	Grid_CopyBaryonFieldToOldBaryonField.o \
 	Grid_CopyOverlappingMassField.o \
 	Grid_CopyParentToGravitatingFieldBoundary.o \
 	InitializeEquilibriumCoolData.o \
 	InitializeLocal.o \
 	InitializeNew.o \
+	InitializePythonInterface.o \
 	InitializeRadiationFieldData.o \
 	InitializeRateData.o \
 	interp1d.o \

File src/enzo/SetDefaultGlobalValues.C

   MetalCooling = FALSE;
   MetalCoolingTable = (char*) "metal_cool.dat";
 
+#ifdef EMBEDDED_PYTHON
+  NumberOfPythonCalls = 0;
+  grid_dictionary = PyDict_New();
+  old_grid_dictionary = PyDict_New();
+  hierarchy_information = PyDict_New();
+  yt_parameter_file = PyDict_New();
+  conversion_factors = PyDict_New();
+#endif
+
   return SUCCESS;
 }

File src/enzo/enzo.C

 /    doing a new run, a restart, an extraction, or a projection.
 /
 ************************************************************************/
+
+ 
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
  
 #ifdef USE_MPI
 #include "mpi.h"
 #endif /* USE_MPI */
  
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <unistd.h>
- 
 #include "svn_version.def"
 #include "performance.h"
 #include "macros_and_parameters.h"
 #include "PhotonCommunication.h"
 #endif
 #undef DEFINE_STORAGE
+#ifdef EMBEDDED_PYTHON
+int InitializePythonInterface(int argc, char **argv);
+#endif
  
 // Function prototypes
  
   fprintf(memtracePtr, "Call evolve hierarchy %8"ISYM"  %16"ISYM" \n", MetaData.CycleNumber, MemInUse);
 #endif
 
+#ifdef EMBEDDED_PYTHON
+  // We initialize our Python interface now
+  if(debug)fprintf(stdout, "Initializing Python interface\n");
+  InitializePythonInterface(argc, argv);
+#endif 
 
 
  

File src/enzo/global_data.h

 EXTERN char GlobalPath[MAX_LINE_LENGTH];
 #endif
 
+#ifdef EMBEDDED_PYTHON
+EXTERN int NumberOfPythonCalls;
+EXTERN PyObject *grid_dictionary;
+EXTERN PyObject *old_grid_dictionary;
+EXTERN PyObject *hierarchy_information;
+EXTERN PyObject *yt_parameter_file;
+EXTERN PyObject *conversion_factors;
+#endif
 /* Multi-species rate equation flag and associated data. */
 
 EXTERN int MetalCooling;

File src/enzo/macros_and_parameters.h

 / MACRO DEFINITIONS AND PARAMETERS
 /
 ************************************************************************/
+#ifdef EMBEDDED_PYTHON
+#ifndef ENZO_PYTHON_IMPORTED
+#define PY_ARRAY_UNIQUE_SYMBOL enzo_ARRAY_API
+#define NO_IMPORT_ARRAY 
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+#endif
+
 
 #include "message.h"
 
 #define HDF5_REAL HDF5_R4
 #define HDF5_FILE_REAL HDF5_FILE_R8
 #endif
+#ifdef EMBEDDED_PYTHON
+#define ENPY_FLOAT NPY_FLOAT
+#endif
 #endif
 
 #ifdef CONFIG_BFLOAT_8
 #define HDF5_REAL HDF5_R8
 #define HDF5_FILE_REAL HDF5_FILE_R8
 #endif
+#ifdef EMBEDDED_PYTHON
+#define ENPY_FLOAT NPY_DOUBLE
+#endif
 #endif
 
 #ifdef CONFIG_PFLOAT_4
 #define FLOATDataType MPI_DOUBLE
 #define HDF5_PREC HDF5_R8
 #define HDF5_FILE_PREC HDF5_FILE_R8
+#ifdef EMBEDDED_PYTHON
+#define ENPY_FLOAT NPY_DOUBLE
+#endif
 #endif
 
 #ifdef CONFIG_PFLOAT_16
 #define FLOATDataType MPI_LONG_DOUBLE
 #define HDF5_PREC HDF5_R16
 #define HDF5_FILE_PREC HDF5_FILE_R16
+#ifdef EMBEDDED_PYTHON
+#define ENPY_FLOAT NPY_LONGDOUBLE
+#endif
 #endif