Merge pull request #45 from natj/phys2

Physics update

CMakeLists.txt

@@ -12,10 +12,12 @@ set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
 # optimization flags
 #set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -Wall")
 #set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -ffast-math -fprofile-generate -fprofile-use")
-#set(CMAKE_CXX_FLAGS_RELEASE "-O2 -march=native -ftree-vectorize -fopt-info-vec")
-set(CMAKE_CXX_FLAGS_RELEASE "-O2 -march=haswell -mtune=haswell") #latest optimization
+set(CMAKE_CXX_FLAGS_RELEASE "-O2 -march=haswell -ftree-vectorize -fopt-info-all-all=all.all")
 #set(CMAKE_CXX_FLAGS_RELEASE "-O2 -march=native -mtune=native") #general optimization
 
+# epyc rome
+#set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=znver1 -mtune=znver1 -mfma -mavx2 -m3dnow -fomit-frame-pointer -fopt-info-all-all=all.all")
+
 
 # AddressSanitizer
 set(CMAKE_CXX_FLAGS_ASAN "-fsanitize=address -fno-optimize-sibling-calls -fsanitize-address-use-after-scope -fno-omit-frame-pointer -g -O1")
@@ -26,6 +28,8 @@ set(CMAKE_CXX_FLAGS_LSAN "-fsanitize=leak -fno-omit-frame-pointer -g -O1")
 # MemorySanitizer
 set(CMAKE_CXX_FLAGS_MSAN "-fsanitize=memory -fno-optimize-sibling-calls -fsanitize-memory-track-origins=2 -fno-omit-frame-pointer -g -O2")
 
+message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
+message(STATUS "Flags: ${CMAKE_CXX_FLAGS_RELEASE}")
 
 
 # set location of python libraries

archs/rusty.sh

@@ -1,13 +1,13 @@
 #modules
-module load slurm/18.08.8
+module load slurm
 module load gcc/7.4.0
 module load openmpi2/2.1.6-hfi
 module load python3/3.7.3
 module load python3-mpi4py/3.7.3-openmpi2
 module load cmake/3.14.5
-module load modules-nix/20200227-37
-module load nix/ffmpeg/3.4.7
-module load nix/fftw/3.3.8
+module load modules-nix
+module load nix/ffmpeg
+module load nix/fftw
 module load lib/hdf5/1.8.21-openmpi2
 
 

bindings/CMakeLists.txt

