Merge pull request #66 from n-claes/feature/matrix_rework

Rework matrix assembly + Hall/viscosity/v01/B01 extensions

.github/workflows/regression_tests.yml

@@ -66,11 +66,12 @@ jobs:
     - name: Run regression tests
       run: |
         cd tests/regression_tests
-        pytest regression.py test* -v --mpl --mpl-results-path=results
+        pytest regression.py -v
 
     - name: Archive failed logs
-      uses: actions/upload-artifact@v1
+      uses: actions/upload-artifact@v2
       if: failure()
       with:
         name: failed_logs
-        path: tests/regression_tests/results
+        path: tests/regression_tests/image_comparisons
+        if-no-files-found: warn

legolas_config.par

@@ -7,14 +7,27 @@
   mesh_accumulation = .false.
 /
 
+&solvelist
+  solver = "QR-invert"
+  ! solver = "arnoldi"
+  ! arpack_mode = "shift-invert"
+  ! number_of_eigenvalues = 10
+  ! which_eigenvalues = "LI"
+  ! sigma = (0.0d0, 0.01d0)
+/
+
+&physicslist
+  ! hall_mhd = .true.
+/
+
 &equilibriumlist
   ! equilibrium_type = "adiabatic_homo"               ! Cartesian
   ! equilibrium_type = "resistive_homo"               ! Cartesian
   ! equilibrium_type = "gravito_acoustic"             ! Cartesian
   ! equilibrium_type = "gravito_mhd"                  ! Cartesian
-  ! equilibrium_type = "resistive_tearing"            ! Cartesian
+  equilibrium_type = "resistive_tearing"            ! Cartesian
   ! equilibrium_type = "resistive_tearing_flow"       ! Cartesian
-  equilibrium_type = "suydam_cluster"               ! cylindrical
+  ! equilibrium_type = "suydam_cluster"               ! cylindrical
   ! equilibrium_type = "rotating_plasma_cylinder"     ! cylindrical
   ! equilibrium_type = "kelvin_helmholtz_cd"          ! cylindrical
   ! equilibrium_type = "internal_kink"                ! cylindrical
@@ -32,6 +45,7 @@
   ! equilibrium_type = "kelvin_helmholtz"
   ! equilibrium_type = "RTI_KHI"
   ! equilibrium_type = "isothermal_atmosphere"
+  ! equilibrium_type = "harris_sheet"
   ! equilibrium_type = "user_defined"
 
   ! boundary_type = "wall"

post_processing/pylbo/_version.py

@@ -2,7 +2,7 @@
 from packaging import version
 
 # The current pylbo version number.
-__version__ = "1.0.1"
+__version__ = "1.1.0"
 
 
 class VersionHandler:

post_processing/pylbo/automation/defaults.py

@@ -28,6 +28,7 @@
         ("basename_logfile", str),
         ("output_folder", str),
         ("logging_level", (int, np.integer)),
