Cleanup

.github/issue_template.md

@@ -0,0 +1,27 @@
+## Issue description
+<!-- A short description of the issue. Also add label(s) indicating if it is a feature request, bug report or something else.
+Feel free to delete sections that are not applicable. -->
+
+## Feature request
+<!-- What would you like to see added to the code?
+If you want start a discussion instead, feel free to open a discussion in the corresponding repository tab or ask us through Slack. -->
+
+
+## Bug report
+**Minimal example for reproduction**
+<!-- Paste your code between the quotation marks below. -->
+```
+
+```
+
+**Actual result**
+<!-- Result obtained when the above piece of code is run. This can be a description, screenshot, console statements etc. -->
+
+**Expected result**
+<!-- Which result do you expect instead? -->
+
+**Version info**
+- Legolas version:
+- Pylbo version:
+- Python version:
+- Version of relevant libraries (matplotlib, numpy, etc.):

.github/pull_request_template.md

@@ -1,42 +1,19 @@
-Hi, thanks for the pull request! To make sure everything stays nice and consistent, here is a default template which you can use.
-It also includes a checklist of small things that are usually forgotten or overlooked. You may remove parts that are not applicable to the current PR
-(including these sentences).
-
-Cheers,
-
-Legolas dev team
-
 ## PR description
-
-> fill or delete
+<!-- A short description of the pull request. Feel free to delete sections that are not applicable.-->
 
 ## New features
+**Legolas**
+-
 
-> fill or delete
+**Pylbo**
+-
 
 ## Bugfixes
+**Legolas**
+-
 
-> fill or delete
-
+**Pylbo**
+-
 ## Checklist
