2.3.1 joss

log-grid/CHANGELOG.md

@@ -5,6 +5,33 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres
 to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [2.3.1]
+
+### Fixed
+
+- Fix build process
+
+## [2.3.0]
+
+### Changed
+
+- Improve docs
+- Improve `Utils/benchmark_c` scripts and visualization
+- Move `Maths` and `Physics` subclasses outside of the `Grid` class
+- Move `ETD4RK` from `Examples` to `Utils`
+
+### Fixed
+
+- Error on a,b not integers in l_params
+- Fix a potential error where `plastic` was not supplied in `l_params`
+- Fix compiler warnings on `convolver.c`
+- Fix signature of interactive functions in `plotLib`
+- Improve stability of windows installation
+
+### Removed
+
+- Remove obsolete `Examples/2D_basilisk/`
+
 ## [2.2.1]
 
 ### Changed
@@ -15,6 +42,7 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 - Fix Gitlab CI to only publish on new master tags
 - Fix docs
+- Fix `poetry install` when compiled files already exist
 
 ## [2.2.0]
 

log-grid/Utils/ETD4RK/README.md

@@ -0,0 +1,9 @@
+# ETD4RK Thresholds
+
+ETD4RK is a solver based on _Exponential time differencing for stiff systems_ (Cox & Matthews 2002), equations (26)-(29).
+
+In its equations intervene differences of functions and their truncated Taylor series, such as $e^x-1-x$.
+
+Due to the floating point precision of python, for small enough $x$, this difference becomes numerically unstable, and should be replaced by the higher order term of the Taylor series.
+
+This folder contains the `thresholds.py` script which showcases both the exact numerical function and the higher-order taylor term, as a means to determine at which threshold we should switch from one to the other. Those thresholds are the ones used in `Framework.ETD4RK`.

log-grid/Utils/ETD4RK/thresholds.py

@@ -0,0 +1,57 @@
+"""see README.md: empirical determination of thresholds between numerical differences and higher-order Taylor terms in ETD4RK"""
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+# res = (e - 1) / (hc/h)
+f = lambda x: (np.exp(x / 2) - 1) / x
+g = lambda x: np.ones_like(x) * 1 / 2
+
+x = np.logspace(-10, 0, int(1e4))
+plt.loglog(x, f(x), label="numerical difference")
+plt.loglog(x, g(x), "--", label="higher-order taylor")
+plt.loglog(x, f(x) - g(x), ":", label="diff")
+plt.axvline(x=1e-6, linestyle="--", label="threshold", color="black")
+plt.xlabel("$x$")
+plt.title("$(e^{x/2}-1)/x$")
+plt.legend()
+plt.show()
+
+# res = (-4 - hc + e * (4 - 3 * hc + hc ** 2)) / ((hc)**3/h)
+f = lambda x: (-4 - x + np.exp(x) * (4 - 3 * x + x**2)) / x**3
+g = lambda x: np.ones_like(x) * 1 / 6
+
+plt.loglog(x, f(x), label="numerical difference")
+plt.loglog(x, g(x), "--", label="higher-order taylor")
+plt.loglog(x, np.abs(f(x) - g(x)), ":", label="diff")
+plt.axvline(x=1e-3, linestyle="--", label="threshold", color="black")
+plt.xlabel("$x$")
+plt.title("$(-4 - x + e^x (4 - 3x + x^2)) / x^3$")
+plt.legend()
+plt.show()
+
+# res = (2 + hc + e * (-2 + hc)) / ((hc)**3/h)
+f = lambda x: (2 + x + np.exp(x) * (-2 + x)) / x**3
+g = lambda x: np.ones_like(x) * 1 / 6
+
+plt.loglog(x, f(x), label="numerical difference")
+plt.loglog(x, g(x), "--", label="higher-order taylor")
+plt.loglog(x, np.abs(f(x) - g(x)), ":", label="diff")
+plt.axvline(x=1e-3, linestyle="--", label="threshold", color="black")
+plt.xlabel("$x$")
+plt.title("$(2 + x + e^x (-2+x)) / x^3$")
+plt.legend()
+plt.show()
+
+# res = (-4 - 3 * hc - hc ** 2 + e * (4 - hc)) / ((hc)**3/h)
+f = lambda x: (-4 - 3 * x - x**2 + np.exp(x) * (4 - x)) / x**3
+g = lambda x: np.ones_like(x) * 1 / 6
+
+plt.loglog(x, f(x), label="numerical difference")
+plt.loglog(x, g(x), "--", label="higher-order taylor")
+plt.loglog(x, np.abs(f(x) - g(x)), ":", label="diff")
+plt.axvline(x=6e-3, linestyle="--", label="threshold", color="black")
+plt.xlabel("$x$")
+plt.title("$(-4-3 x - x^2 + e^x (4 - x)) / x^3$")
+plt.legend()
+plt.show()