@@ -12,10 +12,13 @@ set (TOOLS_FILES
 set (FIELDS_FILES 
      pyfields.c++
      ../em-fields/tile.c++ 
-     ../em-fields/damping_tile.c++
      ../em-fields/propagator/fdtd2.c++ 
+     ../em-fields/propagator/fdtd2_pml.c++ 
      ../em-fields/propagator/fdtd4.c++ 
+     ../em-fields/propagator/fdtd_general.c++ 
      ../em-fields/filters/digital.c++
+     ../em-fields/boundaries/damping_tile.c++
+     ../em-fields/boundaries/conductor.c++
      )
 
 set (VLV_FILES 
@@ -39,6 +42,8 @@ set (PIC_FILES
      ../pic/pushers/boris_rad.c++
      ../pic/pushers/boris_grav.c++
      ../pic/pushers/vay.c++
+     ../pic/pushers/higuera_cary.c++
+     ../pic/pushers/rgca.c++
      ../pic/interpolators/linear.c++
      ../pic/depositers/zigzag.c++
      )

bindings/pyfields.c++

@@ -8,15 +8,19 @@
 #include "../tools/mesh.h"
 
 #include "../em-fields/tile.h"
-#include "../em-fields/damping_tile.h"
 
 #include "../em-fields/propagator/propagator.h"
 #include "../em-fields/propagator/fdtd2.h"
+#include "../em-fields/propagator/fdtd2_pml.h"
 #include "../em-fields/propagator/fdtd4.h"
+#include "../em-fields/propagator/fdtd_general.h"
 
 #include "../em-fields/filters/filter.h"
 #include "../em-fields/filters/digital.h"
 
+#include "../em-fields/boundaries/damping_tile.h"
+#include "../em-fields/boundaries/conductor.h"
+
 #include "../io/writers/writer.h"
 #include "../io/writers/fields.h"
 #include "../io/snapshots/fields.h"
@@ -37,6 +41,8 @@ auto declare_tile(
     py::module& m,
     const std::string& pyclass_name) 
 {
+  std::vector<int> iarr{0,1,2};
+
   return py::class_<
              fields::Tile<D>,
               corgi::Tile<D>, 
@@ -47,9 +53,10 @@ auto declare_tile(
     .def("cycle_yee",           &fields::Tile<D>::cycle_yee)
     .def("clear_current",       &fields::Tile<D>::clear_current)
     .def("deposit_current",     &fields::Tile<D>::deposit_current)
-    .def("update_boundaries",   &fields::Tile<D>::update_boundaries)
     .def("exchange_currents",   &fields::Tile<D>::exchange_currents)
-
+    .def("update_boundaries",   &fields::Tile<D>::update_boundaries,
+            py::arg("grid"),
+            py::arg("iarr")=iarr)
     .def("get_yee",             &fields::Tile<D>::get_yee,
         py::arg("i")=0,
         py::return_value_policy::reference,
@@ -154,6 +161,20 @@ class PyFDTD4 : public FDTD4<D>
 };
 
 
+template<size_t D>
+class PyFDTDGen : public FDTDGen<D>
+{
+  using FDTDGen<D>::FDTDGen;
+
+  void push_e( Tile<D>& tile ) override {
+  PYBIND11_OVERLOAD_PURE( void, FDTDGen<D>, push_e, tile);
+  }
+
+  void push_half_b( Tile<D>& tile ) override {
+  PYBIND11_OVERLOAD_PURE( void, FDTDGen<D>, push_half_b, tile);
+  }
+};
+
 // templated base-class for Propagator; 
 // see https://github.com/pybind/pybind11/blob/master/tests/test_virtual_functions.cpp
 //template <class Base = Propagator<1>, size_t D=1>
@@ -383,13 +404,39 @@ void bind_fields(py::module& m_sub)
     .def(py::init<>())
     .def_readwrite("corr",     &fields::FDTD2<3>::corr);
 
+  // fdtd2 propagator with perfectly matched ouer layer
+  py::class_<fields::FDTD2_pml<3>>(m_3d, "FDTD2_pml", fieldspropag3d)
+    .def(py::init<>())
+    .def_readwrite("cenx",     &fields::FDTD2_pml<3>::cenx)
+    .def_readwrite("ceny",     &fields::FDTD2_pml<3>::ceny)
+    .def_readwrite("cenz",     &fields::FDTD2_pml<3>::cenz)
+    .def_readwrite("radx",     &fields::FDTD2_pml<3>::radx)
+    .def_readwrite("rady",     &fields::FDTD2_pml<3>::rady)
+    .def_readwrite("radz",     &fields::FDTD2_pml<3>::radz)
+    .def_readwrite("rad_lim",  &fields::FDTD2_pml<3>::rad_lim)
+    .def_readwrite("norm_abs", &fields::FDTD2_pml<3>::norm_abs)
+    .def_readwrite("corr",     &fields::FDTD2_pml<3>::corr)
+    .def("push_eb",            &fields::FDTD2_pml<3>::push_eb);
+
 
   // fdtd4 propagator
   py::class_<fields::FDTD4<3>, Propagator<3>, PyFDTD4<3> >(m_3d, "FDTD4")
     .def_readwrite("corr",     &fields::FDTD4<3>::corr)
     .def(py::init<>());
 
 
+  // fdtd general propagator
+  //py::class_<fields::FDTDGen<3>, Propagator<3>, PyFDTDGen<3> >(m_3d, "FDTDGen")
+  py::class_<fields::FDTDGen<3>>(m_3d, "FDTDGen")
+    .def_readwrite("corr",     &fields::FDTDGen<3>::corr)
+    .def_readwrite("CXs",      &fields::FDTDGen<3>::CXs, py::return_value_policy::reference,py::keep_alive<1,0>())
+    .def_readwrite("CYs",      &fields::FDTDGen<3>::CYs, py::return_value_policy::reference,py::keep_alive<1,0>())
+    .def_readwrite("CZs",      &fields::FDTDGen<3>::CZs, py::return_value_policy::reference,py::keep_alive<1,0>())
+    .def("push_e",             &fields::FDTDGen<3>::push_e)
+    .def("push_half_b",        &fields::FDTDGen<3>::push_half_b)
+    .def(py::init<>());
+
+
   //--------------------------------------------------
   // 2D Filter bindings
   py::class_< fields::Filter<2>, PyFilter<2> > fieldsfilter2d(m_2d, "Filter");
@@ -437,6 +484,27 @@ void bind_fields(py::module& m_sub)
     .def("solve",      &fields::Binomial2<3>::solve);
 
 
+  //--------------------------------------------------
+  // EM boundary conditions
+
+  // 3D rotating conductor
+  py::class_<fields::Conductor<3>>(m_3d, "Conductor")
+    .def(py::init<>())
+    .def_readwrite("B0",    &fields::Conductor<3>::B0)
+    .def_readwrite("radius",&fields::Conductor<3>::radius)
+    .def_readwrite("period",&fields::Conductor<3>::period)
+    .def_readwrite("chi",   &fields::Conductor<3>::chi)
+    .def_readwrite("phase", &fields::Conductor<3>::phase)
+    .def_readwrite("cenx",  &fields::Conductor<3>::cenx)
+    .def_readwrite("ceny",  &fields::Conductor<3>::ceny)
+    .def_readwrite("cenz",  &fields::Conductor<3>::cenz)
+    .def_readwrite("delta", &fields::Conductor<3>::delta)
+    .def("insert_em",       &fields::Conductor<3>::insert_em)
+    .def("update_e",        &fields::Conductor<3>::update_e)
+    .def("update_b",        &fields::Conductor<3>::update_b);
+
+
+
   //--------------------------------------------------
   // Snapshot IO 
 

bindings/pypic.c++

@@ -13,6 +13,8 @@ namespace py = pybind11;
 #include "../pic/pushers/boris_rad.h"
 #include "../pic/pushers/boris_grav.h"
 #include "../pic/pushers/vay.h"
+#include "../pic/pushers/higuera_cary.h"
+#include "../pic/pushers/rgca.h"
 
 #include "../pic/interpolators/interpolator.h"
 #include "../pic/interpolators/linear.h"
@@ -144,7 +146,58 @@ auto declare_prtcl_container(
           int nparts = s.size();
           assert(i < nparts);
           return s.Bpart[i + 2*nparts];
-        }, py::return_value_policy::reference);
+        }, py::return_value_policy::reference)
+    .def("__getitem__", [](pic::ParticleContainer<D>& s, const py::tuple& indx)
+      {
+        auto ip = indx[0].cast<int>();
+        auto v = indx[1].cast<int>();
+        int nparts = s.size();
+
+        if (ip >= nparts) throw py::index_error();
+        if ( v >= 12)     throw py::index_error();
+
+        if(v == 0) return s.loc(0, ip);
+        if(v == 1) return s.loc(1, ip);
+        if(v == 2) return s.loc(2, ip);
+
+        if(v == 3) return s.vel(0, ip);
+        if(v == 4) return s.vel(1, ip);
+        if(v == 5) return s.vel(2, ip);
+
+        if(v == 6) return s.Epart[ip + 0*nparts];
+        if(v == 7) return s.Epart[ip + 1*nparts];
+        if(v == 8) return s.Epart[ip + 2*nparts];
+
+        if(v == 9) return s.Bpart[ip + 0*nparts];
+        if(v ==10) return s.Bpart[ip + 1*nparts];
+        if(v ==11) return s.Bpart[ip + 2*nparts];
+
+      })
+    .def("__setitem__", [](pic::ParticleContainer<D>& s, const py::tuple& indx, real_prtcl val)
+      {
+        auto ip = indx[0].cast<int>();
+        auto v = indx[1].cast<int>();
+        int nparts = s.size();
+
+        if (ip >= nparts) throw py::index_error();
+        if ( v >= 12)     throw py::index_error();
+
+        if(v == 0) s.loc(0, ip) = val;
+        if(v == 1) s.loc(1, ip) = val;
+        if(v == 2) s.loc(2, ip) = val;
+
+        if(v == 3) s.vel(0, ip) = val;
+        if(v == 4) s.vel(1, ip) = val;
+        if(v == 5) s.vel(2, ip) = val;
+
+        if(v == 6) s.Epart[ip + 0*nparts] = val;
+        if(v == 7) s.Epart[ip + 1*nparts] = val;
+        if(v == 8) s.Epart[ip + 2*nparts] = val;
+
+        if(v == 9) s.Bpart[ip + 0*nparts] = val;
+        if(v ==10) s.Bpart[ip + 1*nparts] = val;
+        if(v ==11) s.Bpart[ip + 2*nparts] = val;
+      });
 
 }
 