+        ("dry_run", bool),
     ],
     "physicslist": [
         ("mhd_gamma", float),
@@ -47,6 +48,13 @@
         ("use_eta_dropoff", bool),
         ("dropoff_edge_dist", (int, np.integer, float)),
         ("dropoff_width", (int, np.integer, float)),
+        ("viscosity", bool),
+        ("viscosity_value", (int, np.integer, float)),
+        ("hall_mhd", bool),
+        ("hall_dropoff", bool),
+        ("elec_pressure", bool),
+        ("elec_inertia", bool),
+        ("inertia_dropoff", bool),
     ],
     "unitslist": [
         ("cgs_units", bool),

post_processing/pylbo/utilities/datfile_utils.py

@@ -407,8 +407,6 @@ def read_eigenfunction(istream, header, ef_index):
         reals = hdr[::2]
         imags = hdr[1::2]
         ef_values = np.asarray([complex(x, y) for x, y in zip(reals, imags)])
-        # omit very small values
-        ef_values[np.where(np.isclose(ef_values, 0, atol=1e-12))] = 0.0
         eigenfunctions.update({name: ef_values})
     return eigenfunctions
 

post_processing/pylbo/utilities/debug.py

@@ -0,0 +1,136 @@
+import struct
+import numpy as np
+import matplotlib.pyplot as plt
+from pylbo.utilities.datfile_utils import ALIGN
+
+
+def get_debug_coolingcurves(filename):
+    coolcurves = {}
+    with open(filename, "rb") as istream:
+        istream.seek(0)
+        fmt = ALIGN + "i"
+        (size_tables,) = struct.unpack(fmt, istream.read(struct.calcsize(fmt)))
+        fmt = ALIGN + size_tables * "d"
+        coolcurves["T_table"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        coolcurves["L_table"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        fmt = ALIGN + "i"
+        (size_curves,) = struct.unpack(fmt, istream.read(struct.calcsize(fmt)))
+        fmt = ALIGN + size_curves * "d"
+        coolcurves["T_curve"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        coolcurves["L_curve"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        coolcurves["dLdT_curve"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+    return coolcurves
+
+
+def get_debug_atmocurves(filename):
+    atmocurves = {}
+    with open(filename, "rb") as istream:
+        istream.seek(0)
+        fmt = ALIGN + "i"
+        (size_tables,) = struct.unpack(fmt, istream.read(struct.calcsize(fmt)))
+        fmt = ALIGN + size_tables * "d"
+        atmocurves["h_table"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        atmocurves["T_table"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        atmocurves["nh_table"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        fmt = ALIGN + "i"
+        (size_curves,) = struct.unpack(fmt, istream.read(struct.calcsize(fmt)))
+        fmt = ALIGN + size_curves * "d"
+        atmocurves["h_curve"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        atmocurves["T_curve"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        atmocurves["nh_curve"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+        atmocurves["dTdr_curve"] = np.asarray(
+            struct.unpack(fmt, istream.read(struct.calcsize(fmt))),
+            dtype=float,
+        )
+    return atmocurves
+
+
+def plot_debug_coolingcurves(filename):
+    coolcurves = get_debug_coolingcurves(filename)
+
+    fig, (axl, axr) = plt.subplots(1, 2, figsize=(10, 6))
+    axl.loglog(coolcurves["T_table"], coolcurves["L_table"], ".b", label="table")
+    axl.loglog(coolcurves["T_curve"], coolcurves["L_curve"], "-r", label="curve")
+    axl.set_xlabel("temperature [K]")
+    axl.set_ylabel(r"$\Lambda(T)$ [erg s$^{-1}$ cm$^3$]")
+    axl.legend(loc="best")
+    axr.semilogx(
+        coolcurves["T_curve"],
+        coolcurves["dLdT_curve"],
+        color="red",
+        label="legolas 6th O",
+    )
+    axr.semilogx(
+        coolcurves["T_curve"],
+        np.gradient(coolcurves["L_curve"], coolcurves["T_curve"], edge_order=2),
+        label="numpy 2nd O",
+    )
+    axr.set_xlabel("temperature [K]")
+    axr.set_ylabel(r"$\frac{d\Lambda(T)}{dT}$")
+    axr.legend(loc="best")
+    fig.suptitle(f"cooling curves, interpolated at {len(coolcurves['L_curve'])} points")
+    fig.tight_layout()
+    return (fig, (axl, axr))
+
+
+def plot_debug_atmocurves(filename):
+    atmocurves = get_debug_atmocurves(filename)
+
+    fig, (axl, axr) = plt.subplots(1, 2, figsize=(10, 6))
+    axl.semilogy(atmocurves["h_table"], atmocurves["T_table"], ".b", label="table")
+    axl.semilogy(atmocurves["h_curve"], atmocurves["T_curve"], "-r", label="curve")
+    axl.set_xlabel("height [cm]")
+    axl.set_ylabel("temperature [K]")
+    axl.legend(loc="best")
+    axr.plot(
+        atmocurves["h_curve"],
+        atmocurves["dTdr_curve"],
+        color="red",
+        label="legolas 6th O",
+    )
+    axr.plot(
+        atmocurves["h_curve"],
+        np.gradient(atmocurves["T_curve"], atmocurves["h_curve"], edge_order=2),
+        label="numpy 2nd O",
+    )
+    axr.set_xlabel("height [cm]")
+    axr.set_ylabel(r"$\frac{dT}{dr}$")
+    axr.legend(loc="best")
+    fig.suptitle(
+        f"atmosphere curves, interpolated at {len(atmocurves['T_curve'])} points"
+    )
+    fig.tight_layout()
+    return (fig, (axl, axr))

post_processing/pylbo/utilities/toolbox.py

@@ -95,3 +95,56 @@ def transform_to_numpy(obj: any) -> np.ndarray:
     elif isinstance(obj, np.ndarray):
         return obj
     return np.asarray([obj])
+
+
+def solve_cubic_exact(a, b, c, d):
+    """
+    Solves a given cubic polynomial of the form
+    :math:`ax^3 + bx^2 + cx + d = 0` using the analytical cubic root formula
+    instead of the general `numpy.roots` routine.
+    From `StackOverflow <https://math.stackexchange.com/questions
+    15865why-not-write-the-solutions-of-a-cubic-this-way/18873#18873/>`_.
+
+    Parameters
+    ----------
+    a : int, float, complex
+        Cubic coefficient.
+    b : int, float, complex
+        Quadratic coefficient.
+    c : int, float, complex
+        Linear coefficient.
+    d : int, float, complex
+        Constant term
+
+    Returns
+    -------
+    roots : np.ndarray(ndim=3, dtype=complex)
+        The three roots of the cubic polynomial as a Numpy array.
+    """
+
+    if a == 0:
+        raise ValueError("cubic coefficient may not be zero")
+    p = b / a
+    q = c / a
+    r = d / a
+    Aterm = (
+        -2 * p ** 3
+        + 9 * p * q
+        - 27 * r
+        + 3
+        * np.sqrt(3)
+        * np.sqrt(
+            -(p ** 2) * q ** 2
+            + 4 * q ** 3
+            + 4 * p ** 3 * r
+            - 18 * p * q * r
+            + 27 * r ** 2
+        )
+    ) ** (1 / 3) / (3 * 2 ** (1 / 3))
+    Bterm = (-(p ** 2) + 3 * q) / (9 * Aterm)
+    cterm_min = (-1 - np.sqrt(3) * 1j) / 2
+    cterm_pos = (-1 + np.sqrt(3) * 1j) / 2
+    x1 = -p / 3 + Aterm - Bterm
+    x2 = -p / 3 + cterm_min * Aterm - cterm_pos * Bterm
+    x3 = -p / 3 + cterm_pos * Aterm - cterm_min * Bterm
+    return np.array([x1, x2, x3], dtype=complex)

post_processing/pylbo/visualisation/continua.py

@@ -44,7 +44,11 @@ def calculate_continua(ds):
 
     # Thermal continuum calculation
     # wave vector parallel to magnetic field, uses vector projection and scale factor
-    kpara = (k2 * B02 / ds.scale_factor + k3 * B03) / B0
+    if not all(np.isclose(B0, 0, atol=1e-10)):
+        # take care of hydro cases, this will throw a division by 0 warning
+        kpara = (k2 * B02 / ds.scale_factor + k3 * B03) / B0
+    else:
+        kpara = 0
     cs2 = gamma * p / rho  # sound speed
     vA2 = B0 ** 2 / rho  # Alfvén speed
     ci2 = p / rho  # isothermal sound speed
@@ -87,36 +91,45 @@ def calculate_continua(ds):
         thermal = np.empty(shape=len(ds.grid_gauss), dtype=complex)
         slow_min = np.empty_like(thermal)
         slow_plus = np.empty_like(thermal)
+        warning_logged = False
         for idx, _ in enumerate(ds.grid_gauss):
             solutions = np.roots([coeff3[idx], coeff2[idx], coeff1[idx], coeff0[idx]])
             # create mask for purely imaginary solutions
             mask = np.isclose(abs(solutions.real), 0)
-            imag_sol = solutions[mask]
-            # extract slow continuum solutions
-            s_neg, s_pos = np.sort_complex(solutions[np.invert(mask)])
-            slow_min[idx] = s_neg
-            slow_plus[idx] = s_pos
-            # process thermal continuum solution
-            if (imag_sol == 0).all():
-                # if all solutions are zero (no thermal continuum), then append zero
-                w = 0
+            # this should only have 1 True element, the thermal solution
+            if np.count_nonzero(mask) == 1:
+                s_neg, s_pos = np.sort_complex(solutions[np.invert(mask)])
             else:
-                # else there should be ONE value that is nonzero,
-                # this is the thermal continuum value.
-                # Filter out non-zero solution and try to unpack.
-                # If there is more than one non-zero
-                # value unpacking fails, this is a sanity check so we raise an error.
-                try:
-                    (w,) = imag_sol[imag_sol != 0]  # this is an array, so unpack with ,
-                except ValueError:
+                if not warning_logged:
                     pylboLogger.warning(
-                        f"Something went wrong, more than one solution for the thermal "
-                        f"continuum was found (and there should be only one). "
-                        f"Setting value to zero. "
-                        f"Solutions found: {imag_sol[imag_sol != 0]}."
+                        f"encountered at least 1 point where the slow continuum has a "
+                        f"real value close to zero. \nContinuum solutions: {solutions}"
                     )
-                    w = 0
-            thermal[idx] = w
+                    warning_logged = True
+                mask = np.array([False] * 3, dtype=bool)
+                # check if this already happens on the first index
+                if idx == 0:
+                    # if so, assume largest mode is the thermal one and complex sort
+                    # slow+ and slow- solutions
+                    mask[solutions.imag.argmax()] = True
+                    pylboLogger.warning(
+                        f"this happened already at the first gridpoint, so we assumed "
+                        f"{solutions[mask]} to be the thermal mode"
+                    )
+                    s_neg, s_pos = np.sort_complex(solutions[np.invert(mask)])
+                else:
+                    # if not, get solution closest to previous thermal mode
+                    mask[np.abs(solutions - thermal[idx - 1]).argmin()] = True
+                    # retrieve closest slow- solution
+                    slow_solutions = solutions[np.invert(mask)]
+                    idx_neg = np.abs(slow_solutions - slow_min[idx - 1]).argmin()
+                    s_neg = slow_solutions[idx_neg]
+                    # remaining value is slow+ solution, remove slow- from array
+                    (s_pos,) = np.delete(slow_solutions, idx_neg)
+            # extract continuum solutions
+            (thermal[idx],) = solutions[mask]
+            slow_min[idx] = s_neg
+            slow_plus[idx] = s_pos
 
     # get doppler-shifted continua and return
     continua = {

post_processing/pylbo/visualisation/spectra.py

@@ -173,7 +173,7 @@ def add_continua(self, interactive=True):
             min_value = np.min(continuum)
             max_value = np.max(continuum)
             # check if continua are complex
-            if np.all(np.iscomplex(continuum)):
+            if np.any(np.iscomplex(continuum)):
                 item = self.ax.scatter(
                     continuum.real * self.x_scaling,
                     continuum.imag * self.y_scaling,