log-grid/Utils/benchmark_c/analyze_results.py

@@ -35,21 +35,26 @@ def plot_results():
     fig, axs = plt.subplots(N_ax, 1)
     if N_ax == 1:
         axs = [axs]
-    colors = ["blue", "red", "green", "orange"]
+    colors = ["blue", "red", "green", "orange", "teal", "yellow"]
     for env_i, (env, res_env) in enumerate(results.items()):
         ax = axs[env_i]
         for i, (flags, data) in enumerate(res_env.items()):
             ax.plot(np.ones_like(data) * i, data, alpha=0.005, linestyle="", marker="o", markeredgewidth=0, color=colors[env_i])
             ax.semilogy(i, np.mean(data), linestyle="", marker="o", color="black")
             ax.text(i - 0.3, np.mean(data), flags, rotation="vertical")
         ax.plot([], marker="o", color=colors[env_i], label=env or "default")
+        if env_i == 0:
+            ax.set_ylabel("execution time")
+        ax.set_xticks([])
+        ax.set_xlim(-0.5, None)
         ax.legend()
 
     if default_comp in results:
         ax = axs[-1]
         default_i = None
         for i, flags in enumerate(results[default_comp].keys()):
-            if flags.startswith("'DEFAULT"):
+            is_default = flags.startswith("'DEFAULT")
+            if is_default:
                 if default_i is None:
                     default_i = i
                 i = default_i
@@ -58,12 +63,19 @@ def plot_results():
                     data = res_env[flags]
                     ax.semilogy(i, np.mean(data), linestyle="", marker="o", color=colors[env_i])
             text_pos = np.mean(results[default_comp][flags])
-            if flags.startswith("'DEFAULT"):
-                flags = "DEFAULT"
+            if is_default:
+                if flags.endswith("0'"):
+                    flags = "DEFAULT"
+                else:
+                    continue
             ax.text(i - 0.3, text_pos, flags, rotation="vertical")
+        ax.set_xticks([])
+        ax.set_title("Aggregated averages")
+        ax.set_xlim(-0.5, None)
 
     # plt.tight_layout()
     plt.show()
 
 
-plot_results()
+if __name__ == "__main__":
+    plot_results()

log-grid/Utils/benchmark_c/c_optimizer.py

@@ -1,10 +1,32 @@
-"""find fastest compiler settings"""
+"""find fastest compiler settings by brute force
+
+``flags_to_test``: a set of flags to test
+``flags_default``: default flags used before testing flags
+``N_default``: how many times to run only the default flags
+``compilers``: set of compilers to use
+``N_cycles``: how many convolutions to run per benchmark
+``N``: which grid size to choose
+``plot_results``: if True, show results
+
+This will test for each compiler:
+* ``N_default`` times the default flags
+* Once ``flags_default`` + each member of ``flags_to_test`` separately
+
+Individual benchmark configurations are in ``run_benchmark.py``.
+"""
+import re
 import shutil
 from subprocess import call
 
