From 946570241e1eff8cd38776a30ac092f028b5128e Mon Sep 17 00:00:00 2001 From: Matthias Date: Wed, 9 Apr 2025 11:12:18 +0200 Subject: [PATCH 01/16] Change fully to pyproject.toml --- README.md | 2 +- pyproject.toml | 67 +++++++++++++++++++++++++++++++++++++++++++++++--- setup.cfg | 32 ------------------------ setup.py | 6 ----- 4 files changed, 64 insertions(+), 43 deletions(-) delete mode 100644 setup.cfg delete mode 100644 setup.py diff --git a/README.md b/README.md index 44b6317..8b76ac2 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,7 @@ unparallelized in python. As a second addition, now also Feltor's geometries extension is available in python. However, the geometries functions and classes are not re-implemented in python, but -they are bound to python via the [pybind11](https://github.com/pybind11/pybind11) +they are bound to python via the [pybind11](https://github.com/pybind/pybind11) library. As such the corresponding C++ binding code must be compiled in order to generate the module `dg.geo`. ## Installation diff --git a/pyproject.toml b/pyproject.toml index 3db995f..2c1b474 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,65 @@ -[build-system] -requires = ["setuptools", "wheel"] -build-backend = "setuptools.build_meta" +[project] +name = "pyfeltor" +authors = [ + {name = "Matthias Wiesenberger", email = "mattwi@fysik.dtu.dk"}, +] +description = "A subset of Feltor's dg library as a python implementation" +license = {file = "LICENSE"} +readme = "README.md" +requires-python = ">=3.12" # see scientific-python.org/specs/spec-0000 : +dependencies = [ + "scipy >= 1.10.0", + "numpy >= 1.24.0", +] +classifiers = [ + "Intended Audience :: Science/Research", + "Intended Audience :: Education", + "Intended Audience :: Developers", + "Natural Language :: English", + "Operating System :: POSIX :: Linux", + "Programming Language :: Python :: 3", +] +dynamic = ["version"] # build backend fills this field +[project.urls] +Source = "https://github.com/feltor-dev/pyfeltor" +Tracker = "https://github.com/feltor-dev/pyfeltor/issues" +[build-system] # which build system to use +requires = [ +# "wheel", + "setuptools >= 65", # build with setuptools + "setuptools-scm >= 8.0.0", # for automatic version from git tags +] +build-backend = "setuptools.build_meta" [tool.setuptools_scm] -write_to = "pyfeltor/_version.py" +# Activate setuptools-scm + + +[project.optional-dependencies] +# https://docs.astral.sh/ruff/ +# pip install -e .[lint] +lint = [ + "ruff", # a python linter and formatter +] +# Ruff rules +# https://docs.astral.sh/ruff/rules/ +[tool.ruff.lint] +select = [ + "F", # Pyflakes + "E", # pycodestyle error + "W", # pycodestyle warning + "I", # isort + "PT", # flake8-pytest-style + "B", # flake8-bugbear + "A", # flake8-builtins + "C4", # flake8-comprehensions + "T10", # flake8-debugger + # "PTH", # flake8-use-pathlib + "SIM", # flake8-simplify + "TCH", # flake8-type-checking + "UP", # pyupgrade + "FURB", # refurb + "PERF", # perflint + "RUF", # ruff specific + "NPY", # NumPy specific +] diff --git a/setup.cfg b/setup.cfg deleted file mode 100644 index 2dbce6b..0000000 --- a/setup.cfg +++ /dev/null @@ -1,32 +0,0 @@ -[metadata] -name = pyfeltor -url = https://github.com/feltor-dev/pyFeltor -author = Matthias Wiesenberger -author_email = mattwi@fysik.dtu.dk -description = A subset of Feltor's dg library as a python implementation -license = MIT License -python_requires >= 3.8 -long_description = file: README.md -classifiers = - Intended Audience :: Science/Research - Intended Audience :: Education - Intended Audience :: Developers - Natural Language :: English - Operating System :: POSIX :: Linux - Programming Language :: Python :: 3.9 - Programming Language :: Python :: 3.8 - -[options] -setup_requires = - setuptools>=58 - setuptools_scm -install_requires = - scipy >= 1.7.1 - numpy >= 1.21.1 -include_package_data = True -packages = find: - -[options.extras_require] -calc = - numpy >= 1.21.1 -docs = sphinx >= 1.4 diff --git a/setup.py b/setup.py deleted file mode 100644 index 6012477..0000000 --- a/setup.py +++ /dev/null @@ -1,6 +0,0 @@ -#!/usr/bin/env python3 - -from setuptools import setup - -if __name__ == "__main__": - setup() From cd287a7d6ac7e1ef4f696aa3fb7ce1ba427defc3 Mon Sep 17 00:00:00 2001 From: Matthias Date: Tue, 29 Apr 2025 11:02:16 +0200 Subject: [PATCH 02/16] First steps in pyproject toml --- pyproject.toml | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 2c1b474..262c6fa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,13 +26,11 @@ Tracker = "https://github.com/feltor-dev/pyfeltor/issues" [build-system] # which build system to use requires = [ -# "wheel", - "setuptools >= 65", # build with setuptools - "setuptools-scm >= 8.0.0", # for automatic version from git tags + "scikit-build-core >= 0.10", + "pybind11" ] -build-backend = "setuptools.build_meta" -[tool.setuptools_scm] -# Activate setuptools-scm +build-backend = "scikit_build_core.build" +# TODO How to get automatic versioning in scikit? [project.optional-dependencies] From 2f229bd70a99e421767662c2c79e0e37158f0f4c Mon Sep 17 00:00:00 2001 From: Matthias Date: Thu, 12 Jun 2025 12:11:34 +0200 Subject: [PATCH 03/16] Write CMakeLists.txt using the mwiesenberger/feltor cmake branch Works and should be updated when ready README should also be updated when ready --- CMakeLists.txt | 70 +++++++++++++++++++++++++++++++++++++ pyfeltor/dg/geo/solovev.cpp | 1 + pyproject.toml | 8 +++-- 3 files changed, 76 insertions(+), 3 deletions(-) create mode 100644 CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..1b3b808 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,70 @@ +## Inspired from https://github.com/pybind/scikit_build_example/blob/master/CMakeLists.txt + +# Require CMake 3.15+ (matching scikit-build-core) Use new versions of all +# policies up to CMake 3.27 +cmake_minimum_required( VERSION "3.26") + + +# Scikit-build-core sets these values for you, or you can just hard-code the +# name and version. +project( + ${SKBUILD_PROJECT_NAME} + VERSION ${SKBUILD_PROJECT_VERSION} + LANGUAGES CXX) + +include(cmake/CPM.cmake) +set( CCCL_THRUST_DEVICE_SYSTEM "CPP") +CPMAddPackage( + NAME feltor + GITHUB_REPOSITORY "mwiesenberger/feltor" + GIT_TAG cmake + SYSTEM ON + EXCLUDE_FROM_ALL ON + OPTIONS "FELTOR_DG_WITH_MATRIX OFF" "FELTOR_FILE_WITH_NETCDF OFF" +) + + +# Find the module development requirements (requires FindPython from 3.17 or +# scikit-build-core's built-in backport) +# MW: https://cmake.org/cmake/help/latest/module/FindPython.html +find_package(Python REQUIRED COMPONENTS Interpreter Development.Module) +# Install pyprojects pybind11 +find_package(pybind11 CONFIG REQUIRED) + +# Needs to come after feltor because json is already included there +# Should we include guard feltor's CPM? +find_package(pybind11_json CONFIG REQUIRED) + +set( PYFELTOR_EXE + pyfeltor/dg/geo/geometries.cpp + pyfeltor/dg/geo/polynomial.cpp + pyfeltor/dg/geo/solovev.cpp + pyfeltor/dg/geo/guenter.cpp + pyfeltor/dg/geo/circular.cpp + pyfeltor/dg/geo/flux.cpp + pyfeltor/dg/geo/mod.cpp + pyfeltor/dg/geo/toroidal.cpp + pyfeltor/dg/geo/utility.cpp +) + +# Add a library using FindPython's tooling (pybind11 also provides a helper like +# this) +python_add_library(_core MODULE WITH_SOABI ${PYFELTOR_EXE}) +target_compile_features(_core PRIVATE cxx_std_17) +if (MSVC) + # warning level 4 + target_compile_options(_core PRIVATE /W4 /permissive- /arch:AVX2 /O2 /EHsc -DFP_FAST_FMA -D_USE_MATH_DEFINES) +else() + # additional warnings + target_compile_options(_core PRIVATE -Wall -Wextra -Wpedantic -march=native -O2) +endif() + +target_link_libraries(_core PRIVATE pybind11::headers) + + +target_link_libraries(_core PRIVATE feltor::dg::geometries) +target_link_libraries(_core PRIVATE pybind11_json) + + +# The install directory is the output (wheel) directory +install(TARGETS _core DESTINATION scikit_build_example) diff --git a/pyfeltor/dg/geo/solovev.cpp b/pyfeltor/dg/geo/solovev.cpp index 4432046..26b81b4 100644 --- a/pyfeltor/dg/geo/solovev.cpp +++ b/pyfeltor/dg/geo/solovev.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include "dg/geometries/geometries.h" namespace py = pybind11; diff --git a/pyproject.toml b/pyproject.toml index 262c6fa..6fcf185 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,7 +1,7 @@ [project] name = "pyfeltor" authors = [ - {name = "Matthias Wiesenberger", email = "mattwi@fysik.dtu.dk"}, + {name = "Matthias Wiesenberger", email = "mattwi@dtu.dk"}, ] description = "A subset of Feltor's dg library as a python implementation" license = {file = "LICENSE"} @@ -19,7 +19,7 @@ classifiers = [ "Operating System :: POSIX :: Linux", "Programming Language :: Python :: 3", ] -dynamic = ["version"] # build backend fills this field +version = "0.1.0" # TODO How to get automatic versioning in scikit? [project.urls] Source = "https://github.com/feltor-dev/pyfeltor" Tracker = "https://github.com/feltor-dev/pyfeltor/issues" @@ -30,7 +30,9 @@ requires = [ "pybind11" ] build-backend = "scikit_build_core.build" -# TODO How to get automatic versioning in scikit? + +[tool.scikit-build] +build.verbose = true # cmake --build . -v ## only visible if pip install . -v is used [project.optional-dependencies] From ca3a9cfca457edb38e4eb7d88b09ce8656bb7fc8 Mon Sep 17 00:00:00 2001 From: Matthias Date: Mon, 23 Jun 2025 13:31:55 +0200 Subject: [PATCH 04/16] Use released feltor version 8.2 --- CMakeLists.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1b3b808..2bf49fb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,8 +16,8 @@ include(cmake/CPM.cmake) set( CCCL_THRUST_DEVICE_SYSTEM "CPP") CPMAddPackage( NAME feltor - GITHUB_REPOSITORY "mwiesenberger/feltor" - GIT_TAG cmake + GITHUB_REPOSITORY "feltor-dev/feltor" + VERSION 8.2 SYSTEM ON EXCLUDE_FROM_ALL ON OPTIONS "FELTOR_DG_WITH_MATRIX OFF" "FELTOR_FILE_WITH_NETCDF OFF" @@ -32,7 +32,7 @@ find_package(Python REQUIRED COMPONENTS Interpreter Development.Module) find_package(pybind11 CONFIG REQUIRED) # Needs to come after feltor because json is already included there -# Should we include guard feltor's CPM? +# unless we use CPM_USE_LOCAL_PACKAGES find_package(pybind11_json CONFIG REQUIRED) set( PYFELTOR_EXE From 4eef5ad1c04add710d032b6ccbccdd10adb8b0a6 Mon Sep 17 00:00:00 2001 From: Matthias Date: Mon, 23 Jun 2025 13:39:05 +0200 Subject: [PATCH 05/16] Add missing CPM.cmake --- cmake/CPM.cmake | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 cmake/CPM.cmake diff --git a/cmake/CPM.cmake b/cmake/CPM.cmake new file mode 100644 index 0000000..77b333c --- /dev/null +++ b/cmake/CPM.cmake @@ -0,0 +1,24 @@ +# SPDX-License-Identifier: MIT +# +# SPDX-FileCopyrightText: Copyright (c) 2019-2023 Lars Melchior and contributors + +set(CPM_DOWNLOAD_VERSION 0.41.0) +set(CPM_HASH_SUM "e570f03806b9aae2082ca5b950a9e6b3b41ad56972a78a876aedcaad16653116") + +if(CPM_SOURCE_CACHE) + set(CPM_DOWNLOAD_LOCATION "${CPM_SOURCE_CACHE}/cpm/CPM_${CPM_DOWNLOAD_VERSION}.cmake") +elseif(DEFINED ENV{CPM_SOURCE_CACHE}) + set(CPM_DOWNLOAD_LOCATION "$ENV{CPM_SOURCE_CACHE}/cpm/CPM_${CPM_DOWNLOAD_VERSION}.cmake") +else() + set(CPM_DOWNLOAD_LOCATION "${CMAKE_BINARY_DIR}/cmake/CPM_${CPM_DOWNLOAD_VERSION}.cmake") +endif() + +# Expand relative path. This is important if the provided path contains a tilde (~) +get_filename_component(CPM_DOWNLOAD_LOCATION ${CPM_DOWNLOAD_LOCATION} ABSOLUTE) + +file(DOWNLOAD + https://github.com/cpm-cmake/CPM.cmake/releases/download/v${CPM_DOWNLOAD_VERSION}/CPM.cmake + ${CPM_DOWNLOAD_LOCATION} EXPECTED_HASH SHA256=${CPM_HASH_SUM} +) + +include(${CPM_DOWNLOAD_LOCATION}) From 10e9d67f5d2a659d05ef3daad4ae88630f7b889d Mon Sep 17 00:00:00 2001 From: Matthias Date: Mon, 23 Jun 2025 14:45:28 +0200 Subject: [PATCH 06/16] Add pybind11_json via CPM --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2bf49fb..ab0fc68 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,7 +33,7 @@ find_package(pybind11 CONFIG REQUIRED) # Needs to come after feltor because json is already included there # unless we use CPM_USE_LOCAL_PACKAGES -find_package(pybind11_json CONFIG REQUIRED) +CPMAddPackage( "gh:pybind/pybind11_json#0.2.15") set( PYFELTOR_EXE pyfeltor/dg/geo/geometries.cpp From 4d4535e9a3ba9c881a9cc1dafd14df31f4d113bf Mon Sep 17 00:00:00 2001 From: Matthias Date: Mon, 23 Jun 2025 16:34:58 +0200 Subject: [PATCH 07/16] Use pytest-datadir for json files pytest can be run from any directory then --- tests/test_geometries.py | 15 +++++++-------- tests/{ => test_geometries}/enrx_tcv.json | 0 .../geometry_params_Xpoint.json | 0 3 files changed, 7 insertions(+), 8 deletions(-) rename tests/{ => test_geometries}/enrx_tcv.json (100%) rename tests/{ => test_geometries}/geometry_params_Xpoint.json (100%) diff --git a/tests/test_geometries.py b/tests/test_geometries.py index 7edb63a..9c31b3b 100644 --- a/tests/test_geometries.py +++ b/tests/test_geometries.py @@ -1,8 +1,8 @@ import pytest import json import numpy as np -import magneticfielddb as magdb from pyfeltor import dg +import pytest_datadir geo_loaded = True try: @@ -35,9 +35,9 @@ def test_polynomial(): print( psi_) -def test_make_field(): +def test_make_field(datadir): if not geo_exists(): return - with open ("geometry_params_Xpoint.json", "r") as f: + with open (datadir / "geometry_params_Xpoint.json", "r") as f: magparams = json.load(f) mag = dg.geo.createMagneticField( magparams) a = mag.params().a() @@ -54,10 +54,9 @@ def test_make_field(): -def test_q_profile(): +def test_q_profile(datadir): if not geo_exists(): return - #magparams = magdb.select( "COMPASS/compass_1X.json") - with open ("enrx_tcv.json", "r") as f: + with open (datadir / "enrx_tcv.json", "r") as f: magparams = json.load(f) mag = dg.geo.createMagneticField(magparams) qfunctor = dg.geo.SafetyFactor(mag) @@ -68,9 +67,9 @@ def test_q_profile(): psi_values = np.linspace( psipO, 0, 20, endpoint = False) print(qfunctor(psi_values)) -def test_sheath(): +def test_sheath(datadir): if not geo_exists(): return - with open ("enrx_tcv.json", "r") as f: + with open (datadir / "enrx_tcv.json", "r") as f: magparams = json.load(f) mag = dg.geo.createMagneticField(magparams) wall = {"type" : "sol_pfr", "alpha": [0.05,0.05], "boundary" : [1.09,0.97]} diff --git a/tests/enrx_tcv.json b/tests/test_geometries/enrx_tcv.json similarity index 100% rename from tests/enrx_tcv.json rename to tests/test_geometries/enrx_tcv.json diff --git a/tests/geometry_params_Xpoint.json b/tests/test_geometries/geometry_params_Xpoint.json similarity index 100% rename from tests/geometry_params_Xpoint.json rename to tests/test_geometries/geometry_params_Xpoint.json From ebddc5d50963eb415f10804e2f4dc19316b1114d Mon Sep 17 00:00:00 2001 From: Matthias Date: Mon, 23 Jun 2025 16:37:14 +0200 Subject: [PATCH 08/16] Add cmake policy version minimum for cmake > 4.0 --- pyproject.toml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 6fcf185..137b7c0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,9 +33,15 @@ build-backend = "scikit_build_core.build" [tool.scikit-build] build.verbose = true # cmake --build . -v ## only visible if pip install . -v is used +cmake.args = ["-DCMAKE_POLICY_VERSION_MINIMUM=3.5"] # Necessary for nlohmann_json with cmake>4.0 [project.optional-dependencies] +test = [ + "pytest", + "pytest-datadir" # For json files during geometries test +] + # https://docs.astral.sh/ruff/ # pip install -e .[lint] lint = [ From f8bdb94b77edb63e3bc832aac902aef35428c443 Mon Sep 17 00:00:00 2001 From: Matthias Date: Mon, 23 Jun 2025 16:38:10 +0200 Subject: [PATCH 09/16] Fix CMakeLists to build and install *.so to dg/geo --- CMakeLists.txt | 53 +++++++++++++++++++++++++++----------------------- 1 file changed, 29 insertions(+), 24 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ab0fc68..614b61b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -35,36 +35,41 @@ find_package(pybind11 CONFIG REQUIRED) # unless we use CPM_USE_LOCAL_PACKAGES CPMAddPackage( "gh:pybind/pybind11_json#0.2.15") -set( PYFELTOR_EXE - pyfeltor/dg/geo/geometries.cpp - pyfeltor/dg/geo/polynomial.cpp - pyfeltor/dg/geo/solovev.cpp - pyfeltor/dg/geo/guenter.cpp - pyfeltor/dg/geo/circular.cpp - pyfeltor/dg/geo/flux.cpp - pyfeltor/dg/geo/mod.cpp - pyfeltor/dg/geo/toroidal.cpp - pyfeltor/dg/geo/utility.cpp + +# Follow the Makefile approach and compile each cpp into its own so lib + +set( PYFELTOR_LIBS + geometries + polynomial + solovev + guenter + circular + flux + mod + toroidal + utility ) # Add a library using FindPython's tooling (pybind11 also provides a helper like # this) -python_add_library(_core MODULE WITH_SOABI ${PYFELTOR_EXE}) -target_compile_features(_core PRIVATE cxx_std_17) -if (MSVC) - # warning level 4 - target_compile_options(_core PRIVATE /W4 /permissive- /arch:AVX2 /O2 /EHsc -DFP_FAST_FMA -D_USE_MATH_DEFINES) -else() - # additional warnings - target_compile_options(_core PRIVATE -Wall -Wextra -Wpedantic -march=native -O2) -endif() +foreach( lib IN LISTS PYFELTOR_LIBS) + python_add_library(${lib} MODULE WITH_SOABI "pyfeltor/dg/geo/${lib}.cpp") + target_compile_features(${lib} PRIVATE cxx_std_17) + if (MSVC) + # warning level 4 + target_compile_options(${lib} PRIVATE /W4 /permissive- /arch:AVX2 /O2 /EHsc -DFP_FAST_FMA -D_USE_MATH_DEFINES) + else() + # additional warnings + target_compile_options(${lib} PRIVATE -Wall -Wextra -Wpedantic -march=native -O2) + endif() -target_link_libraries(_core PRIVATE pybind11::headers) + target_link_libraries(${lib} PRIVATE pybind11::headers) -target_link_libraries(_core PRIVATE feltor::dg::geometries) -target_link_libraries(_core PRIVATE pybind11_json) + target_link_libraries(${lib} PRIVATE feltor::dg::geometries) + target_link_libraries(${lib} PRIVATE pybind11_json) + # The install directory is the output (wheel) directory + install(TARGETS ${lib} DESTINATION pyfeltor/dg/geo) +endforeach() -# The install directory is the output (wheel) directory -install(TARGETS _core DESTINATION scikit_build_example) From c6b546e6d8371e584b22f4e60ba2889c646f887a Mon Sep 17 00:00:00 2001 From: Matthias Wiesenberger Date: Tue, 24 Jun 2025 14:01:17 +0200 Subject: [PATCH 10/16] Fix flux.cpp py import so it works on Windows Also fix test_geometries so it passes on Windows --- .gitignore | 9 +++++++++ pyfeltor/dg/geo/flux.cpp | 2 -- tests/test_geometries.py | 4 ++-- 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/.gitignore b/.gitignore index 0fa872d..954c84f 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,12 @@ pyfeltor/dg/geo/__pycache__ .eggs dist pyfeltor/_version.py + +uv.lock + +# ignore visual studio .vs table +/.vs + +# CMake output +build +out \ No newline at end of file diff --git a/pyfeltor/dg/geo/flux.cpp b/pyfeltor/dg/geo/flux.cpp index 328a4a2..e39346e 100644 --- a/pyfeltor/dg/geo/flux.cpp +++ b/pyfeltor/dg/geo/flux.cpp @@ -10,8 +10,6 @@ namespace py = pybind11; -// import part of itself!? Does this always work?? Wow -py::object grid = py::module_::import( "pyfeltor.dg").attr("Grid"); PYBIND11_MODULE(flux, m) { m.def( "createSolovevField", &dg::geo::createSolovevField); diff --git a/tests/test_geometries.py b/tests/test_geometries.py index 9c31b3b..b514b03 100644 --- a/tests/test_geometries.py +++ b/tests/test_geometries.py @@ -86,8 +86,8 @@ def test_sheath(datadir): point,RO,ZO = dg.geo.findOpoint( mag.get_psip(), R0, 0) assert( point == 1) - assert( RO == 0.9022324795511738) - assert( ZO == -0.005684720497897021) + assert( abs(RO - 0.9022324795511738) < 1e-10) + assert( abs(ZO - (-0.005684720497897021)) < 1e-10) closed = dg.geo.mod.ClosedFieldlineRegion( mag); assert( closed( 0, R0)) From 139140e17656117e225dfeee81efa5fb149fe8be Mon Sep 17 00:00:00 2001 From: Matthias Date: Tue, 24 Jun 2025 14:10:45 +0200 Subject: [PATCH 11/16] Remove include guard for dg.geo Now it is always included It is not possible to install without any more ... --- pyfeltor/dg/__init__.py | 5 +---- tests/guenter_params.json | 7 ------- tests/test_geometries.py | 18 ------------------ 3 files changed, 1 insertion(+), 29 deletions(-) delete mode 100644 tests/guenter_params.json diff --git a/pyfeltor/dg/__init__.py b/pyfeltor/dg/__init__.py index a87aa10..ba873be 100644 --- a/pyfeltor/dg/__init__.py +++ b/pyfeltor/dg/__init__.py @@ -5,7 +5,4 @@ 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 +from . import geo diff --git a/tests/guenter_params.json b/tests/guenter_params.json deleted file mode 100644 index ed5b61d..0000000 --- a/tests/guenter_params.json +++ /dev/null @@ -1,7 +0,0 @@ -{ - "equilibrium" : "guenter", - "description" : "square", - "I_0" : 10.0, - "R_0" : 3.0, - "comment" : "Sine functions; Guenter et al, J. Comput. Phys. 209 (2005) 354" -} diff --git a/tests/test_geometries.py b/tests/test_geometries.py index b514b03..d3ca102 100644 --- a/tests/test_geometries.py +++ b/tests/test_geometries.py @@ -4,23 +4,8 @@ from pyfeltor import dg import pytest_datadir -geo_loaded = True -try: - from pyfeltor.dg import geo -except ImportError: - geo_loaded = False - - -# Run with pytest-3 -s . to see stdout output -def geo_exists(): - if not geo_loaded: - print( "dg.geo not compiled") - return False - return True - def test_polynomial(): - if not geo_exists(): return c = np.array( [1,2,3,4]) params = {"R_0" : 400, "inverseaspectratio" : 20, "elongation" : 1, "triangularity" : 1, "PP" : 1, "PI" : 1, "description" : "standardX", "M" : 2, "N" : 2, "c" : c.tolist()} @@ -36,7 +21,6 @@ def test_polynomial(): def test_make_field(datadir): - if not geo_exists(): return with open (datadir / "geometry_params_Xpoint.json", "r") as f: magparams = json.load(f) mag = dg.geo.createMagneticField( magparams) @@ -55,7 +39,6 @@ def test_make_field(datadir): def test_q_profile(datadir): - if not geo_exists(): return with open (datadir / "enrx_tcv.json", "r") as f: magparams = json.load(f) mag = dg.geo.createMagneticField(magparams) @@ -68,7 +51,6 @@ def test_q_profile(datadir): print(qfunctor(psi_values)) def test_sheath(datadir): - if not geo_exists(): return with open (datadir / "enrx_tcv.json", "r") as f: magparams = json.load(f) mag = dg.geo.createMagneticField(magparams) From 4bf6c183a9a705b04f5fcc86aa266bf5905d637c Mon Sep 17 00:00:00 2001 From: Matthias Wiesenberger Date: Tue, 24 Jun 2025 14:22:53 +0200 Subject: [PATCH 12/16] Update the README --- README.md | 33 ++------------------------------- 1 file changed, 2 insertions(+), 31 deletions(-) diff --git a/README.md b/README.md index 8b76ac2..cfc3b11 100644 --- a/README.md +++ b/README.md @@ -20,14 +20,13 @@ through python this enables many applications beyond simple simulations diagnostics. The only downside is of course that all functions are unparallelized in python. -As a second addition, now also Feltor's geometries extension is available in python. +As a second addition, Feltor's geometries extension is available in python. However, the geometries functions and classes are not re-implemented in python, but they are bound to python via the [pybind11](https://github.com/pybind/pybind11) library. As such the corresponding C++ binding code must be compiled in order to generate the module `dg.geo`. ## Installation -### The pyfeltor.dg module -> You need python3 to install this module +> The `pyfeltor.dg.geo` part of the module contains python bindings for the underlying C++ [feltor](https://github.com/feltor-dev/feltor) code using [pybind11](https://github.com/pybind/pybind11). During installation the C++ code will be compiled using cmake, which may take a minute or two. The simplest way is to install from the python package index [pypi](https://pypi.org/) via the package manager [pip](https://pip.pypa.io/en/stable/) v23.0 ```bash @@ -45,34 +44,6 @@ cd tests pytest-3 -s . # run all the unittests with output ``` -### The pyfeltor.dg.geo module -Currently, the only way to install this module is via a local, editable install. -Assuming that the pyfeltor.dg module was succesfully installed this way -- the first step is to also install [feltor](https://github.com/feltor-dev/feltor) -following the quick start guide for a base installation. -- Second, instead of jsoncpp we here use the -[nlohmann/json](https://github.com/nlohmann/json) parser available either as a -system package `nlohmann-json3-dev`. -- Next, we follow the first steps guide on [pybind11](https://github.com/pybind11/pybind11) -and install it via `python3 -m pip install pybind11`. -- Further, we install the 'python3-dev', `pybind11-dev` and `pybind11-json-dev` system packages for -the corresponding C++ header files. - -Finally, invoke the Makefile in this repository -```bash -export FELTOR_PATH=path/to/feltor -make -j 4 -``` -Replace `path/to/feltor` with the path to the Feltor library relative to the current -directory. By default `FELTOR_PATH=../feltor`. - -That's it. With the editable install the `pyfeltor.dg.geo` module is now automatically -imported together with `pyfeltor.dg`. -You can test if it works by executing the test -```bash -cd tests -pytest -s test_geometries -``` ## Usage Generally, pyfeltor is built to mimic the `dg` library in feltor. From 54cb98b3d36cd8ebf05320b950dfaaae0d6387c4 Mon Sep 17 00:00:00 2001 From: Matthias Date: Tue, 24 Jun 2025 14:37:29 +0200 Subject: [PATCH 13/16] Add pypi publishing workflow --- .github/workflows/python_publish.yml | 48 ++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 .github/workflows/python_publish.yml diff --git a/.github/workflows/python_publish.yml b/.github/workflows/python_publish.yml new file mode 100644 index 0000000..286beb5 --- /dev/null +++ b/.github/workflows/python_publish.yml @@ -0,0 +1,48 @@ +name: Build/publish Python Package + +on: push + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.x' + - name: Install dependencies + run: | + python -m pip install --upgrade pip + python -m pip install build twine + - name: Build package + run: | + python -m build --sdist --wheel + - name: Check package + run: | + python -m twine check dist/* + - name: Store the distribution packages + uses: actions/upload-artifact@v4 + with: + name: python-package-distributions + path: dist/ + + publish-to-pypi: + if: startsWith(github.ref, 'refs/tags/') # only publish to PyPI on tag pushes + needs: + - build + runs-on: ubuntu-latest + environment: + name: pypi + permissions: + id-token: write # IMPORTANT: mandatory for trusted publishing + steps: + - name: Download distribution packages + uses: actions/download-artifact@v4 + with: + name: python-package-distributions + path: dist/ + - name: Publish to PyPI + uses: pypa/gh-action-pypi-publish@release/v1 From 835e89ca26bf574e440c1da1736b4c239048f816 Mon Sep 17 00:00:00 2001 From: Matthias Date: Tue, 24 Jun 2025 14:41:03 +0200 Subject: [PATCH 14/16] Add github test workflow --- .github/workflows/test.yml | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 .github/workflows/test.yml diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 0000000..ccb4e8f --- /dev/null +++ b/.github/workflows/test.yml @@ -0,0 +1,25 @@ +name: tests + +on: [push, pull_request] + +jobs: + pytest: + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.10", "3.11", "3.12"] + + steps: + - uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + python -m pip install .[test] + - name: Test with pytest + run: | + pytest -vv tests From 86b6f09078cff192dafbba9429956a597e489260 Mon Sep 17 00:00:00 2001 From: Matthias Date: Tue, 24 Jun 2025 14:43:36 +0200 Subject: [PATCH 15/16] Fix python version in test workflow --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index ccb4e8f..d214768 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -8,7 +8,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.10", "3.11", "3.12"] + python-version: ["3.12", "3.13"] steps: - uses: actions/checkout@v4 From fbb50cd12ab55167219443e49797cb5885a94caf Mon Sep 17 00:00:00 2001 From: Matthias Date: Tue, 24 Jun 2025 15:30:19 +0200 Subject: [PATCH 16/16] Pass all ruff checks --- pyfeltor/dg/__init__.py | 14 ++-- pyfeltor/dg/create/__init__.py | 12 ++-- pyfeltor/dg/create/derivative.py | 3 +- pyfeltor/dg/create/dx.py | 50 ++++++------- pyfeltor/dg/create/elliptic.py | 6 +- pyfeltor/dg/create/interpolation.py | 67 ++++++++++-------- pyfeltor/dg/create/operators.py | 5 +- pyfeltor/dg/create/weights.py | 35 ++++----- pyfeltor/dg/evaluation.py | 16 ++--- pyfeltor/dg/geo/__init__.py | 12 +++- pyfeltor/dg/grid.py | 9 ++- pyproject.toml | 4 ++ tests/test_derivative.py | 25 ++++--- tests/test_elliptic.py | 20 +++--- tests/test_evaluation.py | 3 +- tests/test_geometries.py | 106 ++++++++++++++-------------- tests/test_grid.py | 26 +++---- tests/test_interpolation.py | 7 +- 18 files changed, 223 insertions(+), 197 deletions(-) diff --git a/pyfeltor/dg/__init__.py b/pyfeltor/dg/__init__.py index ba873be..178854b 100644 --- a/pyfeltor/dg/__init__.py +++ b/pyfeltor/dg/__init__.py @@ -1,8 +1,12 @@ # docstring displayed by help(pyfeltor) """ The python version of the dg library """ -from . import create -from .enums import bc, direction, inverse_bc, inverse_dir -from .grid import Grid -from .evaluation import evaluate, integrate -from . import geo +from . import create as create +from . import geo as geo +from .enums import bc as bc +from .enums import direction as direction +from .enums import inverse_bc as inverse_bc +from .enums import inverse_dir as inverse_dir +from .evaluation import evaluate as evaluate +from .evaluation import integrate as integrate +from .grid import Grid as Grid diff --git a/pyfeltor/dg/create/__init__.py b/pyfeltor/dg/create/__init__.py index a7b6c3d..33ec8ed 100644 --- a/pyfeltor/dg/create/__init__.py +++ b/pyfeltor/dg/create/__init__.py @@ -1,7 +1,11 @@ # docstring displayed by help(pyfeltor) """ The python version of the dg library, creation """ -from .weights import weights, abscissas, grid_from_abscissas -from .derivative import dx, jump -from .elliptic import elliptic -from .interpolation import interpolation, projection +from .derivative import dx as dx +from .derivative import jump as jump +from .elliptic import elliptic as elliptic +from .interpolation import interpolation as interpolation +from .interpolation import projection as projection +from .weights import abscissas as abscissas +from .weights import grid_from_abscissas as grid_from_abscissas +from .weights import weights as weights diff --git a/pyfeltor/dg/create/derivative.py b/pyfeltor/dg/create/derivative.py index 6a4fc48..326be1a 100644 --- a/pyfeltor/dg/create/derivative.py +++ b/pyfeltor/dg/create/derivative.py @@ -1,6 +1,7 @@ -from .dx import normed, jump_normed import scipy.sparse +from .dx import jump_normed, normed + def dx(dim, grid, bc, direction): deriv = normed(grid.n[dim], grid.N[dim], grid.h()[dim], bc, direction) diff --git a/pyfeltor/dg/create/dx.py b/pyfeltor/dg/create/dx.py index f4a3965..32496fc 100644 --- a/pyfeltor/dg/create/dx.py +++ b/pyfeltor/dg/create/dx.py @@ -1,12 +1,12 @@ -import numpy as np import scipy.sparse + +from ..enums import bc from . import operators as ops -from ..enums import bc, direction def symm(n, N, h, bcx): - l = ops.lilj(n) - r = ops.rirj(n) + ll = ops.lilj(n) + rr = ops.rirj(n) lr = ops.lirj(n) rl = ops.rilj(n) d = ops.pidxpj(n) @@ -19,14 +19,14 @@ def symm(n, N, h, bcx): a_bound_left = a.copy() # left boundary if (bcx == bc.DIR) or (bcx == bc.DIR_NEU): - a_bound_left += 0.5 * (t @ l) + a_bound_left += 0.5 * (t @ ll) elif (bcx == bc.NEU) or (bcx == bc.NEU_DIR): - a_bound_left -= 0.5 * (t @ l) + a_bound_left -= 0.5 * (t @ ll) # right boundary if (bcx == bc.DIR) or (bcx == bc.NEU_DIR): - a_bound_right -= 0.5 * (t @ r) + a_bound_right -= 0.5 * (t @ rr) elif (bcx == bc.NEU) or (bcx == bc.DIR_NEU): - a_bound_right += 0.5 * (t @ r) + a_bound_right += 0.5 * (t @ rr) b = t @ (1.0 / 2.0 * rl) bp = t @ (-1.0 / 2.0 * lr) # pitfall: T*-m^T is NOT -(T*m)^T # transform to XSPACE @@ -94,20 +94,20 @@ def symm(n, N, h, bcx): vals.append(b[i, j]) # sort - rows, cols, vals = zip(*sorted(zip(rows, cols, vals))) + rows, cols, vals = zip(*sorted(zip(rows, cols, vals, strict=False)), strict=False) return scipy.sparse.coo_matrix((vals, (rows, cols))) def plus(n, N, h, bcx): - l = ops.lilj(n) - r = ops.rirj(n) - lr = ops.lirj(n) + ll = ops.lilj(n) + #rr = ops.rirj(n) + #lr = ops.lirj(n) rl = ops.rilj(n) d = ops.pidxpj(n) t = ops.pipj_inv(n) t *= 2.0 / h - a = t @ (-l - d.transpose()) + a = t @ (-ll - d.transpose()) # bcx = PER a_bound_left = a.copy() # PER, NEU, and NEU_DIR a_bound_right = a.copy() # PER, DIR, and NEU_DIR @@ -168,20 +168,20 @@ def plus(n, N, h, bcx): vals.append(b[i, j]) # sort - rows, cols, vals = zip(*sorted(zip(rows, cols, vals))) + rows, cols, vals = zip(*sorted(zip(rows, cols, vals, strict=False)), strict=False) return scipy.sparse.coo_matrix((vals, (rows, cols))) def minus(n, N, h, bcx): - l = ops.lilj(n) - r = ops.rirj(n) + ll = ops.lilj(n) + #rr = ops.rirj(n) lr = ops.lirj(n) - rl = ops.rilj(n) + #rl = ops.rilj(n) d = ops.pidxpj(n) t = ops.pipj_inv(n) t *= 2.0 / h - a = t @ (l + d) + a = t @ (ll + d) # bcx = PER a_bound_right = a.copy() # PER, NEU and DIR_NEU a_bound_left = a.copy() # PER, DIR and DIR_NEU @@ -242,25 +242,25 @@ def minus(n, N, h, bcx): vals.append(a[i, j]) # sort - rows, cols, vals = zip(*sorted(zip(rows, cols, vals))) + rows, cols, vals = zip(*sorted(zip(rows, cols, vals, strict=False)), strict=False) return scipy.sparse.coo_matrix((vals, (rows, cols))) def jump_normed(n, N, h, bcx): - l = ops.lilj(n) - r = ops.rirj(n) + ll = ops.lilj(n) + rr = ops.rirj(n) lr = ops.lirj(n) rl = ops.rilj(n) t = ops.pipj_inv(n) t *= 2.0 / h - a = t @ (l + r) + a = t @ (ll + rr) a_bound_left = a.copy() # DIR and PER if (bcx == bc.NEU) or (bcx == bc.NEU_DIR): - a_bound_left = t @ r + a_bound_left = t @ rr a_bound_right = a.copy() # DIR and PER if (bcx == bc.NEU) or (bcx == bc.DIR_NEU): - a_bound_right = t @ l + a_bound_right = t @ ll b = -t @ rl bp = -t @ lr # transform to XSPACE @@ -328,7 +328,7 @@ def jump_normed(n, N, h, bcx): vals.append(b[i, j]) # sort - rows, cols, vals = zip(*sorted(zip(rows, cols, vals))) + rows, cols, vals = zip(*sorted(zip(rows, cols, vals, strict=False)), strict=False) return scipy.sparse.coo_matrix((vals, (rows, cols))) diff --git a/pyfeltor/dg/create/elliptic.py b/pyfeltor/dg/create/elliptic.py index 80d8277..cfbbffb 100644 --- a/pyfeltor/dg/create/elliptic.py +++ b/pyfeltor/dg/create/elliptic.py @@ -1,8 +1,8 @@ -import numpy as np -from .derivative import dx, jump -from .. import enums import scipy.sparse as sparse +from .. import enums +from .derivative import dx, jump + def elliptic(grid, bcs, directions, sigma, jumpfactor=1): left = [ diff --git a/pyfeltor/dg/create/interpolation.py b/pyfeltor/dg/create/interpolation.py index 30fbb62..3524b31 100644 --- a/pyfeltor/dg/create/interpolation.py +++ b/pyfeltor/dg/create/interpolation.py @@ -1,10 +1,12 @@ -import numpy as np import itertools + +import numpy as np import scipy.sparse -from .weights import weights -from . import operators as ops -from ..enums import bc, direction + +from ..enums import bc from ..evaluation import evaluate +from . import operators as ops +from .weights import weights def shift(grid, x, bcs): @@ -76,24 +78,25 @@ def interpolation(xs, grid, bcs): else: raise Exception("interpolation not implemented for ndim > 3") for it in tuples: - I = 0 # the index to push back - V = 1 # the value to push back + II = 0 # the index to push back + VV = 1 # the value to push back for kk in range(0, grid.ndim): - I = (I * grid.N[kk] + nn[kk]) * grid.n[kk] + it[kk] - V = V * px[kk][it[kk]] + II = (II * grid.N[kk] + nn[kk]) * grid.n[kk] + it[kk] + VV = VV * px[kk][it[kk]] # print( i, I, V) rows.append(i) - cols.append(round(I)) + cols.append(round(II)) if not negative: - vals.append(V) + vals.append(VV) else: - vals.append(-V) + vals.append(-VV) # sort - rows, cols, vals = zip(*sorted(zip(rows, cols, vals))) - return scipy.sparse.coo_matrix((vals, (rows, cols)), shape = (len(xs[0]),grid.size())) + rows, cols, vals = zip(*sorted(zip(rows, cols, vals, strict=False)), strict=False) + return scipy.sparse.coo_matrix((vals, (rows, cols)), shape=(len(xs[0]), + grid.size())) -def projection( grid_new, grid_old): +def projection(grid_new, grid_old): """ Create a projection between two grids This matrix can be applied to vectors defined on the old (fine) grid to obtain @@ -106,26 +109,28 @@ def projection( grid_new, grid_old): """ ndim = grid_old.ndim if grid_new.ndim != ndim: - raise Exception( "Cannot project between grids with different dimensions") - for i in range( 0, ndim): - if grid_old.N[i] % grid_new.N[i] != 0 : - print( f"WARNING you project between incompatible grids!! old N: {grid_old.N[i]} new N {grid_new.N[i]}") - if grid_old.n[i] < grid_new.n[i] : - print( f"WARNING you project between incompatible grids!! old n: {grid_old.n[i]} new n {grid_new.n[i]}") - wf = scipy.sparse.diags(weights( grid_old)) - points = list() - bcs = [bc.PER for i in range(0,ndim)] + raise Exception("Cannot project between grids with different dimensions") + for i in range(0, ndim): + if grid_old.N[i] % grid_new.N[i] != 0: + print(f"WARNING you project between incompatible grids!! old N: \ +{grid_old.N[i]} new N {grid_new.N[i]}") + if grid_old.n[i] < grid_new.n[i]: + print(f"WARNING you project between incompatible grids!! old n: \ +{grid_old.n[i]} new n {grid_new.n[i]}") + wf = scipy.sparse.diags(weights(grid_old)) + points = [] + bcs = [bc.PER for i in range(0, ndim)] if ndim == 1: - points.append( evaluate( lambda x: x, grid_old)) + points.append(evaluate(lambda x: x, grid_old)) elif ndim == 2: - points.append( evaluate( lambda y,x: y, grid_old)) - points.append( evaluate( lambda y,x: x, grid_old)) + points.append(evaluate(lambda y, x: y, grid_old)) + points.append(evaluate(lambda y, x: x, grid_old)) elif ndim == 3: - points.append( evaluate( lambda z,y,x: z, grid_old)) - points.append( evaluate( lambda z,y,x: y, grid_old)) - points.append( evaluate( lambda z,y,x: x, grid_old)) + points.append(evaluate(lambda z, y, x: z, grid_old)) + points.append(evaluate(lambda z, y, x: y, grid_old)) + points.append(evaluate(lambda z, y, x: x, grid_old)) else: raise Exception("Projection not implemented for ndim > 3") - A = interpolation( points, grid_new, bcs ) - vc = scipy.sparse.diags(1./weights( grid_new)) + A = interpolation(points, grid_new, bcs) + vc = scipy.sparse.diags(1. / weights(grid_new)) return vc @ A.transpose() @ wf diff --git a/pyfeltor/dg/create/operators.py b/pyfeltor/dg/create/operators.py index 5ee4e86..9d7f827 100644 --- a/pyfeltor/dg/create/operators.py +++ b/pyfeltor/dg/create/operators.py @@ -26,9 +26,8 @@ def pidxpj(n): op = np.zeros((n, n)) for i in range(0, n): for j in range(0, n): - if i < j: - if (i + j) % 2 != 0: - op[i, j] = 2 + if (i < j) and ((i + j) % 2 != 0): + op[i, j] = 2 return op diff --git a/pyfeltor/dg/create/weights.py b/pyfeltor/dg/create/weights.py index 2c55621..174eab2 100644 --- a/pyfeltor/dg/create/weights.py +++ b/pyfeltor/dg/create/weights.py @@ -1,4 +1,5 @@ import numpy as np + from ..grid import Grid @@ -25,25 +26,27 @@ def abscissas(grid, dimension=0): abscissas[i * grid.n[dim] + k] = i * h + x0 + h / 2.0 * (1.0 + x[k]) return abscissas -def _get_x0_x1_n_N( x): - for n in range(1,21): - if len(x)%n != 0: + +def _get_x0_x1_n_N(x): + for n in range(1, 21): + if len(x) % n != 0: continue - N = len(x)/n + N = len(x) / n (xx, w) = np.polynomial.legendre.leggauss(n) h = x[n] - x[0] - x0 = x[0] - h/2.0 * (1.0+xx[0]) - first_absc = np.zeros(2*n) - last_value = x0 + h*N - h / 2.0 * (1.0 + xx[0]) - for i in range(0,2): - for k in range(0,n): - first_absc[ i*n+k] = i * h + x0 + h / 2.0 * (1.0 + xx[k]) - if abs( last_value - x[-1]) < 1e-10*abs(x[-1]) + 1e-10: - return (x0, x0+h*N, int(n), int(N)) - raise Exception( "Could not determine grid") + x0 = x[0] - h / 2.0 * (1.0 + xx[0]) + first_absc = np.zeros(2 * n) + last_value = x0 + h * N - h / 2.0 * (1.0 + xx[0]) + for i in range(0, 2): + for k in range(0, n): + first_absc[i * n + k] = i * h + x0 + h / 2.0 * (1.0 + xx[k]) + if abs(last_value - x[-1]) < 1e-10 * abs(x[-1]) + 1e-10: + return (x0, x0 + h * N, int(n), int(N)) + raise Exception("Could not determine grid") return (0, 1, 1, 1) -def grid_from_abscissas( xs ) : + +def grid_from_abscissas(xs): """ This function reverse engineers a Grid that corresponds to the given abscissas. This is useful for example to get a Grid that corresponds to a given simulation output @@ -58,6 +61,6 @@ def grid_from_abscissas( xs ) : x1 = np.zeros(ndim) n = np.arange(ndim) N = np.arange(ndim) - for dim in range(0,ndim): + for dim in range(0, ndim): x0[dim], x1[dim], n[dim], N[dim] = _get_x0_x1_n_N(xs[dim]) - return Grid( x0, x1, n, N) + return Grid(x0, x1, n, N) diff --git a/pyfeltor/dg/evaluation.py b/pyfeltor/dg/evaluation.py index 9bc2f89..c241485 100644 --- a/pyfeltor/dg/evaluation.py +++ b/pyfeltor/dg/evaluation.py @@ -1,7 +1,8 @@ -from .create import operators as ops +import numpy as np + from . import create +from .create import operators as ops from .enums import direction -import numpy as np def evaluate(function, grid): @@ -9,12 +10,11 @@ def evaluate(function, grid): function: has to take numpy arrays as arguments, f(x), f(y,x), f(z,y,x) grid: instance of dg.Grid - return: flat np.array with x the fastest varying dimension. Can be reshaped with reshape(grid.shape) + return: flat np.array with x the fastest varying dimension. Can be reshaped + with reshape(grid.shape) """ - xs = [] ndim = grid.ndim - for dim in range(0, ndim): - xs.append(create.abscissas(grid, dim)) + xs = [create.abscissas(grid, dim) for dim in range(0, ndim)] if ndim == 1: return np.array([function(x) for x in xs[0]]) @@ -59,8 +59,8 @@ def integrate(to_integrate, grid, direction=direction.forward): out = np.zeros(grid.size()) for i in range(0, grid.N[0]): for k in range(0, grid.n[0]): - for l in range(0, grid.n[0]): - out[i * n + k] += ninj[k, l] * to_in[i * n + l] + for ll in range(0, grid.n[0]): + out[i * n + k] += ninj[k, ll] * to_in[i * n + ll] out[i * n + k] += constant for k in range(0, grid.n[0]): constant += h * forward[0, k] * to_in[i * n + k] diff --git a/pyfeltor/dg/geo/__init__.py b/pyfeltor/dg/geo/__init__.py index 6324189..023dd06 100644 --- a/pyfeltor/dg/geo/__init__.py +++ b/pyfeltor/dg/geo/__init__.py @@ -2,9 +2,15 @@ """ The python version of the dg library """ +# Import into their own dg.geo.xxx namespace +from . import circular as circular +from . import guenter as guenter +from . import mod as mod +from . import polynomial as polynomial +from . import solovev as solovev +from . import toroidal as toroidal + # Import into the dg.geo namespace +from .flux import * from .geometries import * from .utility import * -from .flux import * -# Import into their own dg.geo.xxx namespace -from . import polynomial, solovev, guenter, toroidal, circular, mod diff --git a/pyfeltor/dg/grid.py b/pyfeltor/dg/grid.py index a6429ee..b260720 100644 --- a/pyfeltor/dg/grid.py +++ b/pyfeltor/dg/grid.py @@ -30,11 +30,10 @@ def __init__(self, x0, x1, n, N): len(N) except TypeError: self.N = np.asarray([N]) - assert ( # lengths must be consistent - (len(self.x0) == len(self.x1)) - and (len(self.x0) == len(self.n)) - and (len(self.x0) == len(self.N)) - ) + # lengths must be consistent + assert(len(self.x0) == len(self.x1)) + assert(len(self.x0) == len(self.n)) + assert(len(self.x0) == len(self.N)) @property def n(self): diff --git a/pyproject.toml b/pyproject.toml index 137b7c0..d8acbb2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,6 +44,7 @@ test = [ # https://docs.astral.sh/ruff/ # pip install -e .[lint] +# uv run ruff check --fix lint = [ "ruff", # a python linter and formatter ] @@ -69,3 +70,6 @@ select = [ "RUF", # ruff specific "NPY", # NumPy specific ] +[tool.ruff.lint.per-file-ignores] +# Ignore `F403` (import violations) in all `__init__.py` files +"__init__.py" = ["F403"] diff --git a/tests/test_derivative.py b/tests/test_derivative.py index 6df7ad2..c794f29 100644 --- a/tests/test_derivative.py +++ b/tests/test_derivative.py @@ -1,7 +1,7 @@ -import pytest -from pyfeltor import dg import numpy as np +from pyfeltor import dg + # Run with pytest-3 -s . to see stdout output @@ -28,11 +28,11 @@ def test_dx(): gNEU = dg.Grid(np.pi / 2.0, 3 * np.pi / 2.0, n, N) gDIR_NEU = dg.Grid(0, np.pi / 2.0, n, N) gNEU_DIR = dg.Grid(np.pi / 2.0, np.pi, n, N) - g = [gPER, gDIR, gNEU, gDIR_NEU, gNEU_DIR] - bcx = [dg.bc.PER, dg.bc.DIR, dg.bc.NEU, dg.bc.DIR_NEU, dg.bc.NEU_DIR] + grids = [gPER, gDIR, gNEU, gDIR_NEU, gNEU_DIR] + bcxs = [dg.bc.PER, dg.bc.DIR, dg.bc.NEU, dg.bc.DIR_NEU, dg.bc.NEU_DIR] print("TEST NORMAL TOPOLOGY: YOU SHOULD SEE CONVERGENCE FOR ALL OUTPUTS!!!") print("COMPARE TO dx_t.cu") - for (g, bcx) in zip(g, bcx): + for (g, bcx) in zip(grids, bcxs, strict=False): print("Boundary condition ", bcx) hs = dg.create.dx(0, g, bcx, dg.direction.centered) hf = dg.create.dx(0, g, bcx, dg.direction.forward) @@ -45,24 +45,24 @@ def test_dx(): null = dg.evaluate(lambda x: 0, g) error = deri - hs.dot(func) print( - f"Distance to true solution (symmetric): {np.sqrt(np.sum( w1d*error**2) )}" + f"Distance to true solution (symmetric): {np.sqrt(np.sum(w1d * error**2))}" ) error = deri - hf.dot(func) - print(f"Distance to true solution (forward): {np.sqrt(np.sum( w1d*error**2) )}") + print(f"Distance to true solution (forward): {np.sqrt(np.sum(w1d * error**2))}") error = deri - hb.dot(func) print( - f"Distance to true solution (backward): {np.sqrt(np.sum( w1d*error**2) )}" + f"Distance to true solution (backward): {np.sqrt(np.sum(w1d * error**2))}" ) error = null - js.dot(func) print( - f"Distance to true solution (jump ): {np.sqrt(np.sum( w1d*error**2) )}" + f"Distance to true solution (jump ): {np.sqrt(np.sum(w1d * error**2))}" ) def test_derivative(): - n, Nx, Ny, Nz = 3, 24, 28, 100 - print(f"On Grid {n} x {Nx} x {Ny} x {Nz}") - bcx, bcy, bcz = dg.bc.DIR, dg.bc.PER, dg.bc.NEU_DIR + n, Nx, Ny = 3, 24, 28 + print(f"On Grid {n} x {Nx} x {Ny}") + bcx, bcy = dg.bc.DIR, dg.bc.PER g2d = dg.Grid([0.1, 0], [2 * np.pi + 0.1, np.pi], [n, n], [Ny, Nx]) w2d = dg.create.weights(g2d) @@ -76,7 +76,6 @@ def test_derivative(): dy2d = dg.evaluate(cosy, g2d) null2 = dg.evaluate(zero, g2d) sol2 = [dx2d, dy2d, null2, null2] - # binary2[4562611930300281864,4553674328256556132,4567083257206218817,4574111364446550002] sol = [ 0.0010775034079703078, 0.00027314872436790867, diff --git a/tests/test_elliptic.py b/tests/test_elliptic.py index 13b850a..092152e 100644 --- a/tests/test_elliptic.py +++ b/tests/test_elliptic.py @@ -1,18 +1,20 @@ -from pyfeltor import dg import numpy as np -import scipy.sparse.linalg import scipy.linalg +import scipy.sparse.linalg + +from pyfeltor import dg def test_elliptic(): amp = 0.9 - pol = lambda y, x: 1 + amp * np.sin(x) * np.sin(y) - rhs = ( - lambda y, x: 2.0 * np.sin(x) * np.sin(y) * (amp * np.sin(x) * np.sin(y) + 1) - - amp * np.sin(x) * np.sin(x) * np.cos(y) * np.cos(y) - - amp * np.cos(x) * np.cos(x) * np.sin(y) * np.sin(y) - ) - sol = lambda y, x: np.sin(x) * np.sin(y) + def pol(y, x): + return 1 + amp * np.sin(x) * np.sin(y) + def rhs(y, x): + return 2.0 * np.sin(x) * np.sin(y) * (amp * np.sin(x) * np.sin(y) + 1)\ + - amp * np.sin(x) * np.sin(x) * np.cos(y) * np.cos(y)\ + - amp * np.cos(x) * np.cos(x) * np.sin(y) * np.sin(y) + def sol(y,x): + return np.sin(x) * np.sin(y) lx, ly = np.pi, 2 * np.pi bcx, bcy = dg.bc.DIR, dg.bc.PER n, Nx, Ny = 3, 64, 64 diff --git a/tests/test_evaluation.py b/tests/test_evaluation.py index 86e6a88..30604e7 100644 --- a/tests/test_evaluation.py +++ b/tests/test_evaluation.py @@ -1,5 +1,5 @@ -import pytest import numpy as np + from pyfeltor import dg # Run with pytest-3 -s . to see stdout output @@ -36,7 +36,6 @@ def test_integration(): n = 3 Nx = 12 Ny = 28 - Nz = 100 g1d = dg.Grid(x0=1, x1=2, n=n, N=Nx) g2d = dg.Grid((0, 0), (2 * np.pi, 2 * np.pi), (n, n), (Ny, Nx)) w1d = dg.create.weights(g1d) diff --git a/tests/test_geometries.py b/tests/test_geometries.py index d3ca102..1e8f51f 100644 --- a/tests/test_geometries.py +++ b/tests/test_geometries.py @@ -1,87 +1,85 @@ -import pytest import json + import numpy as np + from pyfeltor import dg -import pytest_datadir def test_polynomial(): - c = np.array( [1,2,3,4]) - params = {"R_0" : 400, "inverseaspectratio" : 20, "elongation" : 1, "triangularity" : 1, - "PP" : 1, "PI" : 1, "description" : "standardX", "M" : 2, "N" : 2, "c" : c.tolist()} + c = np.array([1, 2, 3, 4]) + params = {"R_0": 400, "inverseaspectratio": 20, "elongation": 1, "triangularity": 1, + "PP": 1, "PI": 1, "description": "standardX", "M": 2, "N": 2, + "c": c.tolist()} pp = dg.geo.polynomial.Parameters(params) - psip = dg.geo.polynomial.Psip( pp) - grid = dg.Grid( x0 = (pp.R_0-pp.a, -pp.a), x1 = (pp.R_0+pp.a, pp.a), n=(3,3), N=(24, 24)) - psi = dg.evaluate( psip, grid) - R = dg.evaluate( lambda R, Z: R, grid) - Z = dg.evaluate( lambda R, Z: Z, grid) - psi_ = psip(R,Z) - print( psi) - print( psi_) + psip = dg.geo.polynomial.Psip(pp) + grid = dg.Grid(x0=(pp.R_0 - pp.a, -pp.a), x1=(pp.R_0 + pp.a, pp.a), + n=(3, 3), N=(24, 24)) + psi = dg.evaluate(psip, grid) + R = dg.evaluate(lambda R, Z: R, grid) + Z = dg.evaluate(lambda R, Z: Z, grid) + psi_ = psip(R, Z) + print(psi) + print(psi_) def test_make_field(datadir): - with open (datadir / "geometry_params_Xpoint.json", "r") as f: + with open(datadir / "geometry_params_Xpoint.json") as f: magparams = json.load(f) - mag = dg.geo.createMagneticField( magparams) + mag = dg.geo.createMagneticField(magparams) a = mag.params().a() R0 = mag.R0() - grid = dg.Grid(x0=(R0-a, -a), x1=(R0+a, +a), n=(3, 3), N=(24, 24)) - psi = dg.evaluate( mag.psip(), grid) - BR = dg.evaluate( dg.geo.BFieldR(mag), grid) - print( mag.R0()) - print( mag.params().a()) - print( mag.params().elongation()) - print( mag.params().getDescription()) - print( mag.params().getEquilibrium()) - print (psi, BR) - + grid = dg.Grid(x0=(R0 - a, -a), x1=(R0 + a, +a), n=(3, 3), N=(24, 24)) + psi = dg.evaluate(mag.psip(), grid) + BR = dg.evaluate(dg.geo.BFieldR(mag), grid) + print(mag.R0()) + print(mag.params().a()) + print(mag.params().elongation()) + print(mag.params().getDescription()) + print(mag.params().getEquilibrium()) + print(psi, BR) def test_q_profile(datadir): - with open (datadir / "enrx_tcv.json", "r") as f: + with open(datadir / "enrx_tcv.json") as f: magparams = json.load(f) mag = dg.geo.createMagneticField(magparams) qfunctor = dg.geo.SafetyFactor(mag) - RO,ZO = mag.R0(), 0 + RO, ZO = mag.R0(), 0 (point, RO, ZO) = dg.geo.findOpoint(mag.get_psip(), RO, ZO) - print( "O-point found at ", RO, ZO) - psipO = mag.psip()(RO,ZO) - psi_values = np.linspace( psipO, 0, 20, endpoint = False) + print("O-point found at ", RO, ZO) + psipO = mag.psip()(RO, ZO) + psi_values = np.linspace(psipO, 0, 20, endpoint=False) print(qfunctor(psi_values)) + def test_sheath(datadir): - with open (datadir / "enrx_tcv.json", "r") as f: + with open(datadir / "enrx_tcv.json") as f: magparams = json.load(f) mag = dg.geo.createMagneticField(magparams) - wall = {"type" : "sol_pfr", "alpha": [0.05,0.05], "boundary" : [1.09,0.97]} - sheath= {"boundary" : 3/32, "alpha" : 2/32} + wall = {"type": "sol_pfr", "alpha": [0.05, 0.05], "boundary": [1.09, 0.97]} + sheath = {"boundary": 3 / 32, "alpha": 2 / 32} wall_f = dg.geo.CylindricalFunctor() transition_f = dg.geo.CylindricalFunctor() sheath_f = dg.geo.CylindricalFunctor() - mod_mag = dg.geo.createModifiedField( magparams, wall, wall_f, transition_f) + dg.geo.createModifiedField(magparams, wall, wall_f, transition_f) a = mag.params().a() R0 = mag.R0() - print(R0-a, -a, R0+a, +a) - grid = dg.Grid(x0=(R0-a, -a), x1=(R0+a, +a), n=(3, 3), N=(24, 24)) - dg.geo.createSheathRegion( sheath, mag, wall_f, grid, sheath_f) - - point,RO,ZO = dg.geo.findOpoint( mag.get_psip(), R0, 0) - assert( point == 1) - assert( abs(RO - 0.9022324795511738) < 1e-10) - assert( abs(ZO - (-0.005684720497897021)) < 1e-10) - - closed = dg.geo.mod.ClosedFieldlineRegion( mag); - assert( closed( 0, R0)) - opened = dg.geo.mod.ClosedFieldlineRegion( mag, False); - assert( not opened( 0, R0)) - sol = dg.geo.mod.SOLRegion( mag, wall_f); - assert( not sol( 0, R0)) - - cutter = dg.geo.ZCutter( 10) - assert( cutter( 1, 11)) - assert( not cutter( 1, 8)) - + print(R0 - a, -a, R0 + a, +a) + grid = dg.Grid(x0=(R0 - a, -a), x1=(R0 + a, +a), n=(3, 3), N=(24, 24)) + dg.geo.createSheathRegion(sheath, mag, wall_f, grid, sheath_f) + point, RO, ZO = dg.geo.findOpoint(mag.get_psip(), R0, 0) + assert (point == 1) + assert (abs(RO - 0.9022324795511738) < 1e-10) + assert (abs(ZO - (-0.005684720497897021)) < 1e-10) + closed = dg.geo.mod.ClosedFieldlineRegion(mag) + assert (closed(0, R0)) + opened = dg.geo.mod.ClosedFieldlineRegion(mag, False) + assert (not opened(0, R0)) + sol = dg.geo.mod.SOLRegion(mag, wall_f) + assert (not sol(0, R0)) + cutter = dg.geo.ZCutter(10) + assert (cutter(1, 11)) + assert (not cutter(1, 8)) diff --git a/tests/test_grid.py b/tests/test_grid.py index c265910..c63fc0e 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -1,5 +1,5 @@ -import pytest import numpy as np + from pyfeltor import dg # Run with pytest-3 -s . to see stdout output @@ -28,17 +28,19 @@ def test_evaluation(): grid = dg.Grid(x0=(0, 0), x1=(np.pi, np.pi), n=(3, 3), N=(24, 24)) weights = dg.create.weights(grid) sine = dg.evaluate(lambda y, x: np.sin(y) * np.sin(x), grid) - print(f"Integral is {np.sum( sine*weights)}") + print(f"Integral is {np.sum(sine * weights)}") + def test_grid_from_abscissas(): print("TEST REVERSE GRID ENGINEER") - grid = dg.Grid(x0=(-10, 0, 20,-10,0,20), x1=(2*np.pi, np.pi, 10*np.pi,2*np.pi, np.pi, 10*np.pi), n=(1, 2, 3,4,5,6), N=(20, 24, 12,13,15,16)) - xs = [dg.create.abscissas(grid,i) for i in range(0,grid.ndim)] - test_grid = dg.create.grid_from_abscissas( xs) - - print( grid.x0, grid.x1, grid.n, grid.N) - print( test_grid.x0, test_grid.x1, test_grid.n, test_grid.N) - assert np.allclose( grid.x0, test_grid.x0) - assert np.allclose( grid.x1, test_grid.x1) - assert np.allclose( grid.n, test_grid.n) - assert np.allclose( grid.N, test_grid.N) + grid = dg.Grid(x0=(-10, 0, 20, -10, 0, 20), x1=(2 * np.pi, np.pi, 10 * np.pi, + 2 * np.pi, np.pi, 10 * np.pi), n=(1, 2, 3, 4, 5, 6), N=(20, 24, 12, 13, 15, 16)) + xs = [dg.create.abscissas(grid, i) for i in range(0, grid.ndim)] + test_grid = dg.create.grid_from_abscissas(xs) + + print(grid.x0, grid.x1, grid.n, grid.N) + print(test_grid.x0, test_grid.x1, test_grid.n, test_grid.N) + assert np.allclose(grid.x0, test_grid.x0) + assert np.allclose(grid.x1, test_grid.x1) + assert np.allclose(grid.n, test_grid.n) + assert np.allclose(grid.N, test_grid.N) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 6b59abb..a4a350e 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -1,6 +1,7 @@ -from pyfeltor import dg import numpy as np +from pyfeltor import dg + def test_interpolation1d(): print("1D INTERPOLATION TEST") @@ -76,9 +77,9 @@ def test_projection(): n, Nx, Ny = 3, 8, 8 g = dg.Grid(x0=[-5 * np.pi, -np.pi], x1=[-4 * np.pi, 0], n=[n, n], N=[Ny, Nx]) - g_fine = dg.Grid( x0=g.x0, x1=g.x1, n=g.n, N=[n * Ny, n * Nx]) + g_fine = dg.Grid(x0=g.x0, x1=g.x1, n=g.n, N=[n * Ny, n * Nx]) - project = dg.create.projection( g, g_fine) + project = dg.create.projection(g, g_fine) vec = dg.evaluate(lambda y, x: np.sin(x) * np.sin(y), g_fine) project = project.dot(vec) projectE = dg.evaluate(lambda y, x: np.sin(x) * np.sin(y), g)