-
-**Source code stuff**
-- [ ] everything is nicely formatted using [Black-like code style](https://black.readthedocs.io/en/stable/the_black_code_style.html)
-- [ ] use-statements should use `use mod_xxx, only:` as much as possible
-- [ ] things that can go in a separate module, should go in a separate module
-- [ ] code additions|changes are also added|changed in the docstrings
-
-**Testing stuff**
-- [ ] all tests pass
-- [ ] in case of new features, new tests are added
-
-**Website stuff**
-- [ ] relevant pages on the website are edited, if relevant
-- [ ] on the edited pages, the `last_modified_at` frontmatter key is updated
-- [ ] all links render as they should, check through a local `bundle exec jekyll serve` in the `docs` directory
-- [ ] additions to the parfile are also added to `parameters_file.md`
-
-**Release stuff**
-(only if this is a merge-in-master release)
-- [ ] bump version number
-- [ ] changelog pages on the website are updated
+- [ ] all tests pass (and for new features, new tests are added)
+- [ ] documentation has been updated (if applicable)

.gitignore

@@ -4,6 +4,7 @@ output/
 build/
 legolas
 test_legolas
+legolas_config.par
 
 *.DS_Store
 *.aux

legolas_config.par

@@ -1,58 +1,58 @@
 &gridlist
-  ! geometry = 'Cartesian'
-  ! geometry = 'cylindrical'
-  ! x_start  = 0.0d0
-  ! x_end    = 1.0d0
-  gridpoints = 51
-  mesh_accumulation = .false.
+  geometry = "cylindrical"
+  x_start  = 0.0d0
+  x_end    = 1.0d0
+  gridpoints = 50
 /
 
 &solvelist
   solver = "QR-invert"
   ! solver = "arnoldi"
   ! arpack_mode = "shift-invert"
-  ! number_of_eigenvalues = 10
+  ! number_of_eigenvalues = 50
   ! which_eigenvalues = "LI"
-  ! sigma = (0.0d0, 0.01d0)
+  ! sigma = (-0.13d0, 0.0d0)
 /
 
 &physicslist
-  ! hall_mhd = .true.
+  ! flow = .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_flow"       ! Cartesian
-  ! equilibrium_type = "suydam_cluster"               ! cylindrical
-  ! equilibrium_type = "rotating_plasma_cylinder"     ! cylindrical
-  ! equilibrium_type = "kelvin_helmholtz_cd"          ! cylindrical
-  ! equilibrium_type = "internal_kink"                ! cylindrical
-  ! equilibrium_type = "RTI_theta_pinch"              ! cylindrical
-  ! equilibrium_type = "discrete_alfven"              ! cylindrical
-  ! equilibrium_type = "MRI_accretion"                ! cylindrical
-  ! equilibrium_type = "interchange_modes"            ! Cartesian
-  ! equilibrium_type = "constant_current_tokamak"     ! cylindrical
-  ! equilibrium_type = "magnetothermal_instabilities"
-  ! equilibrium_type = "photospheric_flux_tube"
-  ! equilibrium_type = "coronal_flux_tube"
-  ! equilibrium_type = "resonant_absorption"
+  ! equilibrium_type = "adiabatic_homo"
+  ! equilibrium_type = "constant_current_tokamak"
+  ! equilibrium_type = "couette_flow"
+  ! equilibrium_type = "discrete_alfven"
   ! equilibrium_type = "gold_hoyle"
-  ! equilibrium_type = "rayleigh_taylor"
+  ! equilibrium_type = "gravito_acoustic"
+  ! equilibrium_type = "gravito_mhd"
+  ! equilibrium_type = "coronal_flux_tube"
+  ! equilibrium_type = "photospheric_flux_tube"
+  ! equilibrium_type = "interchange_modes"
+  ! equilibrium_type = "internal_kink"
   ! equilibrium_type = "kelvin_helmholtz"
+  ! equilibrium_type = "kelvin_helmholtz_cd"
+  ! equilibrium_type = "MRI_accretion"
+  ! equilibrium_type = "magnetothermal_instabilities"
+  ! equilibrium_type = "rayleigh_taylor"
   ! equilibrium_type = "RTI_KHI"
+  ! equilibrium_type = "RTI_theta_pinch"
+  ! equilibrium_type = "resistive_homo"
+  ! equilibrium_type = "resistive_tearing"
+  ! equilibrium_type = "resistive_tearing_flow"
+  ! equilibrium_type = "resonant_absorption"
+  ! equilibrium_type = "rotating_plasma_cylinder"
   ! equilibrium_type = "isothermal_atmosphere"
-  ! equilibrium_type = "harris_sheet"
-  ! equilibrium_type = "user_defined"
+  equilibrium_type = "suydam_cluster"
 
-  ! boundary_type = "wall"
+  ! equilibrium_type = "user_defined"
 /
 
 &savelist
-  write_matrices = .false.
-  write_eigenfunctions = .true.
   show_results = .true.
+  write_eigenfunctions = .true.
+  write_derived_eigenfunctions = .true.
+  ! write_eigenfunction_subset = .true.
+  ! eigenfunction_subset_center = (-0.13d0, 0.0d0)
+  ! eigenfunction_subset_radius = 0.02d0
 /

post_processing/pylbo/visualisation/continua.py

@@ -162,6 +162,8 @@ def __init__(self, interactive):
         super().__init__(interactive)
         self.continua_names = CONTINUA_NAMES
         self._continua_colors = CONTINUA_COLORS
+        self.marker = "."
+        self.markersize = 6
 
     @property
     def continua_colors(self):

post_processing/pylbo/visualisation/legend_interface.py

@@ -1,5 +1,6 @@
 import matplotlib.collections as mpl_collections
 import matplotlib.lines as mpl_lines
+import matplotlib.patches as mpl_patches
 import numpy as np
 from pylbo.utilities.toolbox import add_pickradius_to_item
 
@@ -107,6 +108,10 @@ def make_legend_pickable(self):
             add_pickradius_to_item(item=legend_item, pickradius=self.pickradius)
             # add an attribute to this artist to tell it's from a legend
             setattr(legend_item, "is_legend_item", True)
+            # make sure colourmapping is done properly, only for continua regions
+            if isinstance(legend_item, mpl_patches.Rectangle):
+                legend_item.set_facecolor(drawn_item.get_facecolor())
+                legend_item.set_edgecolor(drawn_item.get_edgecolor())
             # we make the regions invisible until clicked, or set visible as default
             if self._make_visible_by_default:
                 legend_item.set_alpha(self.alpha_point)

post_processing/pylbo/visualisation/spectra.py

@@ -187,45 +187,20 @@ def add_continua(self, interactive=True):
             continuum = self.dataset.continua[key]
             if self._c_handler.check_if_all_zero(continuum=continuum):
                 continue
-            min_value = np.min(continuum)
-            max_value = np.max(continuum)
-            # check if continua are complex
-            if np.any(np.iscomplex(continuum)):
-                item = self.ax.scatter(
-                    continuum.real * self.x_scaling,
-                    continuum.imag * self.y_scaling,
-                    marker=".",
-                    color=color,
-                    linewidth=self.markersize,
-                    alpha=self.alpha / 1.5,
-                    label=key,
-                )
-            elif self._c_handler.check_if_collapsed(continuum=continuum):
-                item = self.ax.scatter(
-                    min_value * self.x_scaling,
-                    0,
-                    marker=self._c_handler.marker,
-                    s=self._c_handler.markersize,
-                    c=color,
-                    alpha=self._c_handler.alpha_point,
-                    label=key,
-                )
-            else:
-                props = dict(
-                    facecolor=colors.to_rgba(color, self._c_handler.alpha_region),
-                    edgecolor=colors.to_rgba(color, self._c_handler.alpha_point),
-                    label=key,
-                )
-                if key == "thermal":
-                    item = self.ax.axhspan(
-                        min_value.imag * self.y_scaling,
-                        max_value * self.y_scaling,
-                        **props,
-                    )
-                else:
-                    item = self.ax.axvspan(
-                        min_value * self.x_scaling, max_value * self.x_scaling, **props
-                    )
+            realvals = continuum.real * self.x_scaling
+            imagvals = continuum.imag * self.y_scaling
+            if self._c_handler.check_if_collapsed(continuum=continuum):
+                realvals = np.min(realvals)
+                imagvals = 0
+            item = self.ax.scatter(
+                realvals,
+                imagvals,
+                marker=self._c_handler.marker,
+                linewidth=self._c_handler.markersize,
+                c=color,
+                alpha=self._c_handler.alpha_point,
+                label=key,
+            )
             self._c_handler.add(item)
         self._c_handler.legend = self.ax.legend(**self._c_handler.legend_properties)
         super().add_continua(self._c_handler.interactive)
@@ -251,11 +226,6 @@ def add_derived_eigenfunctions(self):
             )
         super()._enable_interface(handle=self._def_handler)
 