-flags_all = {"-ftree-vectorize", "-maxv"}
-# flags_default = ["-c", "-g", "-fpic", "-Wall", "-Wextra", "-O3", "-ffast-math", "-march=native", "-fopenmp", "-DUSE_OMP"]
-flags_default = []
+from Utils.benchmark_c.analyze_results import plot_results
+
+flags_to_test = {"-ftree-vectorize", "-maxv", "-march=native"}
+flags_default = ["-c", "-g", "-fpic", "-Wall", "-Wextra", "-O3", "-ffast-math", "-fopenmp", "-DUSE_OMP"]
+N_default = 5  # Number of times to run just the default flags (to assess variability)*
+compilers = {"clang-15", "gcc-11"}
+N_cycles = 1e3
+N = 13
+show_results = True
 
 
 def benchmark_config(flags: list[str], comp="clang") -> None:
@@ -18,50 +40,30 @@ def benchmark_config(flags: list[str], comp="clang") -> None:
     def do_compile():
         """compile the module"""
         shutil.copyfile("../../pyloggrid/LogGrid/convolver_c.c", "convolver_c.c")
+        with open("../../pyloggrid/LogGrid/Makefile", "r", encoding="utf8") as f:
+            makefile_src = f.read()
+
+        makefile_src = re.sub(r"CC\s*=.*\n", f"CC={comp}\n", makefile_src)
+        makefile_src = re.sub(r"\n\s*python setup_.*\.py build_ext --inplace", "", makefile_src)
+        makefile_src = re.sub(r"CFLAGS\s?=.*\n", f"CFLAGS={' '.join(flags_c)}\n", makefile_src)
+
+        with open("makefile", "w", encoding="utf8") as f:
+            f.write(makefile_src)
 
-        setup_str = f"""
-CC={comp}
-CFLAGS=-c -g -fpic -Wall -Wextra -O3 -ffast-math -march=native -fopenmp -DUSE_OMP {" ".join(flags_c)}
-LDFLAGS=-shared -fopenmp
-SOURCES=convolver_c.c
-OBJECTS=$(SOURCES:.c=.o)
-TARGET=convolver_c.so
-
-.PHONY: all clean
-
-all: $(TARGET)
-
-$(TARGET): $(OBJECTS)
-    $(CC) $(LDFLAGS) $(OBJECTS) -o $@
-    rm -f $(OBJECTS) convolver.html convolver.c
-    rm -rf build
-
-.c.o:
-    $(CC) $(CFLAGS) $< -o $@
-
-clean:
-    rm -f $(OBJECTS) $(TARGET) convolver.html convolver.c convolver.*.so
-    rm -rf build
-""".replace(
-            "    ", "	"
-        )
-        with open("makefile", "w") as f:
-            f.write(setup_str)
         call(["make"])
 
     do_compile()
     print(f"compiled {flags}")
-    call(["python", "run_benchmark.py", str(flags) + str(comp)])
+    call(["python", "run_benchmark.py", str(flags) + str(comp), str(N_cycles), str(N)])
 
 
 def run_benchmarks():
     """run the benchmarks"""
-    flags_totest = [f for f in flags_all if f not in flags_default]
+    flags_to_test_filtered = [f for f in flags_to_test if f not in flags_default]
 
     errors = []
-    N_default = 5
-    for env in ["clang-15"]:  # ["gcc", "gcc-12", "clang-15"]:
-        for flags in [[f"DEFAULT_{i}"] for i in range(N_default)] + [[flag] for flag in flags_totest]:
+    for env in compilers:  # ["gcc", "gcc-12", "clang-15"]:
+        for flags in [[f"DEFAULT_{i}"] for i in range(N_default)] + [[flag] for flag in flags_to_test_filtered]:
             try:
                 benchmark_config(flags, env)
             except Exception as e:
@@ -72,4 +74,7 @@ def run_benchmarks():
         print("errors:\n" + "\n".join([f"{f}: {e}" for f, e in errors]))
 
 
