diff --git a/.github/issue_template.md b/.github/issue_template.md
new file mode 100644
index 00000000..13916e48
--- /dev/null
+++ b/.github/issue_template.md
@@ -0,0 +1,27 @@
+## Issue description
+
+
+## Feature request
+
+
+
+## Bug report
+**Minimal example for reproduction**
+
+```
+
+```
+
+**Actual result**
+
+
+**Expected result**
+
+
+**Version info**
+- Legolas version:
+- Pylbo version:
+- Python version:
+- Version of relevant libraries (matplotlib, numpy, etc.):
\ No newline at end of file
diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md
index 8205178b..91922b56 100644
--- a/.github/pull_request_template.md
+++ b/.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
+
## 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)
diff --git a/.gitignore b/.gitignore
index b6662d91..c4be27f1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,7 @@ output/
build/
legolas
test_legolas
+legolas_config.par
*.DS_Store
*.aux
diff --git a/legolas_config.par b/legolas_config.par
index bd099444..8a74e3e2 100644
--- a/legolas_config.par
+++ b/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
/
diff --git a/post_processing/pylbo/visualisation/continua.py b/post_processing/pylbo/visualisation/continua.py
index 02c01121..4c1048e3 100644
--- a/post_processing/pylbo/visualisation/continua.py
+++ b/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):
diff --git a/post_processing/pylbo/visualisation/legend_interface.py b/post_processing/pylbo/visualisation/legend_interface.py
index f22053c2..c831aafd 100644
--- a/post_processing/pylbo/visualisation/legend_interface.py
+++ b/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)
diff --git a/post_processing/pylbo/visualisation/spectra.py b/post_processing/pylbo/visualisation/spectra.py
index 180ba47d..e65ab48c 100644
--- a/post_processing/pylbo/visualisation/spectra.py
+++ b/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):
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ed59885e..248fe69a 100644
--- a/src/CMakeLists.txt
+++ b/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
diff --git a/src/equilibria/smod_equil_internal_kink_instability.f08 b/src/equilibria/smod_equil_internal_kink_instability.f08
index f81c0018..688ea44e 100644
--- a/src/equilibria/smod_equil_internal_kink_instability.f08
+++ b/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,12 +33,13 @@ 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
@@ -46,7 +47,7 @@ module subroutine internal_kink_eq()
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)
@@ -54,15 +55,17 @@ module subroutine internal_kink_eq()
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
diff --git a/src/main.f08 b/src/main.f08
index 96c33a46..e11287ae 100644
--- a/src/main.f08
+++ b/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
diff --git a/src/mod_boundary_conditions.f08 b/src/mod_boundary_conditions.f08
deleted file mode 100644
index c2d49498..00000000
--- a/src/mod_boundary_conditions.f08
+++ /dev/null
@@ -1,282 +0,0 @@
-! =============================================================================
-!> Imposes boundary conditions on the finite element matrices after assembly.
-!! Both natural and essential conditions are handled.
-module mod_boundary_conditions
- use mod_global_variables, only: dp, matrix_gridpts, dim_quadblock, dim_subblock
- implicit none
-
- private
-
- !> boolean to check presence of perpendicular thermal conduction
- logical, save :: kappa_perp_is_zero
-
- public :: apply_boundary_conditions
-
-contains
-
-
- !> Main routine to handle the boundary conditions. This subroutine first performs some checks
- !! and then handles the natural and essential boundary conditions.
- !! @note Perpendicular thermal conduction is hard-checked, that is, not just the
- !! global logical. There must be a value in the corresponding array that is non-zero.
- subroutine apply_boundary_conditions(matrix_A, matrix_B)
- use mod_equilibrium, only: kappa_field
- ! use mod_viscosity, only: viscosity_boundaries
- use mod_global_variables, only: dp_LIMIT
- use mod_equilibrium, only: kappa_field
-
- !> the A-matrix with boundary conditions imposed on exit
- complex(dp), intent(inout) :: matrix_A(matrix_gridpts, matrix_gridpts)
- !> the B-matrix with boundary conditions imposed on exit
- real(dp), intent(inout) :: matrix_B(matrix_gridpts, matrix_gridpts)
- complex(dp) :: quadblock(dim_quadblock, dim_quadblock)
- ! complex(dp) :: quadblock_visc(dim_quadblock, dim_quadblock)
- ! complex(dp) :: quadblock_Hall(dim_quadblock, dim_quadblock)
- integer :: idx_end_left, idx_start_right
- ! real(dp) :: hf
-
- ! check if perpendicular thermal conduction is present
- kappa_perp_is_zero = .true.
- if (any(abs(kappa_field % kappa_perp) > dp_LIMIT)) then
- kappa_perp_is_zero = .false.
- end if
-
- ! end of index first-gridpoint quadblock
- idx_end_left = dim_quadblock
- ! start of index last-gridpoint quadblock
- idx_start_right = matrix_gridpts - dim_quadblock + 1
-
- ! matrix B left-edge quadblock
- quadblock = matrix_B(1:idx_end_left, 1:idx_end_left)
- call essential_boundaries(quadblock, edge='l_edge', matrix='B')
- matrix_B(1:idx_end_left, 1:idx_end_left) = real(quadblock)
- ! matrix B right-edge quadblock
- quadblock = matrix_B(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts)
- call essential_boundaries(quadblock, edge='r_edge', matrix='B')
- matrix_B(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts) = real(quadblock)
-
- ! matrix A left-edge quadblock
- quadblock = matrix_A(1:idx_end_left, 1:idx_end_left)
- call essential_boundaries(quadblock, edge='l_edge', matrix='A')
- call natural_boundaries(quadblock, edge='l_edge')
- ! if (viscosity) then
- ! call viscosity_boundaries(quadblock_visc, edge='l_edge')
- ! quadblock = quadblock + quadblock_visc
- ! end if
- ! if (hall_mhd) then
- ! call hall_boundaries(quadblock_Hall, kappa_perp_is_zero, edge='l_edge')
- ! hf = hall_field % hallfactor(1)
- ! quadblock = quadblock - hf * quadblock_Hall
- ! end if
- matrix_A(1:idx_end_left, 1:idx_end_left) = quadblock
- ! matrix A right-edge quadblock
- quadblock = matrix_A(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts)
- call essential_boundaries(quadblock, edge='r_edge', matrix='A')
- call natural_boundaries(quadblock, edge='r_edge')
- ! if (hall_mhd) then
- ! call hall_boundaries(quadblock_Hall, kappa_perp_is_zero, edge='r_edge')
- ! hf = hall_field % hallfactor(gauss_gridpts)
- ! quadblock = quadblock - hf * quadblock_Hall
- ! end if
- matrix_A(idx_start_right:matrix_gridpts, idx_start_right:matrix_gridpts) = quadblock
- end subroutine apply_boundary_conditions
-
-
- !> Imposes essential boundary conditions on a given matrix, that is,
- !! the ones that have to be implemented explicitly.
- !! For a wall, that means handling \(v1, a2\) and \(a3\) (and \(T1\) if there is
- !! perpendicular thermal conduction).
- !! @note The contribution from the 0 quadratic basis function automatically
- !! zeroes out the odd rows/columns for the quadratic variables on the
- !! left edge. These are explicitly handled. @endnote
- !! @note In the case of the Arnoldi solver, setting boundary conditions to
- !! 1e20 throws off the solver such that it returns completely wrong results.
- !! Setting them to a "normal" value, 0 here, seems to fix it. We now have
- !! an exact match with the QR/QZ solvers in this case.
- !! Eigenvalues will be introduced at 0 instead of at 1, but this does not
- !! affect the final spectrum in any way. @endnote
- !! @warning Throws an error if boundary_type is not known,
- !! or if edge is not known.
- subroutine essential_boundaries(quadblock, edge, matrix)
- use mod_global_variables, only: boundary_type, solver, viscosity, geometry, &
- coaxial
- use mod_logging, only: log_message
-
- !> the quadblock corresponding to the left/right edge
- complex(dp), intent(inout) :: quadblock(dim_quadblock, dim_quadblock)
- !> the edge, either "l_edge" or "r_edge"
- character(len=6), intent(in) :: edge
- !> the matrix, either "A" or "B"
- character, intent(in) :: matrix
-
- complex(dp) :: diagonal_factor
- integer :: i, j, qua_zeroes(5), wall_idx_left(4), wall_idx_right(4)
- integer :: noslip_idx_left(2), noslip_idx_right(2)
-
- if (matrix == 'B') then
- diagonal_factor = (1.0d0, 0.0d0)
- else if (matrix == 'A') then
- if (solver == "arnoldi") then
- diagonal_factor = (0.0d0, 0.0d0)
- else
- diagonal_factor = (1.0d20, 0.0d0)
- end if
- else
- call log_message("essential boundaries: invalid matrix argument", level='error')
- end if
-
- ! explicitly handle zeroed out rows/cols due to 0 term quadratic basis function
- qua_zeroes = [1, 5, 7, 9, 11]
- if (edge == 'l_edge') then
- do i = 1, size(qua_zeroes)
- j = qua_zeroes(i)
- quadblock(j, j) = diagonal_factor
- end do
- end if
-
- ! wall/regularity conditions: handling of v1, a2 and a3 (and T if conduction).
- ! v1, a2 and a3 are cubic elements, so omit non-zero basis functions (odd rows/columns)
- ! T is a quadratic element, so omit even row/columns
- wall_idx_left = [3, 13, 15, 10]
- wall_idx_right = [19, 29, 31, 26]
- ! For a no-slip condition for viscous fluids v2 and v3 should equal the wall's
- ! tangential velocities, here zero. (Quadratic elements, so omit even row/columns)
- noslip_idx_left = [6, 8]
- noslip_idx_right = [22, 24]
-
- select case(boundary_type)
- case('wall')
- if (edge == 'l_edge') then
- ! left regularity/wall conditions
- do i = 1, size(wall_idx_left)
- j = wall_idx_left(i)
- if (j == 10 .and. kappa_perp_is_zero) then
- cycle
- end if
- quadblock(j, :) = (0.0d0, 0.0d0)
- quadblock(:, j) = (0.0d0, 0.0d0)
- quadblock(j, j) = diagonal_factor
- end do
- ! No-slip condition does not apply at left edge for a cylindrical
- ! geometry unless two coaxial walls are used
- if (viscosity .and. (geometry == 'Cartesian' .or. coaxial)) then
- do i = 1, size(noslip_idx_left)
- j = noslip_idx_left(i)
- quadblock(j, :) = (0.0d0, 0.0d0)
- quadblock(:, j) = (0.0d0, 0.0d0)
- quadblock(j, j) = diagonal_factor
- end do
- end if
- else if (edge == 'r_edge') then
- do i = 1, size(wall_idx_right)
- j = wall_idx_right(i)
- if ((j == 26) .and. kappa_perp_is_zero) then
- cycle
- end if
- quadblock(j, :) = (0.0d0, 0.0d0)
- quadblock(:, j) = (0.0d0, 0.0d0)
- quadblock(j, j) = diagonal_factor
- end do
- if (viscosity) then
- do i = 1, size(noslip_idx_right)
- j = noslip_idx_right(i)
- quadblock(j, :) = (0.0d0, 0.0d0)
- quadblock(:, j) = (0.0d0, 0.0d0)
- quadblock(j, j) = diagonal_factor
- end do
- end if
- else
- call log_message("essential boundaries: invalid edge argument", level='error')
- end if
-
- case default
- call log_message( "essential boundaries: invalid boundary_type", level='error')
- end select
- end subroutine essential_boundaries
-
-
- !> Imposes the natural boundary conditions, that is, the ones originating
- !! from the weak formulation. This is only applicable for the A-matrix.
- !! @note For now only the terms of a solid wall boundary are implemented, hence
- !! we only have additional resistive and conductive terms in the energy equation.
- !! If perpendicular thermal conduction is included we have \(T1 = 0\), such that
- !! then there are no surface terms to be added and we simply return.
- !! If free boundary conditions are considered (e.g. vacuum-wall) additional terms
- !! have to be added. @endnote
- !! @warning Throws an error if edge is unknown.
- subroutine natural_boundaries(quadblock, edge)
- use mod_global_variables, only: boundary_type, gauss_gridpts, gamma_1, ic, dim_subblock
- use mod_logging, only: log_message
- use mod_equilibrium, only: B_field, eta_field
- use mod_grid, only: eps_grid, d_eps_grid_dr
- use mod_equilibrium_params, only: k2, k3
-
- !> the quadblock corresponding to the left/right edge
- complex(dp), intent(inout) :: quadblock(dim_quadblock, dim_quadblock)
- !> the edge, either "l_edge" or "r_edge"
- character(len=6), intent(in) :: edge
-
- complex(dp), allocatable :: surface_terms(:)
- real(dp) :: eps, d_eps_dr, eta, B02, dB02, B03, dB03, drB02
- integer, allocatable :: positions(:, :)
- integer :: idx, i
-
- if (.not. boundary_type == 'wall') then
- call log_message('natural boundaries: only wall is implemented!', level='error')
- end if
-
- ! For a wall with perpendicular thermal conduction we also have the condition T1 = 0.
- ! Hence, in that case there are no surface terms to be added, and we simply return.
- ! So for now we only have resistive terms in the energy equation.
- if (.not. kappa_perp_is_zero) then
- return
- end if
-
- if (edge == 'l_edge') then
- idx = 1
- else if (edge == 'r_edge') then
- idx = gauss_gridpts
- else
- call log_message('natural boundaries: wrong edge supplied' // edge, level='error')
- end if
-
- ! retrieve variables at current edge
- eps = eps_grid(idx)
- d_eps_dr = d_eps_grid_dr(idx)
- B02 = B_field % B02(idx)
- dB02 = B_field % d_B02_dr(idx)
- B03 = B_field % B03(idx)
- dB03 = B_field % d_B03_dr(idx)
- drB02 = d_eps_dr * B02 + eps * dB02
- eta = eta_field % eta(idx)
-
- allocate(positions(3, 2))
- allocate(surface_terms(3))
-
- ! surface term for element (5, 6)
- surface_terms(1) = 2.0d0 * ic * gamma_1 * (1.0d0 / eps) * eta * (k3 * drB02 - k2 * dB03)
- positions(1, :) = [5, 6]
- ! surface term for element (5, 7)
- surface_terms(2) = 2.0d0 * ic * gamma_1 * (1.0d0 / eps) * eta * dB03
- positions(2, :) = [5, 7]
- ! surface term for element (5, 8)
- surface_terms(3) = -2.0d0 * ic * gamma_1 * (1.0d0 / eps) * eta * drB02
- positions(3, :) = [5, 8]
-
- ! l_edge: add to bottom-right of 2x2 block, for top-left subblock only
- ! r_edge: add to bottom-right of 2x2 block, for bottom-right subblock only
- if (edge == 'l_edge') then
- positions = 2 * positions
- else if (edge == 'r_edge') then
- positions = 2 * positions + dim_subblock
- end if
-
- do i = 1, size(surface_terms)
- quadblock(positions(i, 1), positions(i, 2)) = quadblock(positions(i, 1), positions(i, 2)) + surface_terms(i)
- end do
-
- deallocate(positions)
- deallocate(surface_terms)
- end subroutine natural_boundaries
-
-end module mod_boundary_conditions
diff --git a/src/mod_hallmhd.f08 b/src/mod_hallmhd.f08
deleted file mode 100644
index ea1ee552..00000000
--- a/src/mod_hallmhd.f08
+++ /dev/null
@@ -1,632 +0,0 @@
-! =============================================================================
-!> This module handles everything related to the Hall term, computing
-!! the normalised Hall factor, the contributions to the A matrix, and natural
-!! boundary conditions, when necessary.
-module mod_hallmhd
- use mod_global_variables, only: dp, dim_quadblock
-
-implicit none
-
-private
-
-public :: set_hallfactor
-public :: get_hallterm
-public :: get_hallterm_Bmat
-public :: hall_boundaries
-
-contains
-
- !> Retrieves the normalised Hall factor as described by Porth et al. (2014),
- !! with a dropoff at the boundary
- subroutine set_hallfactor(hall_field)
- use mod_grid, only: grid_gauss
- use mod_physical_constants, only: dpi
- use mod_global_variables, only: cgs_units, gauss_gridpts, dropoff_edge_dist, &
- dropoff_width, hall_dropoff, inertia_dropoff
- use mod_physical_constants, only: mp_cgs, mp_si, ec_cgs, ec_si, me_cgs, me_si
- use mod_units, only: unit_velocity, unit_length, unit_magneticfield
- use mod_types, only: hall_type
-
- type (hall_type), intent(inout) :: hall_field
-
- real(dp) :: sleft, sright, width, hallval, inertiaval
- real(dp) :: x, shift, stretch, shift2, stretch2
- integer :: i
-
- width = dropoff_width
- if (cgs_units) then
- hallval = (mp_cgs * unit_velocity) / (ec_cgs * unit_length * unit_magneticfield)
- inertiaval = (mp_cgs * me_cgs * unit_velocity**2) / (ec_cgs * unit_length * unit_magneticfield)**2
- else
- hallval = (mp_si * unit_velocity) / (ec_si * unit_length * unit_magneticfield)
- inertiaval = (mp_si * me_si * unit_velocity**2) / (ec_si * unit_length * unit_magneticfield)**2
- end if
-
- sleft = grid_gauss(1) + 0.5d0 * width + dropoff_edge_dist
- sright = grid_gauss(gauss_gridpts) - dropoff_edge_dist - 0.5d0 * width
-
- if (hall_dropoff) then
- shift = hallval * tanh(-dpi) / (tanh(-dpi) - tanh(dpi))
- stretch = hallval / (tanh(dpi) - tanh(-dpi))
-
- do i = 1, gauss_gridpts
- x = grid_gauss(i)
- if (sleft - 0.5d0*width <= x .and. x <= sleft + 0.5d0*width) then
- hall_field % hallfactor(i) = shift + stretch * tanh(2.0d0 * dpi * (x - sleft) / width)
- else if (sleft + 0.5d0*width < x .and. x < sright - 0.5d0*width) then
- hall_field % hallfactor(i) = hallval
- else if (sright - 0.5d0*width <= x .and. x <= sright + 0.5d0*width) then
- hall_field % hallfactor(i) = shift + stretch * tanh(2.0d0 * dpi * (sright - x) / width)
- else
- hall_field % hallfactor(i) = 0.0d0
- end if
- end do
- else
- hall_field % hallfactor = hallval
- end if
-
- if (inertia_dropoff) then
- shift2 = inertiaval * tanh(-dpi) / (tanh(-dpi) - tanh(dpi))
- stretch2 = inertiaval / (tanh(dpi) - tanh(-dpi))
-
- do i = 1, gauss_gridpts
- x = grid_gauss(i)
- if (sleft - 0.5d0*width <= x .and. x <= sleft + 0.5d0*width) then
- hall_field % inertiafactor(i) = shift2 + stretch2 * tanh(2.0d0 * dpi * (x - sleft) / width)
- else if (sleft + 0.5d0*width < x .and. x < sright - 0.5d0*width) then
- hall_field % inertiafactor(i) = inertiaval
- else if (sright - 0.5d0*width <= x .and. x <= sright + 0.5d0*width) then
- hall_field % inertiafactor(i) = shift2 + stretch2 * tanh(2.0d0 * dpi * (sright - x) / width)
- else
- hall_field % inertiafactor(i) = 0.0d0
- end if
- end do
- else
- hall_field % inertiafactor = inertiaval
- end if
- end subroutine set_hallfactor
-
- !> Retrieves the A-matrix Hall contribution for a given quadblock corresponding
- !! to a particular point in the grid and a particular Gaussian weight.
- !! If hall_mhd = True, this routine is called n_gauss times
- !! for every grid interval.
- subroutine get_hallterm(gauss_idx, eps, d_eps_dr, curr_weight, quadblock_Hall, &
- h_quadratic, dh_quadratic_dr, h_cubic, dh_cubic_dr)
- use mod_global_variables, only: dp_LIMIT, elec_pressure
- use mod_equilibrium_params, only: k2, k3
- use mod_equilibrium, only: rho_field, B_field, T_field
- use mod_make_subblock, only: subblock, reset_factors, reset_positions
-
- !> current index in the Gaussian grid
- integer, intent(in) :: gauss_idx
- !> current value for the scale factor epsilon
- real(dp), intent(in) :: eps
- !> current value for the derivative of epsilon
- real(dp), intent(in) :: d_eps_dr
- !> current weight in the Gaussian quadrature
- real(dp), intent(in) :: curr_weight
- !> the A-quadblock for a particular grid interval
- complex(dp), intent(out) :: quadblock_Hall(dim_quadblock, dim_quadblock)
- !> array containing the 4 quadratic basis functions
- real(dp), intent(in) :: h_quadratic(4)
- !> array containing the derivatives of the 4 quadratic basis functions
- real(dp), intent(in) :: dh_quadratic_dr(4)
- !> array containing the 4 cubic basis functions
- real(dp), intent(in) :: h_cubic(4)
- !> array containing the derivatives of the 4 cubic basis functions
- real(dp), intent(in) :: dh_cubic_dr(4)
-
- complex(dp), allocatable :: factors(:)
- integer, allocatable :: positions(:, :)
-
- real(dp) :: eps_inv, rho0_inv
- real(dp) :: rho0, B01, B02, B03, T0
- real(dp) :: drho0, drB02, dB02_r, dB03, dB02, dT0
-
- quadblock_Hall = (0.0d0, 0.0d0)
- eps_inv = 1.0d0 / eps
-
- ! Equilibrium quantities
- rho0 = rho_field % rho0(gauss_idx)
- drho0 = rho_field % d_rho0_dr(gauss_idx)
- rho0_inv = 1.0d0 / rho0
-
- T0 = T_field % T0(gauss_idx)
- dT0 = T_field % d_T0_dr(gauss_idx)
-
- ! To be uncommented when B01 is implemented
- B01 = 0.0d0 !B_field % B01(gauss_idx)
- B02 = B_field % B02(gauss_idx)
- dB02 = B_field % d_B02_dr(gauss_idx)
- B03 = B_field % B03(gauss_idx)
- dB03 = B_field % d_B03_dr(gauss_idx)
-
- ! Calculate derivatives eps*B02, B02/eps
- drB02 = d_eps_dr * B02 + eps * dB02
- dB02_r = eps_inv * dB02 - d_eps_dr * eps_inv**2 * B02
-
- ! Setup of matrix elements
-
- ! Quadratic * Quadratic
- call reset_factors(factors, 3)
- call reset_positions(positions, 3)
- ! A(6, 1)
- factors(1) = - eps_inv * rho0_inv**2 * (B02 * drB02 * eps_inv + B03 * dB03)
- positions(1, :) = [6, 1]
- ! A(6, 5)
- factors(2) = 0.0d0
- positions(2, :) = [6, 5]
- ! A(6, 6)
- factors(3) = - drho0 * rho0_inv**2 * (k2 * B03 * eps_inv - k3 * B02) &
- + k3 * rho0_inv * (drB02 * eps_inv - eps * dB02_r)
- positions(3, :) = [6, 6]
- if (elec_pressure) then
- factors(1) = factors(1) - dT0 * eps_inv * rho0_inv
- factors(2) = factors(2) + drho0 * eps_inv * rho0_inv
- end if
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_quadratic, h_quadratic)
-
- ! d(Quadratic)/dr * Quadratic
- call reset_factors(factors, 3)
- call reset_positions(positions, 3)
- ! A(6, 1)
- factors(1) = 0.0d0
- positions(1, :) = [6, 1]
- ! A(6, 5)
- factors(2) = 0.0d0
- positions(2, :) = [6, 5]
- ! A(6, 6)
- factors(3) = rho0_inv * (k2 * B03 * eps_inv - k3 * B02)
- positions(3, :) = [6, 6]
- if (elec_pressure) then
- factors(1) = factors(1) - T0 * eps_inv * rho0_inv
- factors(2) = factors(2) - eps_inv
- end if
- call subblock(quadblock_Hall, factors, positions, curr_weight, dh_quadratic_dr, h_quadratic)
-
- ! Quadratic * Cubic
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(6, 7)
- factors(1) = - k3 * eps_inv * rho0_inv * (k2 * B02 * eps_inv + k3 * B03)
- positions(1, :) = [6, 7]
- ! A(6, 8)
- factors(2) = k2 * eps_inv * rho0_inv * (k2 * B02 * eps_inv + k3 * B03)
- positions(2, :) = [6, 8]
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_quadratic, h_cubic)
-
- ! Quadratic * d(Cubic)/dr
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(6, 7)
- factors(1) = drho0 * rho0_inv**2 * B03 * eps_inv
- positions(1, :) = [6, 7]
- ! A(6, 8)
- factors(2) = rho0_inv * (eps * dB02_r - drB02 * eps_inv - drho0 * rho0_inv * B02)
- positions(2, :) = [6, 8]
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_quadratic, dh_cubic_dr)
-
- ! d(Quadratic)/dr * d(Cubic)/dr
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(6, 7)
- factors(1) = - B03 * eps_inv * rho0_inv
- positions(1, :) = [6, 7]
- ! A(6, 8)
- factors(2) = B02 * rho0_inv
- positions(2, :) = [6, 8]
- call subblock(quadblock_Hall, factors, positions, curr_weight, dh_quadratic_dr, dh_cubic_dr)
-
- ! Cubic * Quadratic
- call reset_factors(factors, 6)
- call reset_positions(positions, 6)
- ! A(7, 1)
- factors(1) = 0.0d0
- positions(1, :) = [7, 1]
- ! A(7, 5)
- factors(2) = 0.0d0
- positions(2, :) = [7, 5]
- ! A(7, 6)
- factors(3) = - B03 * rho0_inv * ((k2 * eps_inv)**2 + k3**2)
- positions(3, :) = [7, 6]
- ! A(8, 1)
- factors(4) = 0.0d0
- positions(4, :) = [8, 1]
- ! A(8, 5)
- factors(5) = 0.0d0
- positions(5, :) = [8, 5]
- ! A(8, 6)
- factors(6) = B02 * rho0_inv * ((k2 * eps_inv)**2 + k3**2)
- positions(6, :) = [8, 6]
- if (elec_pressure) then
- factors(1) = factors(1) + k2 * T0 * rho0_inv * eps_inv**2
- factors(2) = factors(2) + k2 * eps_inv**2
- factors(4) = factors(4) + k3 * T0 * rho0_inv * eps_inv
- factors(5) = factors(5) + k3 * eps_inv
- end if
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_quadratic)
-
- ! Cubic * Cubic
- call reset_factors(factors, 4)
- call reset_positions(positions, 4)
- ! A(7, 7)
- factors(1) = k3 * drB02 * eps_inv**2 * rho0_inv
- positions(1, :) = [7, 7]
- ! A(7, 8)
- factors(2) = - k2 * drB02 * eps_inv**2 * rho0_inv
- positions(2, :) = [7, 8]
- ! A(8, 7)
- factors(3) = k3 * dB03 * eps_inv * rho0_inv
- positions(3, :) = [8, 7]
- ! A(8, 8)
- factors(4) = - k2 * dB03 * eps_inv * rho0_inv
- positions(4, :) = [8, 8]
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_cubic)
-
- ! Cubic * d(Cubic)/dr
- call reset_factors(factors, 4)
- call reset_positions(positions, 4)
- ! A(7, 7)
- factors(1) = k2 * B03 * eps_inv**2 * rho0_inv
- positions(1, :) = [7, 7]
- ! A(7, 8)
- factors(2) = k3 * B03 * rho0_inv
- positions(2, :) = [7, 8]
- ! A(8, 7)
- factors(3) = - k2 * B02 * eps_inv**2 * rho0_inv
- positions(3, :) = [8, 7]
- ! A(8, 8)
- factors(4) = - k3 * B02 * rho0_inv
- positions(4, :) = [8, 8]
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, dh_cubic_dr)
-
- deallocate(factors)
- deallocate(positions)
-
- if (B01 > dp_LIMIT) then
- call get_hallterm_B01(gauss_idx, eps, d_eps_dr, curr_weight, &
- quadblock_Hall, h_quadratic, h_cubic, dh_cubic_dr)
- end if
- end subroutine get_hallterm
-
- !> Retrieves the A-matrix Hall contribution for a given quadblock corresponding
- !! to a particular point in the grid and a particular Gaussian weight.
- !! If hall_mhd = True, this routine is called n_gauss times
- !! for every grid interval.
- subroutine get_hallterm_Bmat(gauss_idx, eps, d_eps_dr, curr_weight, quadblock_HallB, &
- h_quadratic, h_cubic, dh_cubic_dr)
- use mod_equilibrium_params, only: k2, k3
- use mod_equilibrium, only: rho_field
- use mod_make_subblock, only: subblock, reset_factors, reset_positions
-
- !> current index in the Gaussian grid
- integer, intent(in) :: gauss_idx
- !> current value for the scale factor epsilon
- real(dp), intent(in) :: eps
- !> current value for the derivative of epsilon
- real(dp), intent(in) :: d_eps_dr
- !> current weight in the Gaussian quadrature
- real(dp), intent(in) :: curr_weight
- !> the A-quadblock for a particular grid interval
- complex(dp), intent(out) :: quadblock_HallB(dim_quadblock, dim_quadblock)
- !> array containing the 4 quadratic basis functions
- real(dp), intent(in) :: h_quadratic(4)
- !> array containing the 4 cubic basis functions
- real(dp), intent(in) :: h_cubic(4)
- !> array containing the derivatives of the 4 cubic basis functions
- real(dp), intent(in) :: dh_cubic_dr(4)
-
- complex(dp), allocatable :: factors(:)
- integer, allocatable :: positions(:, :)
-
- real(dp) :: eps_inv, rho0_inv
- real(dp) :: rho0, drho0
-
- quadblock_HallB = (0.0d0, 0.0d0)
- eps_inv = 1.0d0 / eps
-
- ! Equilibrium quantities
- rho0 = rho_field % rho0(gauss_idx)
- drho0 = rho_field % d_rho0_dr(gauss_idx)
- rho0_inv = 1.0d0 / rho0
-
- ! Setup of matrix elements
-
- ! Quadratic * Quadratic
- call reset_factors(factors, 1)
- call reset_positions(positions, 1)
- ! A(6, 6)
- factors(1) = rho0_inv * ((k2 * eps_inv)**2 + k3**2)
- positions(1, :) = [6, 6]
- call subblock(quadblock_HallB, factors, positions, curr_weight, h_quadratic, h_quadratic)
-
- ! Quadratic * d(Cubic)/dr
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(6, 7)
- factors(1) = - k2 * rho0_inv * eps_inv**2
- positions(1, :) = [6, 7]
- ! A(6, 8)
- factors(2) = - k3 * rho0_inv
- positions(2, :) = [6, 8]
- call subblock(quadblock_HallB, factors, positions, curr_weight, h_quadratic, dh_cubic_dr)
-
- ! Cubic * Quadratic
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(7, 6)
- factors(1) = drho0 * k2 * eps_inv * rho0_inv**2
- positions(1, :) = [7, 6]
- ! A(8, 6)
- factors(2) = drho0 * k3 * rho0_inv**2 + d_eps_dr * k3 * eps_inv * rho0_inv
- positions(2, :) = [8, 6]
- call subblock(quadblock_HallB, factors, positions, curr_weight, h_cubic, h_quadratic)
-
- ! d(Cubic)/dr * Quadratic
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(7, 6)
- factors(1) = - k2 * eps_inv * rho0_inv
- positions(1, :) = [7, 6]
- ! A(8, 6)
- factors(2) = - k3 * rho0_inv
- positions(2, :) = [8, 6]
- call subblock(quadblock_HallB, factors, positions, curr_weight, dh_cubic_dr, h_quadratic)
-
- ! Cubic * Cubic
- call reset_factors(factors, 4)
- call reset_positions(positions, 4)
- ! A(7, 7)
- factors(1) = k3**2 * eps_inv * rho0_inv
- positions(1, :) = [7, 7]
- ! A(7, 8)
- factors(2) = - k2 * k3 * eps_inv * rho0_inv
- positions(2, :) = [7, 8]
- ! A(8, 7)
- factors(3) = - k2 * k3 * eps_inv**2 * rho0_inv
- positions(3, :) = [8, 7]
- ! A(8, 8)
- factors(4) = (k2 * eps_inv)**2 * rho0_inv
- positions(4, :) = [8, 8]
- call subblock(quadblock_HallB, factors, positions, curr_weight, h_cubic, h_cubic)
-
- ! Cubic * d(Cubic)/dr
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(7, 7)
- factors(1) = - drho0 * eps_inv * rho0_inv**2
- positions(1, :) = [7, 7]
- ! A(8, 8)
- factors(2) = - rho0_inv * (d_eps_dr * eps_inv + drho0 * rho0_inv)
- positions(2, :) = [8, 8]
- call subblock(quadblock_HallB, factors, positions, curr_weight, h_cubic, dh_cubic_dr)
-
- ! d(Cubic)/dr * d(Cubic)/dr
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(7, 7)
- factors(1) = eps_inv * rho0_inv
- positions(1, :) = [7, 7]
- ! A(8, 8)
- factors(2) = rho0_inv
- positions(2, :) = [8, 8]
- call subblock(quadblock_HallB, factors, positions, curr_weight, dh_cubic_dr, dh_cubic_dr)
-
- deallocate(factors)
- deallocate(positions)
- end subroutine get_hallterm_Bmat
-
- !> Creates a quadblock for the A matrix containing the natural boundary
- !! conditions coming from the Hall term, depending on the supplied edge.
- !! @note Currently, no boundary conditions are applied due to the use of a
- !! dropoff profile for the Hall parameter, which is zero at the
- !! edges. @endnote
- subroutine hall_boundaries(quadblock_Hall, kappa_perp_is_zero, edge)
- use mod_global_variables, only: gauss_gridpts, dim_subblock, elec_pressure
- use mod_logging, only: log_message
- use mod_equilibrium, only: rho_field, B_field, T_field
- use mod_grid, only: eps_grid
- use mod_equilibrium_params, only: k2, k3
-
- !> the quadblock corresponding to the left/right edge
- complex(dp), intent(out) :: quadblock_Hall(dim_quadblock, dim_quadblock)
- !> the edge, either "l_edge" or "r_edge"
- character(len=6), intent(in) :: edge
- !> boolean indicating if perpendicular thermal conduction is included or not
- logical, intent(in) :: kappa_perp_is_zero
-
- complex(dp), allocatable :: surface_terms(:)
- real(dp) :: eps, eps_inv, rho0_inv
- real(dp) :: rho0, B02, B03, T0
- integer, allocatable :: positions(:, :)
- integer :: idx, i
-
- quadblock_Hall = (0.0d0, 0.0d0)
-
- if (edge == 'l_edge') then
- idx = 1
- else if (edge == 'r_edge') then
- idx = gauss_gridpts
- else
- call log_message('Hall boundaries: wrong edge supplied' // edge, level='error')
- end if
-
- ! retrieve variables at current edge
- eps = eps_grid(idx)
- eps_inv = 1.0d0 / eps
- rho0 = rho_field % rho0(idx)
- rho0_inv = 1.0d0 / rho0
- T0 = T_field % T0(idx)
- B02 = B_field % B02(idx)
- B03 = B_field % B03(idx)
-
- allocate(positions(4, 2))
- allocate(surface_terms(4))
-
- ! surface term for element (6, 5)
- surface_terms(1) = 0.0d0
- positions(1, :) = [6, 5]
- ! surface term for element (6, 6)
- surface_terms(2) = - rho0_inv * (k2 * B03 * eps_inv - k3 * B02)
- positions(2, :) = [6, 6]
- ! surface term for element (6, 7)
- surface_terms(3) = B03 * eps_inv * rho0_inv
- positions(3, :) = [6, 7]
- ! surface term for element (6, 8)
- surface_terms(4) = - B02 * rho0_inv
- positions(4, :) = [6, 8]
-
- ! T1 is zero at the wall if perpendicular thermal conduction is included
- if (elec_pressure .and. kappa_perp_is_zero) then
- surface_terms(2) = surface_terms(2) + eps_inv
- end if
-
- ! l_edge: add to bottom-right of 2x2 block, for top-left subblock only
- ! r_edge: add to bottom-right of 2x2 block, for bottom-right subblock only
- if (edge == 'l_edge') then
- positions = 2 * positions
- else if (edge == 'r_edge') then
- positions = 2 * positions + dim_subblock
- end if
-
- do i = 1, 2
- quadblock_Hall(positions(i, 1), positions(i, 2)) = quadblock_Hall(positions(i, 1), positions(i, 2)) + surface_terms(i)
- end do
- do i = 3, size(surface_terms)
- quadblock_Hall(positions(i, 1), positions(i, 2)) = quadblock_Hall(positions(i, 1), positions(i, 2)) + surface_terms(i)
- end do
-
- deallocate(positions)
- deallocate(surface_terms)
- end subroutine hall_boundaries
-
- !> Retrieves the B01 (constant) contributions to the A-matrix Hall
- !! contribution for a given quadblock corresponding
- !! to a particular point in the grid and a particular Gaussian weight.
- !! If hall_mhd = True, this routine is called n_gauss times
- !! for every grid interval.
- !! @note Usage of a constant B01 component is currently not implemented. This
- !! subroutine supports a future extension of the code. @endnote
- subroutine get_hallterm_B01(gauss_idx, eps, d_eps_dr, curr_weight, &
- quadblock_Hall, h_quadratic, h_cubic, dh_cubic_dr)
- use mod_global_variables, only: ic
- use mod_equilibrium_params, only: k2, k3
- use mod_equilibrium, only: rho_field, B_field
- use mod_make_subblock, only: subblock, reset_factors, reset_positions
-
- !> current index in the Gaussian grid
- integer, intent(in) :: gauss_idx
- !> current value for the scale factor epsilon
- real(dp), intent(in) :: eps
- !> current value for the derivative of epsilon
- real(dp), intent(in) :: d_eps_dr
- !> current weight in the Gaussian quadrature
- real(dp), intent(in) :: curr_weight
- !> the A-quadblock for a particular grid interval
- complex(dp), intent(inout):: quadblock_Hall(dim_quadblock, dim_quadblock)
- !> array containing the 4 quadratic basis functions
- real(dp), intent(in) :: h_quadratic(4)
- !> array containing the 4 cubic basis functions
- real(dp), intent(in) :: h_cubic(4)
- !> array containing the derivatives of the 4 cubic basis functions
- real(dp), intent(in) :: dh_cubic_dr(4)
-
- complex(dp), allocatable :: factors(:)
- integer, allocatable :: positions(:, :)
-
- real(dp) :: eps_inv, rho0_inv
- real(dp) :: rho0, B01, B02, B03
- real(dp) :: drho0, drB02, dB02_r, dB03, dB02
-
- eps_inv = 1.0d0 / eps
-
- ! Equilibrium quantities
- rho0 = rho_field % rho0(gauss_idx)
- drho0 = rho_field % d_rho0_dr(gauss_idx)
- rho0_inv = 1.0d0 / rho0
-
- ! To be uncommented when B01 is implemented
- B01 = 0.0d0 !B_field % B01(gauss_idx)
- B02 = B_field % B02(gauss_idx)
- dB02 = B_field % d_B02_dr(gauss_idx)
- B03 = B_field % B03(gauss_idx)
- dB03 = B_field % d_B03_dr(gauss_idx)
-
- ! Calculate derivatives eps*B02, B02/eps
- drB02 = d_eps_dr * B02 + eps * dB02
- dB02_r = eps_inv * dB02 - d_eps_dr * eps_inv**2 * B02
-
- ! Setup of matrix elements
-
- ! Cubic * Quadratic
- call reset_factors(factors, 4)
- call reset_positions(positions, 4)
- ! A(7, 1)
- factors(1) = - ic * eps_inv**2 * rho0_inv**2 * B01 * drB02
- positions(1, :) = [7, 1]
- ! A(7, 6)
- factors(2) = ic * k3 * B01 * rho0_inv * (d_eps_dr * eps_inv + drho0 * rho0_inv)
- positions(2, :) = [7, 6]
- ! A(8, 1)
- factors(3) = - ic * eps_inv * rho0_inv**2 * B01 * dB03
- positions(3, :) = [8, 1]
- ! A(8, 6)
- factors(4) = - ic * drho0 * rho0_inv**2 * k2 * B01 * eps_inv
- positions(4, :) = [8, 6]
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_quadratic)
-
- ! d(Cubic)/dr * Quadratic
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(7, 6)
- factors(1) = - ic * k3 * B01 * rho0_inv
- positions(1, :) = [7, 6]
- ! A(8, 6)
- factors(2) = ic * k2 * B01 * eps_inv * rho0_inv
- positions(2, :) = [8, 6]
- call subblock(quadblock_Hall, factors, positions, curr_weight, dh_cubic_dr, h_quadratic)
-
- ! Cubic * Cubic
- call reset_factors(factors, 4)
- call reset_positions(positions, 4)
- ! A(7, 7)
- factors(1) = - ic * k2 * k3 * B01 * eps_inv**2 * rho0_inv
- positions(1, :) = [7, 7]
- ! A(7, 8)
- factors(2) = ic * k2**2 * B01 * eps_inv**2 * rho0_inv
- positions(2, :) = [7, 8]
- ! A(8, 7)
- factors(3) = - ic * k3**2 * B01 * eps_inv * rho0_inv
- positions(3, :) = [8, 7]
- ! A(8, 8)
- factors(4) = ic * k2 * k3 * B01 * eps_inv * rho0_inv
- positions(4, :) = [8, 8]
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, h_cubic)
-
- ! Cubic * d(Cubic)/dr
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(7, 8)
- factors(1) = ic * B01 * rho0_inv * (d_eps_dr * eps_inv + drho0 * rho0_inv)
- positions(1, :) = [7, 8]
- ! A(8, 7)
- factors(2) = ic * B01 * eps_inv * drho0 * rho0_inv**2
- positions(2, :) = [8, 7]
- call subblock(quadblock_Hall, factors, positions, curr_weight, h_cubic, dh_cubic_dr)
-
- ! d(Cubic)/dr * d(Cubic)/dr
- call reset_factors(factors, 2)
- call reset_positions(positions, 2)
- ! A(7, 8)
- factors(1) = - ic * B01 * rho0_inv
- positions(1, :) = [7, 8]
- ! A(8, 7)
- factors(2) = - ic * B01 * eps_inv * rho0_inv
- positions(2, :) = [8, 7]
- call subblock(quadblock_Hall, factors, positions, curr_weight, dh_cubic_dr, dh_cubic_dr)
-
- deallocate(factors)
- deallocate(positions)
- end subroutine get_hallterm_B01
-
-end module mod_hallmhd
diff --git a/tests/pylbo_tests/test_visualisations_spectra.py b/tests/pylbo_tests/test_visualisations_spectra.py
index 20b5990a..d9d29c64 100644
--- a/tests/pylbo_tests/test_visualisations_spectra.py
+++ b/tests/pylbo_tests/test_visualisations_spectra.py
@@ -33,6 +33,16 @@ def test_spectrum_plot_eigenfunctions(ds_v112):
assert getattr(p, "ef_ax", None) is not None
+def test_spectrum_plot_continua_and_eigenfunctions(ds_v112):
+ p = pylbo.plot_spectrum(ds_v112)
+ p.add_continua()
+ p.add_eigenfunctions()
+ labels = p.ax.get_legend_handles_labels()[1]
+ assert len(labels) > 2
+ for attr in ("c_handler", "ef_handler", "ef_ax"):
+ assert getattr(p, attr, None) is not None
+
+
def test_spectrum_plot_eigenfunctions_notpresent(ds_v100):
p = pylbo.plot_spectrum(ds_v100)
with pytest.raises(EigenfunctionsNotPresent):