update mpi api, debug needed

Makefile.Darwin

@@ -2,13 +2,17 @@ OBJDIR = build
 SHLIB_EXT = so
 SHLIB_FLAGS = -bundle -undefined dynamic_lookup
 
-CFLAGS += -Wall -O3 -ffast-math -march=native -fopenmp
-CXXFLAGS += -std=c++11 -Wall -O3 -ffast-math -march=native -fopenmp
+CFLAGS += -Wall -O3 -ffast-math -march=native
+CXXFLAGS += -std=c++11 -Wall -O3 -ffast-math -march=native
 LA_LIBS = -O3 -ffast-math -framework Accelerate
 
 LUA_INC = -I/usr/local/include
-LUA_LIB = -L/usr/local/bin -llua
+LUA_LIB = -L/usr/local/bin -llua -lm
 
 #### change to proper compilers
-CC = gcc
-CXX = g++
+CC = cc
+CXX = c++
+
+#### change to proper MPI compilers
+MPICC = mpicc
+MPICXX = mpicxx

Makefile.Linux

@@ -8,9 +8,13 @@ CXXFLAGS += -std=c++11 -O3 -Wall -march=native -fcx-limited-range -fPIC -fopenmp
 LA_LIBS = -lopenblas -llapack
 
 LUA_INC = -I/usr/local/include
-LUA_LIB = -llua -ldl
+LUA_LIB = -llua -ldl -lm
 
 
 #### change to proper compilers
 CC = gcc
-CXX = g++
+CXX = g++
+
+#### change to proper MPI compilers
+MPICC = mpicc
+MPICXX = mpic++

Makefile.Msys2

@@ -13,4 +13,8 @@ LUA_MODULE_LIB = -DLUA_BUILD_AS_DLL -Wl,-Bdynamic -llua.dll -Wl,-Bstatic
 
 #### change to proper compilers
 CC = gcc
-CXX = g++
+CXX = g++
+
+#### change to proper MPI compilers
+MPICC = mpicc
+MPICXX = mpic++

Makefile.common

@@ -49,14 +49,14 @@ $(OBJDIR)/mesh: $(OBJDIR)/MESH_o/main_lua.o $(OBJDIR)/$(MESH_LIBNAME)
 	$(CXX) $(CXXFLAGS) $(CPPFLAGS) $(ARMAFLAG) $< -o $@ -L$(OBJDIR) -lmesh $(LA_LIBS) $(LUA_LIB)
 
 ### for MPI only
-$(OBJDIR)/MESH_o/buffer.o: src/luawrapper/buffer.c objdir
-	$(CC) -c -DHAVE_MPI $(CFLAGS) $(CPPFLAGS) $(LUA_INC) $(ARMAFLAG) $< -o $@ $(LA_LIBS)
-$(OBJDIR)/MESH_o/lua_mpi.o: src/luawrapper/lua_mpi.c src/luawrapper/mpifuncs.c objdir
-	$(CC) -c -DHAVE_MPI $(CFLAGS) $(CPPFLAGS) $(LUA_INC) $(ARMAFLAG) $< -o $@ $(LA_LIBS)
+$(OBJDIR)/MESH_o/buffer_mpi.o: src/luawrapper/buffer_mpi.c objdir
+	$(MPICC) -c -DHAVE_MPI $(CFLAGS) $(CPPFLAGS) $(LUA_INC) $(ARMAFLAG) $< -o $@ $(LA_LIBS)
+$(OBJDIR)/MESH_o/lua_mpi.o: src/luawrapper/lua_mpi.c src/luawrapper/funcs_mpi.c objdir
+	$(MPICC) -c -DHAVE_MPI $(CFLAGS) $(CPPFLAGS) $(LUA_INC) $(ARMAFLAG) $< -o $@ $(LA_LIBS)
 $(OBJDIR)/MESH_o/main_lua_mpi.o: src/main_lua.cpp objdir
-	$(CXX) -c -DHAVE_MPI $(CXXFLAGS) $(CPPFLAGS) $(LUA_INC) $(ARMAFLAG) $< -o $@ $(LA_LIBS)
-$(OBJDIR)/meshMPI: $(OBJDIR)/MESH_o/main_lua_mpi.o $(OBJDIR)/MESH_o/lua_mpi.o $(OBJDIR)/MESH_o/buffer.o  $(OBJDIR)/$(MESH_LIBNAME)
-	$(CXX) $(CXXFLAGS) $(CPPFLAGS) $(ARMAFLAG) $< -o $@ -L$(OBJDIR) -lmesh $(LA_LIBS) $(LUA_LIB)
+	$(MPICXX) -c -DHAVE_MPI $(CXXFLAGS) $(CPPFLAGS) $(LUA_INC) $(ARMAFLAG) $< -o $@ $(LA_LIBS)
+$(OBJDIR)/meshMPI: $(OBJDIR)/MESH_o/main_lua_mpi.o $(OBJDIR)/MESH_o/lua_mpi.o $(OBJDIR)/MESH_o/buffer_mpi.o  $(OBJDIR)/$(MESH_LIBNAME)
+	$(MPICXX) $(CXXFLAGS) $(CPPFLAGS) $(ARMAFLAG) $< -o $@ -L$(OBJDIR) -lmesh $(LA_LIBS) $(LUA_LIB)
 
 clean:
 	rm -rf $(OBJDIR)