@@ -275,7 +328,14 @@ void bind_pic(py::module& m_sub)
   py::class_< pic::Pusher<2,3>> picpusher2d(m_2d, "Pusher");
   picpusher2d
     .def(py::init<>())
-    .def("solve", &pic::Pusher<2,3>::solve);
+    .def_readwrite("bx_ext",  &pic::Pusher<2,3>::bx_ext)
+    .def_readwrite("by_ext",  &pic::Pusher<2,3>::by_ext)
+    .def_readwrite("bz_ext",  &pic::Pusher<2,3>::bz_ext)
+    .def_readwrite("ex_ext",  &pic::Pusher<2,3>::ex_ext)
+    .def_readwrite("ey_ext",  &pic::Pusher<2,3>::ey_ext)
+    .def_readwrite("ez_ext",  &pic::Pusher<2,3>::ez_ext)
+    .def("solve", py::overload_cast<pic::Tile<2>&>(     &pic::Pusher<2,3>::solve))
+    .def("solve", py::overload_cast<pic::Tile<2>&, int>(&pic::Pusher<2,3>::solve));
 
   // Boris pusher
   py::class_<pic::BorisPusher<2,3>>(m_2d, "BorisPusher", picpusher2d)
@@ -304,12 +364,24 @@ void bind_pic(py::module& m_sub)
   py::class_<pic::VayPusher<2,3>>(m_2d, "VayPusher", picpusher2d)
     .def(py::init<>());
 
