Allocate default_material in the MPI shared_memory when copy Numpy arrays to meep

configure.ac

@@ -398,7 +398,7 @@ CPPFLAGS="$CPPFLAGS $GUILE_CPPFLAGS"
 
 AC_MSG_CHECKING([for libctl dir])
 if test x != x"$LIBCTL_DIR" -a -r "$LIBCTL_DIR/share/libctl/base/ctl.scm"; then
-	LIBCTL_DIR="$LIBCTL_DIR/share/libctl"
+	LIBCTL_DIR="$LIBCTL_DIR/share/libctl"
 fi
 if test x != x"$LIBCTL_DIR" -a ! -r "$LIBCTL_DIR/base/ctl.scm"; then
 	LIBCTL_DIR=""

python/Makefile.am

@@ -111,7 +111,7 @@ TESTS_ENVIRONMENT = export PYTHONPATH=$(abs_top_builddir)/python:$$PYTHONPATH;
 # the low-level interface between python codes and the libmeep C++ library
 ######################################################################
 BUILT_SOURCES = meep-python.cxx __init__.py meep.py
-EXTRA_DIST = $(BUILT_SOURCES) meep.i meep-python.hpp typemap_utils.cpp materials.py examples tests
+EXTRA_DIST = $(BUILT_SOURCES) meep.i meep-python.hpp typemap_utils.cpp shared_memory.hpp materials.py examples tests
 CLEANFILES = .coverage
 MAINTAINERCLEANFILES = $(BUILT_SOURCES)
 

python/shared_memory.hpp

@@ -0,0 +1,114 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+
+
+#ifdef HAVE_MPI
+
+  #include <mpi.h>
+  void MPI_configuration(MPI_Comm *allcomm, MPI_Comm *nodecomm, int*node_size, int*node_rank, int*world_size, int*world_rank)
+      {
+        // Read MPI properties
+        *allcomm = MPI_COMM_WORLD;
+        int MPI_initialized_already;
+
+        MPI_Initialized(&MPI_initialized_already);
+        if (!MPI_initialized_already)
+            MPI_Init(NULL, NULL);
+
+
+        MPI_Comm_size(*allcomm, world_size);
+        MPI_Comm_rank(*allcomm, world_rank);
+
+        // Create node-local communicator
+        MPI_Comm_split_type(*allcomm, MPI_COMM_TYPE_SHARED, *world_rank, MPI_INFO_NULL, nodecomm);
+        MPI_Comm_size(*nodecomm, node_size);
+        MPI_Comm_rank(*nodecomm, node_rank);
+    }
+
+  template<class T>
+  T* SharedMemoryArray(size_t arraysize) {
+      // Array properties
+      T *array, *localarray;
+      int localarraysize;
+
+      // MPI configuration:
+      MPI_Comm allcomm, nodecomm;
+      int node_size, node_rank;
+      int world_size, world_rank;
+      MPI_configuration(&allcomm, &nodecomm, &node_size, &node_rank, &world_size, &world_rank);
+
+      // Shared memmory vars
+      int *model; 
+      int i, flag;
+      int windisp;
+      MPI_Aint winsize;
+      MPI_Win winarray;
+
+      // Only rank 0 on a node actually allocates memory
+      localarraysize = 0;
+      if (node_rank == 0) localarraysize = arraysize;
+
+      MPI_Win_allocate_shared(localarraysize*sizeof(T), sizeof(T),
+              MPI_INFO_NULL, nodecomm, &localarray, &winarray);
+
+      MPI_Win_get_attr(winarray, MPI_WIN_MODEL, &model, &flag);
+
+      // need to get local pointer valid for array on rank 0
+
+      array = localarray;
+
+      if (node_rank != 0)
+      {
+          MPI_Win_shared_query(winarray, 0, &winsize, &windisp, &array);
+      }
+
+      printf("POINTER TO SHARED ARRAY IN PROCESS %d:   = %p,  SHARED= %p\n", world_rank , (void *) &array, (void *) array);
+
+
+      MPI_Win_fence(0, winarray);
+
+      return array;
+
+    }
+
+  template<class T>
+  void memcpy_shared(T* dest, const T* src, size_t count ) 
+  {
+    // MPI configuration:
+    MPI_Comm allcomm, nodecomm;
+    int node_size, node_rank;
+    int world_size, world_rank;
+    MPI_configuration(&allcomm, &nodecomm, &node_size, &node_rank, &world_size, &world_rank);
+
+    if (node_rank == 0) {
+        memcpy(dest, src, count);
+        std::cout << "copying the values in the shared memory ["<<node_rank << "/" <<node_size<<"] \n" ;
+    }
+    MPI_Barrier(MPI_COMM_WORLD);
+  }
+
+#else
+
+  template<class T>
+  T* SharedMemoryArray(size_t arraysize) {
+    return new T[arraysize];
+    }
+
+  template<class T>
+  void memcpy_shared(T* dest, const T* src, size_t count ) {
+    std::cout << "copying the values on every process \n" ;
+    memcpy(dest, src, count);
+    }
+
+#endif
+
+
+
+
+
+
+
+
+
+