examples/Arxiv_1701_02986/main.lua

@@ -0,0 +1,55 @@
+-- this is a lua realization of the corresponding cpp file
+constants = Constants();
+
+function thetaDerivative(omega, T)
+  local theta = constants.h_bar * omega / (math.exp(constants.h_bar * omega/constants.k_B/T) - 1);
+  return math.pow(theta, 2) * math.exp(constants.h_bar * omega/constants.k_B/T) /constants.k_B / math.pow(T, 2);
+end
+
+
+s1 = SimulationPlanar.new();
+s1:AddMaterial("Si", "Si.txt");
+s1:AddMaterial("Vacuum", "Vacuum.txt");
+s1:AddLayer("SiBottom", 0, "Si");
+s1:AddLayer("VacGap", 20e-9, "Vacuum");
+s1:AddLayerCopy("SiTop", "SiBottom");
+
+s1:SetSourceLayer("SiBottom");
+s1:SetProbeLayer("VacGap");
+
+s1:SetKParallelIntegral(500);
+s1:OptUseQuadgk();
+s1:BuildRCWA();
+s1:IntegrateKParallel();
+phi = s1:GetPhi();
+omega = s1:GetOmega();
+for i = 1, s1:GetNumOfOmega() do
+  local E = constants.h_bar * omega[i] / constants.q;
+  local spectrum = thetaDerivative(omega[i], 300) * phi[i] * constants.q / constants.h_bar;
+  print(spectrum/1e5);
+end
+
+f = 0.98;
+width = math.sqrt(f * 50e-9 * 50e-9);
+s2 = SimulationPattern.new();
+s2:SetPeriodicity(50e-9, 50e-9);
+s2:SetGx(10);
+s2:SetGy(10);
+
+s2:AddMaterial("Si", "Si.txt");
+s2:AddMaterial("Vacuum", "Vacuum.txt");
+s2:AddLayer("SiBottom", 0, "Si");
+s2:SetLayerPatternRectangle("SiBottom", "Vacuum", {25e-9, 25e-9}, {width, width});
+s2:AddLayer("VacGap", 20e-9, "Vacuum");
+s2:AddLayerCopy("SiTop", "SiBottom");
+
+s2:SetSourceLayer("SiBottom");
+s2:SetProbeLayer("VacGap");
+s2:GetSysInfo();
+
+s2:OptPrintIntermediate();
+s2:SetKxIntegralSym(20, 60);
+s2:SetKyIntegralSym(20, 60);
+s2:BuildRCWA();
+s2:IntegrateKxKy();
+

examples/MPI_example/main.lua

@@ -0,0 +1,66 @@
+s = SimulationGrating.new();
+s:SetPeriodicity(1e-6);
+s:AddMaterial("Au", "fullGold.txt");
+s:AddMaterial("Vacuum", "fullVacuum.txt");
+
+s:AddLayer("BottomAir", 0, "Vacuum");
+s:AddLayer("GoldSubstrateBottom", 0.5e-6, "Au");
+s:AddLayer("GoldGratingBottom", 5e-6,"Au");
+s:SetLayerPatternGrating("GoldGratingBottom", "Vacuum", 0.5e-6, 0.2e-6);
+
+s:AddLayer("VacGap", 1e-6, "Vacuum");
+s:AddLayerCopy("GoldGratingTop", "GoldGratingBottom");
+s:AddLayerCopy("GoldSubstrateTop", "GoldSubstrateBottom");
+s:AddLayerCopy("TopAir", "BottomAir");
+
+s:SetSourceLayer("GoldSubstrateBottom");
+s:SetSourceLayer("GoldGratingBottom");
+s:SetProbeLayer("VacGap");
+
+s:OptPrintIntermediate();
+s:SetThread(1);
+s:SetGx(50);
+s:SetKxIntegralSym(500);
+s:SetKyIntegralSym(200, 5);
+s:BuildRCWA();
+
+-- start MPI
+local sizeb = buffer.new_buffer(buffer.sizeof(buffer.int))
+local rankb = buffer.new_buffer(buffer.sizeof(buffer.int))
+
+MPI.Init()
+MPI.Comm_rank(MPI.COMM_WORLD, rankb)
+MPI.Comm_size(MPI.COMM_WORLD, sizeb)
+
+local size = buffer.get_typed(sizeb, buffer.int, 0)
+local rank = buffer.get_typed(rankb, buffer.int, 0)
+
+if rank == 0 then
+  s:GetSysInfo();
+end
+
+status = MPI.Status()
+numOfOmega = s:GetNumOfOmega()
+
+if rank == 0 then
+  print(rank);
+  s:IntegrateKxKyMPI(rank, size);
+  phi_master = s:GetPhi();
+  for i = 1,size - 1 do
+    local phi_local = buffer.new_buffer(buffer.sizeof(numOfOmega * buffer.double));
+    MPI.Recv(phi_local, numOfOmega, MPI.DOUBLE, 0, MPI.COMM_WORLD, status);
+    for j = 1, numOfOmega do
+      phi_master[j] = phi_master[j] + phi_local[j];
+    end
+  end
+  for i = 1,numOfOmega do
+    print(string.format("%e", phi_master[i]));
+  end
+else
+  print(rank);
+  s:IntegrateKxKyMPI(rank, size);
+  local phi_slave = s:GetPhi();
+  MPI.Send(phi_slave, numOfOmega, MPI.DOUBLE, 0, 0, MPI.COMM_WORLD);
+end
+
+MPI.Finalize();