+  // Hiuera-Cary pusher
+  py::class_<pic::HigueraCaryPusher<2,3>>(m_2d, "HigueraCaryPusher", picpusher2d)
+    .def(py::init<>());
+
 
   // 3D version
   py::class_< pic::Pusher<3,3>> picpusher3d(m_3d, "Pusher");
   picpusher3d
     .def(py::init<>())
-    .def("solve", &pic::Pusher<3,3>::solve);
+    .def_readwrite("bx_ext",  &pic::Pusher<3,3>::bx_ext)
+    .def_readwrite("by_ext",  &pic::Pusher<3,3>::by_ext)
+    .def_readwrite("bz_ext",  &pic::Pusher<3,3>::bz_ext)
+    .def_readwrite("ex_ext",  &pic::Pusher<3,3>::ex_ext)
+    .def_readwrite("ey_ext",  &pic::Pusher<3,3>::ey_ext)
+    .def_readwrite("ez_ext",  &pic::Pusher<3,3>::ez_ext)
+    //.def("solve", &pic::Pusher<3,3>::solve);
+    .def("solve", py::overload_cast<pic::Tile<3>&>(     &pic::Pusher<3,3>::solve))
+    .def("solve", py::overload_cast<pic::Tile<3>&, int>(&pic::Pusher<3,3>::solve));
 
   // Boris pusher
   py::class_<pic::BorisPusher<3,3>>(m_3d, "BorisPusher", picpusher3d)
@@ -319,6 +391,7 @@ void bind_pic(py::module& m_sub)
   py::class_<pic::BorisPusherDrag<3,3>>(m_3d, "BorisDragPusher", picpusher3d)
     .def_readwrite("drag", &pic::BorisPusherDrag<3,3>::drag)
     .def_readwrite("temp", &pic::BorisPusherDrag<3,3>::temp)
+    .def_readwrite("freezing_factor", &pic::BorisPusherDrag<3,3>::freezing_factor)
     .def(py::init<>());
 
   // Boris pusher with radiative pressure force
@@ -337,6 +410,14 @@ void bind_pic(py::module& m_sub)
   py::class_<pic::VayPusher<3,3>>(m_3d, "VayPusher", picpusher3d)
     .def(py::init<>());
 
+  // Higuera-Cary
+  py::class_<pic::HigueraCaryPusher<3,3>>(m_3d, "HigueraCaryPusher", picpusher3d)
+    .def(py::init<>());
+
+  // reduced guiding center approximation
+  py::class_<pic::rGCAPusher<3,3>>(m_3d, "rGCAPusher", picpusher3d)
+    .def(py::init<>());
+
 
   //--------------------------------------------------
 
@@ -421,12 +502,14 @@ void bind_pic(py::module& m_sub)
 
   py::class_<h5io::TestPrtclWriter<2>>(m_2d, "TestPrtclWriter")
     .def(py::init<const std::string&, int, int, int, int, int, int, int, int, int>())
-    .def("write",   &h5io::TestPrtclWriter<2>::write);
+    .def("write",   &h5io::TestPrtclWriter<2>::write)
+    .def_readwrite("ispc", &h5io::TestPrtclWriter<2>::ispc);
 
   // 3D test particles
   py::class_<h5io::TestPrtclWriter<3>>(m_3d, "TestPrtclWriter")
     .def(py::init<const std::string&, int, int, int, int, int, int, int, int, int>())
-    .def("write",   &h5io::TestPrtclWriter<3>::write);
+    .def("write",   &h5io::TestPrtclWriter<3>::write)
+    .def_readwrite("ispc", &h5io::TestPrtclWriter<3>::ispc);
 
   //--------------------------------------------------
   // physical moments of distribution

bindings/pyvlv.c++

@@ -10,7 +10,7 @@
 
 #include "../vlasov/tile.h"
 #include "../em-fields/tile.h"
-#include "../em-fields/damping_tile.h"
+#include "../em-fields/boundaries/damping_tile.h"
 #include "../vlasov/amr/mesh.h"
 #include "../vlasov/amr/refiner.h"
 #include "../vlasov/amr/operators.h"

definitions.h

@@ -1,5 +1,8 @@
 #pragma once
 
+// math constants
+#define PI 3.14159265358979323846 
+
 
 // Old type definitions
 // FIXME: remove completely and switch to the new one below.