python/typemap_utils.cpp

@@ -35,6 +35,8 @@
 #define SWIG_PYTHON_THREAD_SCOPED_BLOCK SWIG_PYTHON_THREAD_BEGIN_BLOCK
 #endif
 
+#include "shared_memory.hpp"
+
 PyObject *py_callback = NULL;
 PyObject *py_callback_v3 = NULL;
 PyObject *py_amp_func = NULL;
@@ -522,8 +524,10 @@ static int pymaterial_grid_to_material_grid(PyObject *po, material_data *md) {
   PyArrayObject *pao = (PyArrayObject *)po_dp;
   if (!PyArray_Check(pao)) { meep::abort("MaterialGrid weights failed to init."); }
   if (!PyArray_ISCARRAY(pao)) { meep::abort("Numpy array weights must be C-style contiguous."); }
-  md->weights = new double[PyArray_SIZE(pao)];
-  memcpy(md->weights, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
+  md->weights = SharedMemoryArray<double>(PyArray_SIZE(pao));
+  memcpy_shared(md->weights, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
+  // md->weights = new double[PyArray_SIZE(pao)];
+  // memcpy(md->weights, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
 
   // if needed, combine sus structs to main object
   PyObject *py_e_sus_m1 = PyObject_GetAttrString(po_medium1, "E_susceptibilities");
@@ -608,8 +612,10 @@ static int pymaterial_to_material(PyObject *po, material_type *mt) {
     md = new material_data();
     md->which_subclass = material_data::MATERIAL_FILE;
     md->epsilon_dims[0] = md->epsilon_dims[1] = md->epsilon_dims[2] = 1;
-    md->epsilon_data = new double[PyArray_SIZE(pao)];
-    memcpy(md->epsilon_data, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
+    md->epsilon_data = SharedMemoryArray<double>(PyArray_SIZE(pao));
+    memcpy_shared(md->epsilon_data, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
+    //md->epsilon_data = new double[PyArray_SIZE(pao)];
+    //memcpy(md->epsilon_data, (double *)PyArray_DATA(pao), PyArray_SIZE(pao) * sizeof(double));
 
     for (int i = 0; i < PyArray_NDIM(pao); ++i) {
       md->epsilon_dims[i] = (size_t)PyArray_DIMS(pao)[i];

src/meepgeom.cpp

@@ -107,6 +107,7 @@ void material_gc(material_type m) {
   m->medium_2.H_susceptibilities.clear();
 }
 
+
 void material_free(material_type m) {
   if (!m) return;
 
@@ -119,11 +120,17 @@ void material_free(material_type m) {
 
   // NOTE: We do not delete the user_data field here since it is an opaque/void
   // object so will assume that the caller keeps track of its lifetime.
-  delete[] m->epsilon_data;
-  m->epsilon_data = NULL;
+
+  //MPI_Win_free(&md->epsilon_adress);
+  //MPI_Win_free(&md->weights_adress);
+
+  // delete[] m->epsilon_data;
+  //m->epsilon_data = NULL;
+  //delete[] m->epsilon_adress;
 
-  delete[] m->weights;
-  m->weights = NULL;
+  // delete[] m->weights;
+  //m->weights = NULL;
+  //delete[] m->weights_adress;
   delete m;
 }