examples/PRApplied_6_024014_2016/diagonal.cpp

@@ -25,6 +25,5 @@ int main(){
   s->setThread(4);
   s->buildRCWA();
   s->integrateKParallel();
-  s->outputPhi("output.txt");
   return 0;
 }

examples/PRApplied_6_024014_2016/diagonal.lua

@@ -19,5 +19,3 @@ s:SetThread(4);
 s:OptUseQuadgk();
 s:BuildRCWA();
 s:IntegrateKParallel();
-
-s:OutputPhi("output_lua.txt");

examples/Simple_Planar_Far/test.cpp

@@ -21,6 +21,10 @@ int main(){
   s->setThread(4);
   s->buildRCWA();
   s->integrateKParallel();
-  s->outputPhi("test_output.txt");
+
+  double* phi = s->getPhi();
+  for(int i = 0; i < s->getNumOfOmega(); i++){
+    std::cout << phi[i] << std::endl;
+  }
   return 0;
 }

examples/Simple_Planar_Far/test.lua

@@ -15,4 +15,7 @@ s:OptUseQuadgk();
 s:BuildRCWA();
 s:SetKParallelIntegral(1);
 s:IntegrateKParallel();
-s:OutputPhi("test_lua.txt");
+phi = s:GetPhi();
+for i = 1, s:GetNumOfOmega() do
+  print(phi[i])
+end

examples/iterate_over_gap/iterate.cpp

@@ -26,7 +26,6 @@ int main(){
     s->setLayerThickness("VacGap", i * 1e-9);
     s->buildRCWA();
     s->integrateKParallel();
-    s->outputPhi("test" + std::to_string(i) + ".txt");
   }
   return 0;
 }

examples/iterate_over_gap/iterate.lua

@@ -21,7 +21,6 @@ for i = 10, 100, 10 do
   s:SetLayerThickness("VacGap", i * 1e-9);
   s:BuildRCWA();
   s:IntegrateKParallel();
-  s:OutputPhi("test_lua"..tostring(i)..".txt");
 end
 
 

src/Mesh.cpp

@@ -767,6 +767,9 @@ namespace MESH{
       throw UTILITY::ValueException("Set integration range first!");
     }
     resultArray_ = new double[numOfKx_ * numOfKy_ * numOfOmega_];
+    for(int i = 0; i < numOfKx_ * numOfKy_ * numOfOmega_; i++){
+      resultArray_[i] = 0;
+    }
 
     EMatricesVec_.resize(numOfOmega_);
     grandImaginaryMatrixVec_.resize(numOfOmega_);
@@ -1215,7 +1218,6 @@ namespace MESH{
 
     if(parallel){
       #if defined(_OPENMP)
-      std::cout << "here" << std::endl;
         #pragma omp parallel for num_threads(numOfThread_)
       #endif
       for(int i = start; i < end; i++){

src/luawrapper/lua_mpi.c

@@ -74,7 +74,7 @@ static int _MPI_Init(lua_State *L)
   lua_pushnumber(L, ret);
   return 1;
 }
-#include "mpifuncs.c"
+#include "funcs_mpi.c"
 
 
 static void register_constants(lua_State *L)

src/main_lua.cpp

@@ -29,7 +29,7 @@ extern "C"
 #ifdef HAVE_MPI
 #include <mpi.h>
 #include "luawrapper/lua_mpi.c"
-#include "luawrapper/buffer.c"
+#include "luawrapper/buffer_mpi.c"
 #endif
 
 using namespace MESH;
@@ -724,8 +724,8 @@ lua_State* new_MESH_lua_State(){
 
   luaL_openlibs(L);
   #ifdef HAVE_MPI
-    luaL_requiref(L, "MPI", luaopen_MPI, 0); lua_pop(L, 1);
-    luaL_requiref(L, "buffer", luaopen_buffer, 0); lua_pop(L, 1);
+    luaL_requiref(L, "MPI", &luaopen_MPI, 1); lua_pop(L, 1);
+    luaL_requiref(L, "buffer", &luaopen_buffer, 1); lua_pop(L, 1);
   #endif
   return L;
 }