-run_benchmarks()
+if __name__ == "__main__":
+    run_benchmarks()
+    if show_results:
+        plot_results()

log-grid/Utils/benchmark_c/run_benchmark.py

@@ -1,33 +1,34 @@
-"""run a benchmark for current build"""
+"""
+run a benchmark for current build
+
+Ran as ``run_benchmark.py flags N_cycles N``
+"""
 import sys
 import time
+from pathlib import Path
 
 import numpy as np
 import numpy.ctypeslib as npct
 from matplotlib import pyplot as plt
 
-convolver_c = npct.load_library("convolver_c", ".")  # ../../Utils/
-import ctypes as ct
-
-
-def define_arguments():
-    """Convenience function for defining the arguments of the functions inside the imported module."""
-
-    # Define the arguments accepted by the C functions. This is not strictly necessary,
-    # but it is good practice for avoiding segmentation faults.
-    npflags = ["C_CONTIGUOUS"]  # Require a C contiguous array in memory
-    ui32_1d_type = npct.ndpointer(dtype=np.uint32, ndim=1, flags=npflags)
-    complex_1d_type = npct.ndpointer(dtype=complex, ndim=1, flags=npflags)
+# noinspection PyProtectedMember
+from pyloggrid.LogGrid.Grid import _setup_convolver_c
 
-    args = [ui32_1d_type, ct.c_uint32, complex_1d_type, complex_1d_type, complex_1d_type]
-    convolver_c.convolve.argtypes = args
+convolver_c = npct.load_library("convolver_c", ".")
 
+_setup_convolver_c(convolver_c)
 
-define_arguments()
+N_cycles = float(sys.argv[2])
+N = int(sys.argv[3])
 
 
-def do_benchmark(N=13, N_cycles=1e3, show=True):
-    """benchmark"""
+def do_benchmark(N: int = 13, N_cycles: float = 1e4, show: bool = True) -> None:
+    """do a benchmark
+    Args:
+        N: the grid size
+        N_cycles: the number of convolutions to do
+        show: if True, plot a hist of results
+    """
     from pyloggrid.LogGrid.Grid import Grid
 
     def do_cycle(params: dict) -> None:
@@ -71,7 +72,9 @@ def init_newparam() -> dict:
     flags = sys.argv[1]
     if flags.startswith("[]"):  # default, can be run several timesq
         flags = flags.replace("[]", f"[DEFAULT_{np.random.randint(0, 9999)}]")
+    Path("results").mkdir(exist_ok=True)
     np.save(f"results/flags={flags}", cycles_times)
 
 
-do_benchmark(show=False, N_cycles=1e4, N=15)
+if __name__ == "__main__":
+    do_benchmark(show=False, N_cycles=N_cycles, N=N)

log-grid/Utils/install_clang15.sh

@@ -14,7 +14,7 @@ sudo cp ./lib/libomp.so /usr/lib
 sudo update-alternatives --install /usr/bin/clang-15 clang-15 ~/llvm-project-15.0.7.src/build/bin/clang 100 --slave /usr/bin/clang++-15 clang++-15 ~/llvm-project-15.0.7.src/build/bin/clang++
 sudo update-alternatives --install /usr/bin/clang clang ~/llvm-project-15.0.7.src/build/bin/clang 100
 
-# # If Ubuntu<20: Setup CMake because the system one is too old
+# # If Ubuntu<20: Setup CMake because the system's default is too old
 #wget https://github.com/Kitware/CMake/releases/download/v3.23.4/cmake-3.23.4-linux-x86_64.tar.gz
 #echo "3fbcbff85043d63a8a83c8bdf8bd5b1b2fd5768f922de7dc4443de7805a2670d  cmake-3.23.4-linux-x86_64.tar.gz" | sha256sum -c
 #tar -xf cmake-3.23.4-linux-x86_64.tar.gz