First attempt at dg.geo

using pybind11. Documentation and more testing(?) required
So far the Module 3 and parts of 5 are implemented

.gitignore

@@ -3,7 +3,9 @@ tests/__pycache__
 pyfeltor/__pycache__
 pyfeltor/dg/__pycache__
 pyfeltor/dg/create/__pycache__
+pyfeltor/dg/geo/__pycache__
 *.swp
+*.so
 .eggs
 dist
 pyfeltor/_version.py

Makefile

@@ -0,0 +1,27 @@
+device=cpu
+FELTOR_PATH=../feltor
+
+#configure machine
+include $(FELTOR_PATH)/config/default.mk
+include $(FELTOR_PATH)/config/*.mk
+include $(FELTOR_PATH)/config/devices/devices.mk
+JSONLIB=-DDG_USE_JSONHPP
+
+INCLUDE+=-I$(FELTOR_PATH)/inc/
+
+PRE=pyfeltor/dg/geo
+
+all: geometries polynomial solovev guenter circular flux mod toroidal utility
+
+# the cpp file can have different name, the so file needs the name of the python module
+
+%: $(PRE)/%.cpp
+	$(CC) $(OPT) $(CFLAGS) -shared -fPIC $$(python3 -m pybind11 --includes) $< -o $(PRE)/$@$$(python3-config --extension-suffix) $(JSONLIB) $(INCLUDE) -g
+
+install:
+	python3 -m pip install -e .
+
+.PHONY: clean
+
+clean:
+	rm -f $(PRE)/geometries/$$(python3-config --extension-suffix)

pyfeltor/dg/__init__.py

@@ -5,3 +5,7 @@
 from .enums import bc, direction, inverse_bc, inverse_dir
 from .grid import Grid
 from .evaluation import evaluate, integrate
+try:
+    from . import geo
+except ImportError:
+    pass

pyfeltor/dg/geo/__init__.py

@@ -0,0 +1,8 @@
+# docstring displayed by help(pyfeltor)
+""" The python version of the dg library
+"""
+
+from .geometries import *
+from .utility import *
+from .flux import *
+from . import polynomial, solovev, guenter, toroidal, circular, mod

pyfeltor/dg/geo/circular.cpp

@@ -0,0 +1,29 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+
+#include "dg/geometries/geometries.h"
+namespace py = pybind11;
+namespace poly = dg::geo::circular;
+
+
+PYBIND11_MODULE(circular, m) { // name of the python module, not the C++ file
+    m.doc() = "pybind11 binding guenter field"; // optional module docstring
+    py::class_<poly::Psip>(m,"Psip")
+        .def(py::init<double,double,double>( ))
+        .def( "__call__", py::vectorize([]( poly::Psip& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::PsipR>(m,"PsipR")
+        .def(py::init<double,double>( ))
+        .def( "__call__", py::vectorize([]( poly::PsipR& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::PsipZ>(m,"PsipZ")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::PsipZ& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    m.def( "createPsip", &dg::geo::circular::createPsip);
+    m.def( "createIpol", &dg::geo::circular::createIpol);
+}

pyfeltor/dg/geo/flux.cpp

@@ -0,0 +1,27 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+
+#include <pybind11_json/pybind11_json.hpp>
+
+#include "dg/algorithm.h"
+#include "dg/file/json_utilities.h"
+#include "dg/geometries/geometries.h"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(flux, m) {
+    m.def( "createSolovevField", &dg::geo::createSolovevField);
+    m.def( "createPolynomialField", &dg::geo::createPolynomialField);
+    m.def( "createGuenterField", &dg::geo::createGuenterField);
+    m.def( "createToroidalField", &dg::geo::createToroidalField);
+    m.def( "createCircularField", &dg::geo::createCircularField);
+    m.def( "createModifiedField", []( const nlohmann::json& gs, const nlohmann::json&
+                jsmod, dg::geo::CylindricalFunctor& wall,
+                dg::geo::CylindricalFunctor& transition)
+            { return dg::geo::createModifiedField( gs, jsmod, wall, transition);});
+    m.def( "createWallRegion", []( const nlohmann::json& gs, const nlohmann::json&
+                jsmod) { return dg::geo::createWallRegion( gs, jsmod);});
+    m.def( "createMagneticField", []( const nlohmann::json& json){
+            return dg::geo::createMagneticField( json); });
+}

pyfeltor/dg/geo/geometries.cpp

@@ -0,0 +1,213 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+
+#include <pybind11_json/pybind11_json.hpp>
+
+#include "dg/algorithm.h"
+#include "dg/file/json_utilities.h"
+#include "dg/geometries/geometries.h"
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(geometries, m) { // name of the python module, not the C++ file
+    py::enum_<dg::geo::description>( m, "description")
+        .value("standardO", dg::geo::description::standardO)
+        .value("standardX", dg::geo::description::standardX)
+        .value("doubleX", dg::geo::description::doubleX)
+        .value("none", dg::geo::description::none)
+        .value("square", dg::geo::description::square)
+        .value("centeredX", dg::geo::description::centeredX)
+        ;
+    py::enum_<dg::geo::equilibrium>( m, "equilibrium")
+        .value("solovev", dg::geo::equilibrium::solovev)
+        .value("taylor", dg::geo::equilibrium::taylor)
+        .value("polynomial", dg::geo::equilibrium::polynomial)
+        .value("guenter", dg::geo::equilibrium::guenter)
+        .value("toroidal", dg::geo::equilibrium::toroidal)
+        .value("circular", dg::geo::equilibrium::circular)
+        ;
+    py::enum_<dg::geo::modifier>( m, "modifier")
+        .value("none", dg::geo::modifier::none)
+        .value("heaviside", dg::geo::modifier::heaviside)
+        .value("sol_pfr", dg::geo::modifier::sol_pfr)
+        ;
+    py::class_<dg::geo::MagneticFieldParameters>(m,"MagneticFieldParameters")
+        .def(py::init<>())
+        .def(py::init<double,double,double,dg::geo::equilibrium,dg::geo::modifier,dg::geo::description>())
+        .def("a",&dg::geo::MagneticFieldParameters::a)
+        .def("elongation",&dg::geo::MagneticFieldParameters::elongation)
+        .def("triangularity",&dg::geo::MagneticFieldParameters::triangularity)
+        .def("getEquilibrium",&dg::geo::MagneticFieldParameters::getEquilibrium)
+        .def("getModifier",&dg::geo::MagneticFieldParameters::getModifier)
+        .def("getDescription",&dg::geo::MagneticFieldParameters::getDescription)
+        ;
+    py::class_<dg::geo::TokamakMagneticField>(m,"TokamakMagneticField")
+        .def( "R0", &dg::geo::TokamakMagneticField::R0)
+        .def( "psip", &dg::geo::TokamakMagneticField::psip)
+        .def( "psipR", &dg::geo::TokamakMagneticField::psipR)
+        .def( "psipZ", &dg::geo::TokamakMagneticField::psipZ)
+        .def( "psipRR", &dg::geo::TokamakMagneticField::psipRR)
+        .def( "psipRZ", &dg::geo::TokamakMagneticField::psipRZ)
+        .def( "psipZZ", &dg::geo::TokamakMagneticField::psipZZ)
+        .def( "ipol", &dg::geo::TokamakMagneticField::ipol)
+        .def( "ipolR", &dg::geo::TokamakMagneticField::ipolR)
+        .def( "ipolZ", &dg::geo::TokamakMagneticField::ipolZ)
+        .def( "get_psip", &dg::geo::TokamakMagneticField::get_psip)
+        .def( "get_ipol", &dg::geo::TokamakMagneticField::get_ipol)
+        .def( "params", &dg::geo::TokamakMagneticField::params)
+    ;
+    py::class_<dg::geo::LaplacePsip>(m,"LaplacePsip")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::LaplacePsip& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::Bmodule>(m,"Bmodule")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::Bmodule& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::InvB>(m,"InvB")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::InvB& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::LnB>(m,"LnB")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::LnB& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BR>(m,"BR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BZ>(m,"BZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::CurvatureNablaBR>(m,"CurvatureNablaBR")
+        .def(py::init<const dg::geo::TokamakMagneticField&,int>())
+        .def( "__call__", py::vectorize([]( dg::geo::CurvatureNablaBR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::CurvatureNablaBZ>(m,"CurvatureNablaBZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&,int>())
+        .def( "__call__", py::vectorize([]( dg::geo::CurvatureNablaBZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::CurvatureKappaR>(m,"CurvatureKappaR")
+        .def(py::init<const dg::geo::TokamakMagneticField&,int>())
+        .def( "__call__", py::vectorize([]( dg::geo::CurvatureKappaR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::CurvatureKappaZ>(m,"CurvatureKappaZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&,int>())
+        .def( "__call__", py::vectorize([]( dg::geo::CurvatureKappaZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::DivCurvatureKappa>(m,"DivCurvatureKappa")
+        .def(py::init<const dg::geo::TokamakMagneticField&,int>())
+        .def( "__call__", py::vectorize([]( dg::geo::DivCurvatureKappa& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::DivCurvatureNablaB>(m,"DivCurvatureNablaB")
+        .def(py::init<const dg::geo::TokamakMagneticField&,int>())
+        .def( "__call__", py::vectorize([]( dg::geo::DivCurvatureNablaB& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueCurvatureNablaBR>(m,"TrueCurvatureNablaBR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueCurvatureNablaBR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueCurvatureNablaBZ>(m,"TrueCurvatureNablaBZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueCurvatureNablaBZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueCurvatureNablaBP>(m,"TrueCurvatureNablaBP")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueCurvatureNablaBP& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueCurvatureKappaR>(m,"TrueCurvatureKappaR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueCurvatureKappaR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueCurvatureKappaZ>(m,"TrueCurvatureKappaZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueCurvatureKappaZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueCurvatureKappaP>(m,"TrueCurvatureKappaP")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueCurvatureKappaP& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueDivCurvatureKappa>(m,"TrueDivCurvatureKappa")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueDivCurvatureKappa& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::TrueDivCurvatureNablaB>(m,"TrueDivCurvatureNablaB")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::TrueDivCurvatureNablaB& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::GradLnB>(m,"GradLnB")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::GradLnB& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::Divb>(m,"Divb")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::Divb& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BFieldP>(m,"BFieldP")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BFieldP& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BFieldR>(m,"BFieldR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BFieldR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BFieldZ>(m,"BFieldZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BFieldZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BFieldT>(m,"BFieldT")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BFieldT& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatR>(m,"BHatR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatZ>(m,"BHatZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatP>(m,"BHatP")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatP& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatRR>(m,"BHatRR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatRR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatRZ>(m,"BHatRZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatRZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatZR>(m,"BHatZR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatZR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatZZ>(m,"BHatZZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatZZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatPR>(m,"BHatPR")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatPR& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::BHatPZ>(m,"BHatPZ")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::BHatPZ& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::DivVVP>(m,"DivVVP")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::DivVVP& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::RhoP>(m,"RhoP")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::RhoP& my,
+                        double R, double Z){ return my(R,Z);}));
+    py::class_<dg::geo::Hoo>(m,"Hoo")
+        .def(py::init<const dg::geo::TokamakMagneticField&>())
+        .def( "__call__", py::vectorize([]( dg::geo::Hoo& my,
+                        double R, double Z){ return my(R,Z);}));
+
+}

pyfeltor/dg/geo/guenter.cpp

@@ -0,0 +1,59 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+
+#include "dg/geometries/geometries.h"
+namespace py = pybind11;
+namespace poly = dg::geo::guenter;
+
+
+PYBIND11_MODULE(guenter, m) { // name of the python module, not the C++ file
+    m.doc() = "pybind11 binding guenter field"; // optional module docstring
+    py::class_<poly::Psip>(m,"Psip")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::Psip& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::PsipR>(m,"PsipR")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::PsipR& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::PsipZ>(m,"PsipZ")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::PsipZ& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::PsipRR>(m,"PsipRR")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::PsipRR& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::PsipRZ>(m,"PsipRZ")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::PsipRZ& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::PsipZZ>(m,"PsipZZ")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::PsipZZ& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::Ipol>(m,"Ipol")
+        .def(py::init<double>( ))
+        .def( "__call__", py::vectorize([]( poly::Ipol& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::IpolR>(m,"IpolR")
+        .def(py::init<>( ))
+        .def( "__call__", py::vectorize([]( poly::IpolR& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    py::class_<poly::IpolZ>(m,"IpolZ")
+        .def(py::init<>( ))
+        .def( "__call__", py::vectorize([]( poly::IpolZ& my,
+                        double R, double Z){ return my(R,Z);}))
+        ;
+    m.def( "createPsip", &dg::geo::guenter::createPsip);
+    m.def( "createIpol", &dg::geo::guenter::createIpol);
+}