-    def _ensure_smooth_continuum(self, continuum):
-        # TODO: this method should split the continuum into multiple parts, such that
-        #       regions that lie far apart are not connected by a line.
-        pass
-
 
 class MultiSpectrumPlot(SpectrumFigure):
     """
@@ -432,9 +402,8 @@ def add_continua(self, interactive=True):
             self._c_handler.continua_names, self._c_handler.continua_colors
         ):
             # we skip duplicates if eigenvalues are squared
-            if self.use_squared_omega:
-                if key in ("slow-", "alfven-"):
-                    continue
+            if self.use_squared_omega and key in ("slow-", "alfven-"):
+                continue
             # retrieve continuum, calculate region boundaries
             continuum = self.dataseries.continua[key] ** self._w_pow
             if self.use_real_parts:
@@ -488,7 +457,7 @@ def add_derived_eigenfunctions(self):
             self._def_handler = DerivedEigenfunctionHandler(
                 self.dataseries, self._def_ax, self.ax
             )
-        super()._enable_interface(handle=self._pp_handler)
+        super()._enable_interface(handle=self._def_handler)
 
 
 class MergedSpectrumPlot(SpectrumFigure):

src/CMakeLists.txt

@@ -73,8 +73,6 @@ set (sources
         matrices/smod_conduction_matrix.f08
         matrices/smod_viscosity_matrix.f08
         matrices/smod_hall_matrix.f08
-        # physics/mod_viscosity.f08
-        # mod_boundary_conditions.f08
         solvers/mod_solvers.f08
         solvers/mod_arpack_type.f08
         solvers/smod_qr_invert.f08

src/equilibria/smod_equil_internal_kink_instability.f08

@@ -24,7 +24,7 @@
   module subroutine internal_kink_eq()
     use mod_equilibrium_params, only: cte_rho0, cte_v03, cte_p0, alpha
 
-    real(dp)      :: r, x
+    real(dp)      :: r, x, a0
     real(dp)      :: J0, J1, J2, DJ0, DJ1
     integer       :: i
 
@@ -33,36 +33,39 @@ module subroutine internal_kink_eq()
     )
     call initialise_grid()
 
+    a0 = x_end
     if (use_defaults) then ! LCOV_EXCL_START
       flow = .true.
       cte_rho0 = 1.0d0
       cte_v03  = 1.0d0
       cte_p0 = 9.0d0
-      alpha = 5.0d0 / x_end
+      alpha = 5.0d0 / a0
 
       k2 = 1.0d0
       k3 = 0.16d0 * alpha
     end if ! LCOV_EXCL_STOP
 
     do i = 1, gauss_gridpts
       r = grid_gauss(i)
-      x = r / x_end
+      x = r / a0
 
       J0   = bessel_jn(0, alpha * x)
       J1   = bessel_jn(1, alpha * x)
       J2   = bessel_jn(2, alpha * x)
       DJ0  = -alpha * J1
       DJ1  = alpha * (0.5d0 * J0 - 0.5d0 * J2)
 
-      rho_field % rho0(i) = cte_rho0 * (1.0d0 - x**2)
-      v_field % v03(i)    = cte_v03 * (1.0d0 - x**2)
+      rho_field % rho0(i) = cte_rho0 * (1.0d0 - x**2 / a0**2)
+      v_field % v03(i)    = cte_v03 * (1.0d0 - x**2 / a0**2)
       B_field % B02(i)    = J1
       B_field % B03(i)    = J0
       B_field % B0(i)     = sqrt((B_field % B02(i))**2 + (B_field % B03(i))**2)
       T_field % T0(i)     = cte_p0 / (rho_field % rho0(i))
 
-      rho_field % d_rho0_dr(i) = -2.0d0 * cte_rho0 * x
-      T_field % d_T0_dr(i)     = 2.0d0 * x * cte_p0 / (cte_rho0 * (1.0d0 - x**2)**2)
+      rho_field % d_rho0_dr(i) = -2.0d0 * cte_rho0 * x / a0
+      T_field % d_T0_dr(i)     = 2.0d0 * x * cte_p0 / ( &
+        a0**2 * cte_rho0 * (1.0d0 - x**2)**2 &
+      )
       v_field % d_v03_dr(i)    = -2.0d0 * cte_v03 * x
       B_field % d_B02_dr(i)    = DJ1
       B_field % d_B03_dr(i)    = DJ0

src/main.f08

@@ -62,7 +62,7 @@ subroutine initialisation()
     use mod_input, only: read_parfile, get_parfile
     use mod_equilibrium, only: initialise_equilibrium, set_equilibrium
     use mod_logging, only: print_logo
-    use mod_global_variables, only: viscosity, hall_mhd
+    use mod_global_variables, only: hall_mhd
 
     character(len=str_len)  :: parfile
     integer   :: nb_evs
@@ -89,14 +89,10 @@ subroutine initialisation()
     call set_equilibrium()
 
     ! TODO: remove this warning when fully tested
-    if (viscosity) then
-      call log_message( &
-        "using viscous MHD, note that this is not yet fully tested!", level="warning" &
-      )
-    end if
     if (hall_mhd) then
       call log_message( &
-        "using Hall MHD, note that this does not yet work properly!", level="warning" &
+        "using Hall MHD, note that some features are not yet **fully** tested!", &
+        level="warning" &
       )
     end if