diff --git a/post_processing/pylbo/automation/defaults.py b/post_processing/pylbo/automation/defaults.py index f3be3c33..03f2defd 100644 --- a/post_processing/pylbo/automation/defaults.py +++ b/post_processing/pylbo/automation/defaults.py @@ -56,6 +56,7 @@ ("viscosity", bool), ("viscosity_value", (int, np.integer, float)), ("hall_mhd", bool), + ("hall_substitution", bool), ("hall_dropoff", bool), ("elec_inertia", bool), ("inertia_dropoff", bool), diff --git a/src/boundaries/smod_natural_bounds_hall.f08 b/src/boundaries/smod_natural_bounds_hall.f08 index eaeda08a..703262d4 100644 --- a/src/boundaries/smod_natural_bounds_hall.f08 +++ b/src/boundaries/smod_natural_bounds_hall.f08 @@ -42,53 +42,105 @@ end procedure add_natural_hall_Bterms module procedure add_natural_hall_terms - use mod_global_variables, only: hall_mhd, viscosity, viscosity_value + use mod_global_variables, only: hall_mhd, hall_substitution, electron_fraction, & + viscosity, viscosity_value use mod_equilibrium, only: hall_field real(dp) :: eps, deps - real(dp) :: rho - real(dp) :: eta_H, mu + real(dp) :: rho, T0, B01, B02, B03 + real(dp) :: eta_H, mu, efrac - if (.not. (hall_mhd .and. viscosity)) then + if (.not. hall_mhd) then return end if eps = eps_grid(grid_idx) deps = d_eps_grid_dr(grid_idx) + rho = rho_field % rho0(grid_idx) + T0 = T_field % T0(grid_idx) + B01 = B_field % B01 + B02 = B_field % B02(grid_idx) + B03 = B_field % B03(grid_idx) + eta_H = hall_field % hallfactor(grid_idx) mu = viscosity_value - - ! ==================== Quadratic * Cubic ==================== - call reset_factor_positions(new_size=1) - ! H(6, 2) - factors(1) = -eta_H * ic * mu * deps / (eps * rho) - positions(1, :) = [6, 2] - call subblock(quadblock, factors, positions, weight, h_quad, h_cubic) - - ! ==================== Quadratic * dCubic ==================== - call reset_factor_positions(new_size=1) - ! H(6, 2) - factors(1) = 4.0d0 * eta_H * ic * mu / (3.0d0 * rho) - positions(1, :) = [6, 2] - call subblock(quadblock, factors, positions, weight, h_quad, dh_cubic) - - ! ==================== Cubic * dQuadratic ==================== - call reset_factor_positions(new_size=2) - ! H(7, 3) - factors(1) = eta_H * ic * mu * eps / rho - positions(1, :) = [7, 3] - ! H(8, 4) - factors(2) = eta_H * ic * mu / rho - positions(2, :) = [8, 4] - call subblock(quadblock, factors, positions, weight, h_cubic, dh_quad) - - ! ==================== Cubic * Quadratic ==================== - call reset_factor_positions(new_size=1) - ! H(8, 4) - factors(1) = -eta_H * ic * mu * deps / (eps * rho) - positions(1, :) = [8, 4] - call subblock(quadblock, factors, positions, weight, h_cubic, h_quad) + efrac = electron_fraction + + ! Hall by substitution of the momentum equation + if (hall_substitution) then + if (viscosity) then + ! ==================== Quadratic * Cubic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 2) + factors(1) = -eta_H * ic * mu * deps / (eps * rho) + positions(1, :) = [6, 2] + call subblock(quadblock, factors, positions, weight, h_quad, h_cubic) + + ! ==================== Quadratic * dCubic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 2) + factors(1) = 4.0d0 * eta_H * ic * mu / (3.0d0 * rho) + positions(1, :) = [6, 2] + call subblock(quadblock, factors, positions, weight, h_quad, dh_cubic) + + ! ==================== Cubic * dQuadratic ==================== + call reset_factor_positions(new_size=2) + ! H(7, 3) + factors(1) = eta_H * ic * mu * eps / rho + positions(1, :) = [7, 3] + ! H(8, 4) + factors(2) = eta_H * ic * mu / rho + positions(2, :) = [8, 4] + call subblock(quadblock, factors, positions, weight, h_cubic, dh_quad) + + ! ==================== Cubic * Quadratic ==================== + call reset_factor_positions(new_size=1) + ! H(8, 4) + factors(1) = -eta_H * ic * mu * deps / (eps * rho) + positions(1, :) = [8, 4] + call subblock(quadblock, factors, positions, weight, h_cubic, h_quad) + end if + + ! Hall without substitution, only E redefinition up to a gradient + else + ! ==================== Quadratic * Quadratic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 6) + factors(1) = eta_H * (k2 * B03 - eps * k3 * B02) / rho + positions(1, :) = [6, 6] + call subblock(quadblock, factors, positions, weight, h_quad, h_quad) + + ! ==================== Quadratic * dCubic ==================== + call reset_factor_positions(new_size=2) + ! H(6, 7) + factors(1) = -eta_H * B03 / rho + positions(1, :) = [6, 7] + ! H(6, 8) + factors(2) = eta_H * eps * B02 / rho + positions(2, :) = [6, 8] + call subblock(quadblock, factors, positions, weight, h_quad, dh_cubic) + + ! ==================== Cubic * Quadratic ==================== + call reset_factor_positions(new_size=2) + ! H(7, 6) + factors(1) = -eta_H * ic * B01 * eps * k3 / rho + positions(1, :) = [7, 6] + ! H(8, 6) + factors(2) = eta_H * ic * B01 * k2 / rho + positions(2, :) = [8, 6] + call subblock(quadblock, factors, positions, weight, h_cubic, h_quad) + + ! ==================== Cubic * dCubic ==================== + call reset_factor_positions(new_size=2) + ! H(7, 8) + factors(1) = eta_H * ic * B01 * eps / rho + positions(1, :) = [7, 8] + ! H(8, 7) + factors(2) = -eta_H * ic * B01 / rho + positions(2, :) = [8, 7] + call subblock(quadblock, factors, positions, weight, h_cubic, dh_cubic) + end if end procedure add_natural_hall_terms diff --git a/src/dataIO/mod_input.f08 b/src/dataIO/mod_input.f08 index 0c655fca..d72cc734 100644 --- a/src/dataIO/mod_input.f08 +++ b/src/dataIO/mod_input.f08 @@ -48,7 +48,8 @@ subroutine read_parfile(parfile) resistivity, use_fixed_resistivity, fixed_eta_value, & use_eta_dropoff, dropoff_edge_dist, dropoff_width, & viscosity, viscous_heating, viscosity_value, incompressible, & - hall_mhd, hall_dropoff, elec_inertia, inertia_dropoff, electron_fraction + hall_mhd, hall_substitution, hall_dropoff, & + elec_inertia, inertia_dropoff, electron_fraction namelist /unitslist/ & cgs_units, unit_density, unit_temperature, unit_magneticfield, unit_length, & mean_molecular_weight diff --git a/src/dataIO/mod_logging.f08 b/src/dataIO/mod_logging.f08 index ba93ed94..acbe112b 100644 --- a/src/dataIO/mod_logging.f08 +++ b/src/dataIO/mod_logging.f08 @@ -223,6 +223,7 @@ subroutine print_console_info() call log_message("hall mhd : " // str(hall_mhd)) if (hall_mhd) then + call log_message(" by substitution : " // str(hall_substitution)) call log_message(" electron fraction : " // str(electron_fraction)) call log_message(" electron inertia : " // str(elec_inertia)) end if diff --git a/src/equilibria/smod_equil_couette_flow.f08 b/src/equilibria/smod_equil_couette_flow.f08 index 49dd9b66..3e61ed43 100644 --- a/src/equilibria/smod_equil_couette_flow.f08 +++ b/src/equilibria/smod_equil_couette_flow.f08 @@ -1,5 +1,18 @@ ! ============================================================================= -!> This submodule defines a steady plane Couette flow in a Cartesian geometry. +!> This submodule defines a steady plane Couette flow in a Cartesian geometry +!! with flow and viscosity. +!! @note Default values are given by +!! +!! - k2 = 0 +!! - k3 = 1 +!! - cte_rho0 = 1 +!! - cte_T0 = 1 +!! - cte_v02 = 0 +!! - cte_v03 = 1 +!! - viscosity = True +!! - viscosity_value = 1e-3 +!! +!! and can all be changed in the parfile. @endnote submodule (mod_equilibrium) smod_equil_couette_flow implicit none diff --git a/src/equilibria/smod_equil_harris_sheet.f08 b/src/equilibria/smod_equil_harris_sheet.f08 index e852aef9..f24ee596 100644 --- a/src/equilibria/smod_equil_harris_sheet.f08 +++ b/src/equilibria/smod_equil_harris_sheet.f08 @@ -3,7 +3,20 @@ !! by a Harris sheet. !! !! This equilibrium is taken from Shi et al. (2020), Oblique tearing mode -!! instability: guide field and Hall effect +!! instability: guide field and Hall effect. _The Astrophysical Journal_, 902:142. +!! [DOI](https://doi.org/10.3847/1538-4357/abb6fa). +!! @note Default values are given by +!! +!! - k2 = 0.12 +!! - k3 = 0 +!! - cte_rho0 = 1 +!! - cte_T0 = 1 +!! - cte_B02 = 1 +!! - cte_B03 = 0 : guide field parameter. +!! - alpha = 1 : used to set the width of the current sheet. +!! - eq_bool = False : if True, an alternative force-free Harris sheet is used. +!! +!! and can all be changed in the parfile. @endnote submodule (mod_equilibrium) smod_equil_harris_sheet implicit none diff --git a/src/main.f08 b/src/main.f08 index cbecece6..7f710903 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, hall_field use mod_logging, only: print_logo - use mod_global_variables, only: hall_mhd, x_end, x_start + use mod_global_variables, only: hall_mhd, hall_substitution, elec_inertia, x_end, x_start character(len=str_len) :: parfile integer :: nb_evs @@ -103,6 +103,16 @@ subroutine initialisation() end if end if + if (hall_mhd .and. solver == "arnoldi") then + if (elec_inertia) then + call log_message("electron inertia not compatible with Arnoldi solver", level="error") + else if (hall_substitution) then + call log_message("Hall by substitution not compatible with Arnoldi solver, & + & changing to direct computation", level="warning") + hall_substitution = .false. + end if + end if + ! Arnoldi solver needs this, since it always calculates an orthonormal basis if (write_eigenfunctions .or. solver == "arnoldi") then call log_message("allocating eigenvector arrays", level="debug") diff --git a/src/matrices/smod_hall_matrix.f08 b/src/matrices/smod_hall_matrix.f08 index 5a827cfb..b96132d1 100644 --- a/src/matrices/smod_hall_matrix.f08 +++ b/src/matrices/smod_hall_matrix.f08 @@ -5,7 +5,7 @@ contains module procedure add_hall_bmatrix_terms - use mod_global_variables, only: elec_inertia + use mod_global_variables, only: elec_inertia, hall_substitution real(dp) :: eps, deps real(dp) :: rho, drho @@ -17,28 +17,29 @@ rho = rho_field % rho0(gauss_idx) drho = rho_field % d_rho0_dr(gauss_idx) eta_H = hall_field % hallfactor(gauss_idx) + eta_e = hall_field % inertiafactor(gauss_idx) WVop = get_wv_operator(gauss_idx) - ! ==================== Quadratic * Cubic ==================== - call reset_factor_positions(new_size=1) - ! B_H(6, 2) - factors(1) = eta_H - positions(1, :) = [6, 2] - call subblock(quadblock, factors, positions, current_weight, h_quad, h_cubic) - - ! ==================== Cubic * Quadratic ==================== - call reset_factor_positions(new_size=2) - ! B_H(7, 3) - factors(1) = eta_H * eps - positions(1, :) = [7, 3] - ! B_H(8, 4) - factors(2) = eta_H - positions(2, :) = [8, 4] - call subblock(quadblock, factors, positions, current_weight, h_cubic, h_quad) + if (hall_substitution) then + ! ==================== Quadratic * Cubic ==================== + call reset_factor_positions(new_size=1) + ! B_H(6, 2) + factors(1) = eta_H + positions(1, :) = [6, 2] + call subblock(quadblock, factors, positions, current_weight, h_quad, h_cubic) - if (elec_inertia) then - eta_e = hall_field % inertiafactor(gauss_idx) + ! ==================== Cubic * Quadratic ==================== + call reset_factor_positions(new_size=2) + ! B_H(7, 3) + factors(1) = eta_H * eps + positions(1, :) = [7, 3] + ! B_H(8, 4) + factors(2) = eta_H + positions(2, :) = [8, 4] + call subblock(quadblock, factors, positions, current_weight, h_cubic, h_quad) + end if + if (elec_inertia) then ! ==================== Quadratic * Quadratic ==================== call reset_factor_positions(new_size=1) ! B_H(6, 6) @@ -117,12 +118,14 @@ module procedure add_hall_matrix_terms - use mod_global_variables, only: viscosity, viscosity_value, electron_fraction - use mod_equilibrium, only: v_field, rho_field, T_field + use mod_global_variables, only: viscosity, viscosity_value, & + hall_substitution, electron_fraction + use mod_equilibrium, only: v_field, rho_field, T_field, B_field real(dp) :: eps, deps real(dp) :: v01, v02, v03, dv01, dv02, dv03, ddv01, ddv02, ddv03 - real(dp) :: rho, drho, dT0 + real(dp) :: rho, drho, T0, dT0, B01, B02, B03, dB02, dB03 + real(dp) :: drB02, dB02_r, Fop_plus, Gop_plus, Gop_min, WVop real(dp) :: eta_H, mu, efrac eps = eps_grid(gauss_idx) @@ -140,172 +143,311 @@ rho = rho_field % rho0(gauss_idx) drho = rho_field % d_rho0_dr(gauss_idx) + T0 = T_field % T0(gauss_idx) dT0 = T_field % d_T0_dr(gauss_idx) + B01 = B_field % B01 + 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 = deps * B02 + eps * dB02 + dB02_r = dB02 / eps - B02 * deps / eps**2 + + Fop_plus = get_F_operator(gauss_idx, which="plus") + Gop_plus = get_G_operator(gauss_idx, which="plus") + Gop_min = get_G_operator(gauss_idx, which="minus") + WVop = get_wv_operator(gauss_idx) + eta_H = hall_field % hallfactor(gauss_idx) mu = viscosity_value efrac = electron_fraction - ! ==================== Quadratic * Cubic ==================== - call reset_factor_positions(new_size=1) - ! H(6, 2) - factors(1) = eta_H * (k2 * v02 / eps + k3 * v03) - positions(1, :) = [6, 2] - call subblock(quadblock, factors, positions, current_weight, h_quad, h_cubic) - - ! ==================== Quadratic * Quadratic ==================== - call reset_factor_positions(new_size=3) - ! H(6, 1) - factors(1) = -eta_H * (1.0d0-efrac) * dT0 / rho - positions(1, :) = [6, 1] - ! H(6, 3) - factors(2) = -2.0d0 * eta_H * deps * v02 - positions(2, :) = [6, 3] - ! H(6, 5) - factors(3) = eta_H * (1.0d0-efrac) * drho / rho - positions(3, :) = [6, 5] - call subblock(quadblock, factors, positions, current_weight, h_quad, h_quad) - - ! ==================== Cubic * Cubic ==================== - call reset_factor_positions(new_size=2) - ! H(7, 2) - factors(1) = -eta_H * (dv02 - v02 * deps / eps) - positions(1, :) = [7, 2] - ! H(8, 2) - factors(2) = -eta_H * dv03 - positions(2, :) = [8, 2] - call subblock(quadblock, factors, positions, current_weight, h_cubic, h_cubic) - - ! ==================== Cubic * Quadratic ==================== - call reset_factor_positions(new_size=2) - ! H(7, 3) - factors(1) = eta_H * (k2 * v02 + eps * k3 * v03) - positions(1, :) = [7, 3] - ! H(8, 4) - factors(2) = eta_H * (k2 * v02 / eps + k3 * v03) - positions(2, :) = [8, 4] - call subblock(quadblock, factors, positions, current_weight, h_cubic, h_quad) - - if (viscosity) then + ! Hall by substitution of the momentum equation + if (hall_substitution) then ! ==================== Quadratic * Cubic ==================== call reset_factor_positions(new_size=1) ! H(6, 2) - factors(1) = -eta_H * ic * mu * ((drho / rho + 1.0d0 / eps) * deps / eps & - + (k2 / eps)**2 + k3**2) / rho + factors(1) = eta_H * (k2 * v02 / eps + k3 * v03) positions(1, :) = [6, 2] call subblock(quadblock, factors, positions, current_weight, h_quad, h_cubic) - ! ==================== Quadratic * dCubic ==================== - call reset_factor_positions(new_size=1) - ! H(6, 2) - factors(1) = eta_H * ic * mu * (4.0d0 * drho / rho - deps / eps) / (3.0d0 * rho) - positions(1, :) = [6, 2] - call subblock(quadblock, factors, positions, current_weight, h_quad, dh_cubic) - - ! ==================== dQuadratic * Cubic ==================== - call reset_factor_positions(new_size=1) - ! H(6, 2) - factors(1) = eta_H * ic * mu * deps / (eps * rho) - positions(1, :) = [6, 2] - call subblock(quadblock, factors, positions, current_weight, dh_quad, h_cubic) - - ! ==================== dQuadratic * dCubic ==================== - call reset_factor_positions(new_size=1) - ! H(6, 2) - factors(1) = -4.0d0 * eta_H * ic * mu / (3.0d0 * rho) - positions(1, :) = [6, 2] - call subblock(quadblock, factors, positions, current_weight, dh_quad, dh_cubic) - ! ==================== Quadratic * Quadratic ==================== call reset_factor_positions(new_size=3) ! H(6, 1) - factors(1) = eta_H * mu * 4.0d0 * (ddv01 + deps * (dv01 - v01 / eps) / eps) / (3.0d0 * rho**2) + factors(1) = -eta_H * (1.0d0-efrac) * dT0 / rho positions(1, :) = [6, 1] ! H(6, 3) - factors(2) = 7.0d0 * eta_H * ic * mu * deps * k2 / (3.0d0 * eps * rho) + factors(2) = -2.0d0 * eta_H * deps * v02 positions(2, :) = [6, 3] - ! H(6, 4) - factors(3) = eta_H * ic * mu * k3 * deps / (3.0d0 * eps * rho) - positions(3, :) = [6, 4] + ! H(6, 5) + factors(3) = eta_H * (1.0d0-efrac) * drho / rho + positions(3, :) = [6, 5] call subblock(quadblock, factors, positions, current_weight, h_quad, h_quad) - ! ==================== Quadratic * dQuadratic ==================== - call reset_factor_positions(new_size=2) - ! H(6, 3) - factors(1) = -eta_H * ic * mu * k2 / (3.0d0 * rho) - positions(1, :) = [6, 3] - ! H(6, 4) - factors(2) = -eta_H * ic * mu * k3 / (3.0d0 * rho) - positions(2, :) = [6, 4] - call subblock(quadblock, factors, positions, current_weight, h_quad, dh_quad) - ! ==================== Cubic * Cubic ==================== - call reset_factor_positions(new_size=1) - ! H(7, 2) - factors(1) = 2.0d0 * eta_H * ic * mu * k2 * deps / (eps**2 * rho) - positions(1, :) = [7, 2] - call subblock(quadblock, factors, positions, current_weight, h_cubic, h_cubic) - - ! ==================== Cubic * dCubic ==================== call reset_factor_positions(new_size=2) ! H(7, 2) - factors(1) = eta_H * ic * mu * k2 / (3.0d0 * eps * rho) + factors(1) = -eta_H * (dv02 - v02 * deps / eps) positions(1, :) = [7, 2] ! H(8, 2) - factors(2) = eta_H * ic * mu * k3 / (3.0d0 * rho) + factors(2) = -eta_H * dv03 positions(2, :) = [8, 2] - call subblock(quadblock, factors, positions, current_weight, h_cubic, dh_cubic) + call subblock(quadblock, factors, positions, current_weight, h_cubic, h_cubic) ! ==================== Cubic * Quadratic ==================== - call reset_factor_positions(new_size=6) - ! H(7, 1) - factors(1) = -ic * eta_H * mu * (ddv02 + deps * (dv02 - v02 / eps) / eps) / rho**2 - positions(1, :) = [7, 1] - ! H(7, 3) - factors(2) = -eta_H * ic * mu * (4.0d0 * k2**2 / (3.0d0 * eps) & - + eps * k3**2 + deps / eps) / rho - positions(2, :) = [7, 3] - ! H(7, 4) - factors(3) = -eta_H * ic * mu * k2 * k3 / (3.0d0 * eps * rho) - positions(3, :) = [7, 4] - ! H(8, 1) - factors(4) = -ic * eta_H * mu * (ddv03 + deps * dv03 / eps) / rho**2 - positions(4, :) = [8, 1] - ! H(8, 3) - factors(5) = -eta_H * ic * mu * k2 * k3 / (3.0d0 * rho) - positions(5, :) = [8, 3] - ! H(8, 4) - factors(6) = -eta_H * ic * mu * ((k2 / eps)**2 + 4.0d0 * k3**2 / 3.0d0 & - + drho * deps / (eps * rho)) / rho - positions(6, :) = [8, 4] - call subblock(quadblock, factors, positions, current_weight, h_cubic, h_quad) - - ! ==================== Cubic * dQuadratic ==================== call reset_factor_positions(new_size=2) ! H(7, 3) - factors(1) = eta_H * ic * mu * eps * drho / rho**2 + factors(1) = eta_H * (k2 * v02 + eps * k3 * v03) positions(1, :) = [7, 3] ! H(8, 4) - factors(2) = eta_H * ic * mu * drho / rho**2 + factors(2) = eta_H * (k2 * v02 / eps + k3 * v03) positions(2, :) = [8, 4] - call subblock(quadblock, factors, positions, current_weight, h_cubic, dh_quad) + call subblock(quadblock, factors, positions, current_weight, h_cubic, h_quad) + + if (viscosity) then + ! ==================== Quadratic * Cubic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 2) + factors(1) = -eta_H * ic * mu * ((drho / rho + 1.0d0 / eps) * deps / eps & + + (k2 / eps)**2 + k3**2) / rho + positions(1, :) = [6, 2] + call subblock(quadblock, factors, positions, current_weight, h_quad, h_cubic) + + ! ==================== Quadratic * dCubic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 2) + factors(1) = eta_H * ic * mu * (4.0d0 * drho / rho - deps / eps) / (3.0d0 * rho) + positions(1, :) = [6, 2] + call subblock(quadblock, factors, positions, current_weight, h_quad, dh_cubic) + + ! ==================== dQuadratic * Cubic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 2) + factors(1) = eta_H * ic * mu * deps / (eps * rho) + positions(1, :) = [6, 2] + call subblock(quadblock, factors, positions, current_weight, dh_quad, h_cubic) + + ! ==================== dQuadratic * dCubic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 2) + factors(1) = -4.0d0 * eta_H * ic * mu / (3.0d0 * rho) + positions(1, :) = [6, 2] + call subblock(quadblock, factors, positions, current_weight, dh_quad, dh_cubic) + + ! ==================== Quadratic * Quadratic ==================== + call reset_factor_positions(new_size=3) + ! H(6, 1) + factors(1) = eta_H * mu * 4.0d0 * (ddv01 + deps * (dv01 - v01 / eps) / eps) / (3.0d0 * rho**2) + positions(1, :) = [6, 1] + ! H(6, 3) + factors(2) = 7.0d0 * eta_H * ic * mu * deps * k2 / (3.0d0 * eps * rho) + positions(2, :) = [6, 3] + ! H(6, 4) + factors(3) = eta_H * ic * mu * k3 * deps / (3.0d0 * eps * rho) + positions(3, :) = [6, 4] + call subblock(quadblock, factors, positions, current_weight, h_quad, h_quad) + + ! ==================== Quadratic * dQuadratic ==================== + call reset_factor_positions(new_size=2) + ! H(6, 3) + factors(1) = -eta_H * ic * mu * k2 / (3.0d0 * rho) + positions(1, :) = [6, 3] + ! H(6, 4) + factors(2) = -eta_H * ic * mu * k3 / (3.0d0 * rho) + positions(2, :) = [6, 4] + call subblock(quadblock, factors, positions, current_weight, h_quad, dh_quad) + + ! ==================== Cubic * Cubic ==================== + call reset_factor_positions(new_size=1) + ! H(7, 2) + factors(1) = 2.0d0 * eta_H * ic * mu * k2 * deps / (eps**2 * rho) + positions(1, :) = [7, 2] + call subblock(quadblock, factors, positions, current_weight, h_cubic, h_cubic) + + ! ==================== Cubic * dCubic ==================== + call reset_factor_positions(new_size=2) + ! H(7, 2) + factors(1) = eta_H * ic * mu * k2 / (3.0d0 * eps * rho) + positions(1, :) = [7, 2] + ! H(8, 2) + factors(2) = eta_H * ic * mu * k3 / (3.0d0 * rho) + positions(2, :) = [8, 2] + call subblock(quadblock, factors, positions, current_weight, h_cubic, dh_cubic) + + ! ==================== Cubic * Quadratic ==================== + call reset_factor_positions(new_size=6) + ! H(7, 1) + factors(1) = -ic * eta_H * mu * (ddv02 + deps * (dv02 - v02 / eps) / eps) / rho**2 + positions(1, :) = [7, 1] + ! H(7, 3) + factors(2) = -eta_H * ic * mu * (4.0d0 * k2**2 / (3.0d0 * eps) & + + eps * k3**2 + deps / eps) / rho + positions(2, :) = [7, 3] + ! H(7, 4) + factors(3) = -eta_H * ic * mu * k2 * k3 / (3.0d0 * eps * rho) + positions(3, :) = [7, 4] + ! H(8, 1) + factors(4) = -ic * eta_H * mu * (ddv03 + deps * dv03 / eps) / rho**2 + positions(4, :) = [8, 1] + ! H(8, 3) + factors(5) = -eta_H * ic * mu * k2 * k3 / (3.0d0 * rho) + positions(5, :) = [8, 3] + ! H(8, 4) + factors(6) = -eta_H * ic * mu * ((k2 / eps)**2 + 4.0d0 * k3**2 / 3.0d0 & + + drho * deps / (eps * rho)) / rho + positions(6, :) = [8, 4] + call subblock(quadblock, factors, positions, current_weight, h_cubic, h_quad) + + ! ==================== Cubic * dQuadratic ==================== + call reset_factor_positions(new_size=2) + ! H(7, 3) + factors(1) = eta_H * ic * mu * eps * drho / rho**2 + positions(1, :) = [7, 3] + ! H(8, 4) + factors(2) = eta_H * ic * mu * drho / rho**2 + positions(2, :) = [8, 4] + call subblock(quadblock, factors, positions, current_weight, h_cubic, dh_quad) + + ! ==================== dCubic * dQuadratic ==================== + call reset_factor_positions(new_size=2) + ! H(7, 3) + factors(1) = -eta_H * ic * mu * eps / rho + positions(1, :) = [7, 3] + ! H(8, 4) + factors(2) = -eta_H * ic * mu / rho + positions(2, :) = [8, 4] + call subblock(quadblock, factors, positions, current_weight, dh_cubic, dh_quad) + + ! ==================== dCubic * Quadratic ==================== + call reset_factor_positions(new_size=1) + ! H(8, 4) + factors(1) = eta_H * ic * mu * deps / (eps * rho) + positions(1, :) = [8, 4] + call subblock(quadblock, factors, positions, current_weight, dh_cubic, h_quad) + end if + + ! Hall without substitution, only E redefinition up to a gradient + else + ! ==================== Quadratic * Quadratic ==================== + call reset_factor_positions(new_size=3) + ! H(6, 1) + factors(1) = eta_H * ((B02 * drB02 / eps + B03 * dB03) / rho + efrac * dT0) / rho + positions(1, :) = [6, 1] + ! H(6, 5) + factors(2) = -eta_H * efrac * drho / rho + positions(2, :) = [6, 5] + ! H(6, 6) + factors(3) = -eta_H * (eps * drho * Gop_min / rho + deps * Gop_plus) / rho + positions(3, :) = [6, 6] + call subblock(quadblock, factors, positions, current_weight, h_quad, h_quad) + + ! ==================== dQuadratic * Quadratic ==================== + call reset_factor_positions(new_size=1) + ! H(6, 6) + factors(1) = eta_H * eps * Gop_min / rho + positions(1, :) = [6, 6] + call subblock(quadblock, factors, positions, current_weight, dh_quad, h_quad) - ! ==================== dCubic * dQuadratic ==================== + ! ==================== Quadratic * Cubic ==================== call reset_factor_positions(new_size=2) - ! H(7, 3) - factors(1) = -eta_H * ic * mu * eps / rho - positions(1, :) = [7, 3] - ! H(8, 4) - factors(2) = -eta_H * ic * mu / rho - positions(2, :) = [8, 4] - call subblock(quadblock, factors, positions, current_weight, dh_cubic, dh_quad) + ! H(6, 7) + factors(1) = eta_H * k3 * Fop_plus / rho + positions(1, :) = [6, 7] + ! H(6, 8) + factors(2) = -eta_H * k2 * Fop_plus / rho + positions(2, :) = [6, 8] + call subblock(quadblock, factors, positions, current_weight, h_quad, h_cubic) + + ! ==================== Quadratic * dCubic ==================== + call reset_factor_positions(new_size=2) + ! H(6, 7) + factors(1) = eta_H * B03 * (deps / eps - drho / rho) / rho + positions(1, :) = [6, 7] + ! H(6, 8) + factors(2) = eta_H * B02 * (deps + eps * drho / rho) / rho + positions(2, :) = [6, 8] + call subblock(quadblock, factors, positions, current_weight, h_quad, dh_cubic) + + ! ==================== dQuadratic * dCubic ==================== + call reset_factor_positions(new_size=2) + ! H(6, 7) + factors(1) = eta_H * B03 / rho + positions(1, :) = [6, 7] + ! H(6, 8) + factors(2) = -eta_H * eps * B02 / rho + positions(2, :) = [6, 8] + call subblock(quadblock, factors, positions, current_weight, dh_quad, dh_cubic) + + ! ==================== Cubic * Quadratic ==================== + call reset_factor_positions(new_size=4) + ! H(7, 1) + factors(1) = eta_H * ic * B01 * drB02 / rho**2 + positions(1, :) = [7, 1] + ! H(7, 6) + factors(2) = eta_H * (WVop * B03 - ic * B01 * eps * k3 * drho / rho) / rho + positions(2, :) = [7, 6] + ! H(8, 1) + factors(3) = eta_H * ic * eps * B01 * dB03 / rho**2 + positions(3, :) = [8, 1] + ! H(8, 6) + factors(4) = -eta_H * (WVop * B02 + ic * B01 * k2 * (deps / eps - drho / rho)) / rho + positions(4, :) = [8, 6] + call subblock(quadblock, factors, positions, current_weight, h_cubic, h_quad) ! ==================== dCubic * Quadratic ==================== - call reset_factor_positions(new_size=1) - ! H(8, 4) - factors(1) = eta_H * ic * mu * deps / (eps * rho) - positions(1, :) = [8, 4] + call reset_factor_positions(new_size=2) + ! H(7, 6) + factors(1) = eta_H * ic * B01 * eps * k3 / rho + positions(1, :) = [7, 6] + ! H(8, 6) + factors(2) = -eta_H * ic * B01 * k2 / rho + positions(2, :) = [8, 6] call subblock(quadblock, factors, positions, current_weight, dh_cubic, h_quad) + + ! ==================== Cubic * Cubic ==================== + call reset_factor_positions(new_size=4) + ! H(7, 7) + factors(1) = eta_H * k3 * (ic * B01 * k2 - drB02) / (eps * rho) + positions(1, :) = [7, 7] + ! H(7, 8) + factors(2) = -eta_H * k2 * (ic * B01 * k2 - drB02) / (eps * rho) + positions(2, :) = [7, 8] + ! H(8, 7) + factors(3) = eta_H * k3 * (ic * B01 * k3 - dB03) / rho + positions(3, :) = [8, 7] + ! H(8, 8) + factors(4) = -eta_H * k2 * (ic * B01 * k3 - dB03) / rho + positions(4, :) = [8, 8] + call subblock(quadblock, factors, positions, current_weight, h_cubic, h_cubic) + + ! ==================== Cubic * dCubic ==================== + call reset_factor_positions(new_size=4) + ! H(7, 7) + factors(1) = -eta_H * k2 * B03 / (eps * rho) + positions(1, :) = [7, 7] + ! H(7, 8) + factors(2) = eta_H * eps * (ic * B01 * drho / rho - k3 * B03) / rho + positions(2, :) = [7, 8] + ! H(8, 7) + factors(3) = eta_H * (ic * B01 * (deps / eps - drho / rho) + k2 * B02 / eps) / rho + positions(3, :) = [8, 7] + ! H(8, 8) + factors(4) = eta_H * eps * k3 * B02 / rho + positions(4, :) = [8, 8] + call subblock(quadblock, factors, positions, current_weight, h_cubic, dh_cubic) + + ! ==================== dCubic * dCubic ==================== + call reset_factor_positions(new_size=2) + ! H(7, 8) + factors(1) = -eta_H * ic * B01 * eps / rho + positions(1, :) = [7, 8] + ! H(8, 7) + factors(2) = eta_H * ic * B01 / rho + positions(2, :) = [8, 7] + call subblock(quadblock, factors, positions, current_weight, dh_cubic, dh_cubic) end if end procedure add_hall_matrix_terms diff --git a/src/mod_global_variables.f08 b/src/mod_global_variables.f08 index 5feafdec..f89db59c 100644 --- a/src/mod_global_variables.f08 +++ b/src/mod_global_variables.f08 @@ -71,19 +71,21 @@ module mod_global_variables real(dp) :: dropoff_width !> boolean for viscosity, defaults to False logical, save :: viscosity - !> boolean to include viscoous heating, defaults to False + !> boolean to include viscous heating, defaults to False logical, save :: viscous_heating - !> defines the fixed value for either the dynamic or kinematic viscosity, defaults to 0 + !> defines the fixed value for the dynamic viscosity, defaults to 0 real(dp) :: viscosity_value !> boolean to use Hall MHD, defaults to False logical, save :: hall_mhd - !> boolean to use dropoff profile for Hall parameter, defaults to False + !> boolean to use substitution for Hall elements (as presented in the paper), defaults to True + logical, save :: hall_substitution + !> boolean to use dropoff profile for Hall parameter, defaults to False logical, save :: hall_dropoff !> boolean to use electron inertia in Ohm's law, defaults to False logical, save :: elec_inertia - !> boolean to use dropoff profile for inertia parameter, defaults to False + !> boolean to use dropoff profile for inertia parameter, defaults to False logical, save :: inertia_dropoff - !> fraction of number of electrons to number of all particles (between 0 and 1) + !> fraction of number of electrons to number of all particles (between 0 and 1), defaults to 0.5 real(dp) :: electron_fraction !> defines the geometry of the problem, defaults depend on chosen equilibrium character(len=str_len) :: geometry @@ -216,10 +218,11 @@ subroutine initialise_globals() viscous_heating = .false. viscosity_value = 0.0d0 hall_mhd = .false. + hall_substitution = .true. hall_dropoff = .false. elec_inertia = .false. inertia_dropoff = .false. - electron_fraction = 0.0d0 + electron_fraction = 0.5d0 !! grid variables ! do not initialise these three so they MUST be set in submodules/parfile diff --git a/src/physics/mod_hall.f08 b/src/physics/mod_hall.f08 index 745a1d4d..0794d947 100644 --- a/src/physics/mod_hall.f08 +++ b/src/physics/mod_hall.f08 @@ -1,6 +1,6 @@ ! ============================================================================= -!> !> Module containing Hall-related routines, calculates -!! and sets the Hall and inertia term factors based on the specified profiles. +!> Module containing Hall-related routines. +!! Sets the Hall and electron inertia factors based on normalisations and specified profiles. module mod_hall use mod_global_variables, only: dp, dim_quadblock @@ -13,8 +13,8 @@ module mod_hall contains !> Retrieves the normalised Hall factor as described by Porth et al. (2014), - !! with a dropoff at the boundary, if desired. Additionally defines the - !! inertia term factor if included, with a dropoff profile, if desired. + !! with a dropoff at the boundary, if desired. Additionally, defines the + !! electron inertia factor if included, with a dropoff profile, if desired. subroutine set_hall_factors(hall_field) use mod_grid, only: grid_gauss use mod_physical_constants, only: dpi diff --git a/tests/regression_tests/answers/answer_MRI_accretion.dat b/tests/regression_tests/answers/answer_MRI_accretion.dat index 1a5576d1..b0fcf0ab 100644 Binary files a/tests/regression_tests/answers/answer_MRI_accretion.dat and b/tests/regression_tests/answers/answer_MRI_accretion.dat differ diff --git a/tests/regression_tests/answers/answer_RTI_KHI.dat b/tests/regression_tests/answers/answer_RTI_KHI.dat index 3b239596..034e9773 100644 Binary files a/tests/regression_tests/answers/answer_RTI_KHI.dat and b/tests/regression_tests/answers/answer_RTI_KHI.dat differ diff --git a/tests/regression_tests/answers/answer_RTI_theta_pinch_HD.dat b/tests/regression_tests/answers/answer_RTI_theta_pinch_HD.dat index ce8e1b6d..e23c1a04 100644 Binary files a/tests/regression_tests/answers/answer_RTI_theta_pinch_HD.dat and b/tests/regression_tests/answers/answer_RTI_theta_pinch_HD.dat differ diff --git a/tests/regression_tests/answers/answer_RTI_theta_pinch_MHD.dat b/tests/regression_tests/answers/answer_RTI_theta_pinch_MHD.dat index 94ba0534..c7131f85 100644 Binary files a/tests/regression_tests/answers/answer_RTI_theta_pinch_MHD.dat and b/tests/regression_tests/answers/answer_RTI_theta_pinch_MHD.dat differ diff --git a/tests/regression_tests/answers/answer_adiabatic_homo.dat b/tests/regression_tests/answers/answer_adiabatic_homo.dat index 451c6fa3..4e53c3ac 100644 Binary files a/tests/regression_tests/answers/answer_adiabatic_homo.dat and b/tests/regression_tests/answers/answer_adiabatic_homo.dat differ diff --git a/tests/regression_tests/answers/answer_adiabatic_homo_ef_subset.dat b/tests/regression_tests/answers/answer_adiabatic_homo_ef_subset.dat index bb7ea9c2..4d7db1bd 100644 Binary files a/tests/regression_tests/answers/answer_adiabatic_homo_ef_subset.dat and b/tests/regression_tests/answers/answer_adiabatic_homo_ef_subset.dat differ diff --git a/tests/regression_tests/answers/answer_couette_flow_QR.dat b/tests/regression_tests/answers/answer_couette_flow_QR.dat index 302bbdd6..fef5652a 100644 Binary files a/tests/regression_tests/answers/answer_couette_flow_QR.dat and b/tests/regression_tests/answers/answer_couette_flow_QR.dat differ diff --git a/tests/regression_tests/answers/answer_discrete_alfven.dat b/tests/regression_tests/answers/answer_discrete_alfven.dat index b352b061..67942f05 100644 Binary files a/tests/regression_tests/answers/answer_discrete_alfven.dat and b/tests/regression_tests/answers/answer_discrete_alfven.dat differ diff --git a/tests/regression_tests/answers/answer_elecinertia_multi.pickle b/tests/regression_tests/answers/answer_elecinertia_multi.pickle new file mode 100644 index 00000000..0088d1e3 Binary files /dev/null and b/tests/regression_tests/answers/answer_elecinertia_multi.pickle differ diff --git a/tests/regression_tests/answers/answer_fluxtube_coronal.dat b/tests/regression_tests/answers/answer_fluxtube_coronal.dat index 113ff5a9..1e0802f6 100644 Binary files a/tests/regression_tests/answers/answer_fluxtube_coronal.dat and b/tests/regression_tests/answers/answer_fluxtube_coronal.dat differ diff --git a/tests/regression_tests/answers/answer_fluxtube_photospheric.dat b/tests/regression_tests/answers/answer_fluxtube_photospheric.dat index a8a341fa..dab9f6d0 100644 Binary files a/tests/regression_tests/answers/answer_fluxtube_photospheric.dat and b/tests/regression_tests/answers/answer_fluxtube_photospheric.dat differ diff --git a/tests/regression_tests/answers/answer_gold_hoyle.dat b/tests/regression_tests/answers/answer_gold_hoyle.dat index f0e3867d..6185ae4d 100644 Binary files a/tests/regression_tests/answers/answer_gold_hoyle.dat and b/tests/regression_tests/answers/answer_gold_hoyle.dat differ diff --git a/tests/regression_tests/answers/answer_interchange_modes.dat b/tests/regression_tests/answers/answer_interchange_modes.dat index 72b7bbe7..284d120b 100644 Binary files a/tests/regression_tests/answers/answer_interchange_modes.dat and b/tests/regression_tests/answers/answer_interchange_modes.dat differ diff --git a/tests/regression_tests/answers/answer_internal_kink.dat b/tests/regression_tests/answers/answer_internal_kink.dat index e9b8aa31..0a36434a 100644 Binary files a/tests/regression_tests/answers/answer_internal_kink.dat and b/tests/regression_tests/answers/answer_internal_kink.dat differ diff --git a/tests/regression_tests/answers/answer_kelvin_helmholtz.dat b/tests/regression_tests/answers/answer_kelvin_helmholtz.dat index 8e808006..8f968845 100644 Binary files a/tests/regression_tests/answers/answer_kelvin_helmholtz.dat and b/tests/regression_tests/answers/answer_kelvin_helmholtz.dat differ diff --git a/tests/regression_tests/answers/answer_kh_cd.dat b/tests/regression_tests/answers/answer_kh_cd.dat index 8304123b..3127900e 100644 Binary files a/tests/regression_tests/answers/answer_kh_cd.dat and b/tests/regression_tests/answers/answer_kh_cd.dat differ diff --git a/tests/regression_tests/answers/answer_magnetothermal.dat b/tests/regression_tests/answers/answer_magnetothermal.dat index 0cc23c75..558319e3 100644 Binary files a/tests/regression_tests/answers/answer_magnetothermal.dat and b/tests/regression_tests/answers/answer_magnetothermal.dat differ diff --git a/tests/regression_tests/answers/answer_magnetothermal_arnoldi_si.dat b/tests/regression_tests/answers/answer_magnetothermal_arnoldi_si.dat index 499d5105..b1a9d9a7 100644 Binary files a/tests/regression_tests/answers/answer_magnetothermal_arnoldi_si.dat and b/tests/regression_tests/answers/answer_magnetothermal_arnoldi_si.dat differ diff --git a/tests/regression_tests/answers/answer_quasimodes.dat b/tests/regression_tests/answers/answer_quasimodes.dat index 52122e7a..54a26b20 100644 Binary files a/tests/regression_tests/answers/answer_quasimodes.dat and b/tests/regression_tests/answers/answer_quasimodes.dat differ diff --git a/tests/regression_tests/answers/answer_rayleigh_taylor.dat b/tests/regression_tests/answers/answer_rayleigh_taylor.dat index 8845ba93..ebcc1231 100644 Binary files a/tests/regression_tests/answers/answer_rayleigh_taylor.dat and b/tests/regression_tests/answers/answer_rayleigh_taylor.dat differ diff --git a/tests/regression_tests/answers/answer_resistive_homo.dat b/tests/regression_tests/answers/answer_resistive_homo.dat index b3a981d0..1987d1e2 100644 Binary files a/tests/regression_tests/answers/answer_resistive_homo.dat and b/tests/regression_tests/answers/answer_resistive_homo.dat differ diff --git a/tests/regression_tests/answers/answer_resistive_homo_arnoldi_si.dat b/tests/regression_tests/answers/answer_resistive_homo_arnoldi_si.dat index 5d5af2b4..432ca43f 100644 Binary files a/tests/regression_tests/answers/answer_resistive_homo_arnoldi_si.dat and b/tests/regression_tests/answers/answer_resistive_homo_arnoldi_si.dat differ diff --git a/tests/regression_tests/answers/answer_resistive_homo_ef_subset.dat b/tests/regression_tests/answers/answer_resistive_homo_ef_subset.dat index 2bbb4e62..94071759 100644 Binary files a/tests/regression_tests/answers/answer_resistive_homo_ef_subset.dat and b/tests/regression_tests/answers/answer_resistive_homo_ef_subset.dat differ diff --git a/tests/regression_tests/answers/answer_resistive_tearing.dat b/tests/regression_tests/answers/answer_resistive_tearing.dat index a298dc46..7057536a 100644 Binary files a/tests/regression_tests/answers/answer_resistive_tearing.dat and b/tests/regression_tests/answers/answer_resistive_tearing.dat differ diff --git a/tests/regression_tests/answers/answer_resistive_tearing_flow.dat b/tests/regression_tests/answers/answer_resistive_tearing_flow.dat index 106bd73a..ad8fad94 100644 Binary files a/tests/regression_tests/answers/answer_resistive_tearing_flow.dat and b/tests/regression_tests/answers/answer_resistive_tearing_flow.dat differ diff --git a/tests/regression_tests/answers/answer_rotating_cylinder.dat b/tests/regression_tests/answers/answer_rotating_cylinder.dat index 690170ca..ff7f7169 100644 Binary files a/tests/regression_tests/answers/answer_rotating_cylinder.dat and b/tests/regression_tests/answers/answer_rotating_cylinder.dat differ diff --git a/tests/regression_tests/answers/answer_suydam.dat b/tests/regression_tests/answers/answer_suydam.dat index 97af992a..68b55a38 100644 Binary files a/tests/regression_tests/answers/answer_suydam.dat and b/tests/regression_tests/answers/answer_suydam.dat differ diff --git a/tests/regression_tests/answers/answer_taylor_couette_QR.dat b/tests/regression_tests/answers/answer_taylor_couette_QR.dat index 23b81720..5b400c0c 100644 Binary files a/tests/regression_tests/answers/answer_taylor_couette_QR.dat and b/tests/regression_tests/answers/answer_taylor_couette_QR.dat differ diff --git a/tests/regression_tests/answers/answer_tokamak.dat b/tests/regression_tests/answers/answer_tokamak.dat index 4e9f8119..3de898e0 100644 Binary files a/tests/regression_tests/answers/answer_tokamak.dat and b/tests/regression_tests/answers/answer_tokamak.dat differ diff --git a/tests/regression_tests/regression.py b/tests/regression_tests/regression.py index a74812d8..0b02bb4f 100644 --- a/tests/regression_tests/regression.py +++ b/tests/regression_tests/regression.py @@ -47,6 +47,8 @@ from regression_tests.test_suydam import suydam_setup from regression_tests.test_taylor_couette import taylor_couette_setup from regression_tests.test_tokamak import tokamak_setup +from regression_tests.test_hall_adiabatic_homo import hall_adiabatic_homo +from regression_tests.test_hall_harris_sheet import hall_harris_sheet from regression_tests.test_multi_constant_current import constant_current_setup from regression_tests.test_multi_coronal_fluxtube import coronal_tube_setup @@ -55,9 +57,7 @@ from regression_tests.test_multi_interchange import interchange_setup from regression_tests.test_multi_iso_atmo import iso_atmo_beta_half_setup from regression_tests.test_multi_photospheric_fluxtube import photospheric_tube_setup - -from regression_tests.test_hall_adiabatic_homo import hall_adiabatic_homo -from regression_tests.test_hall_harris_sheet import hall_harris_sheet +from regression_tests.test_multi_elecinertia import elecinertia_setup # retrieve test setups tests_to_run = [ @@ -100,6 +100,7 @@ interchange_setup, iso_atmo_beta_half_setup, photospheric_tube_setup, + elecinertia_setup, ] # configure test setup @@ -331,6 +332,10 @@ def test_generate_multi_spectra_images(series_test, series_answer, setup, imaged p.set_y_scaling(setup.get("y_scaling", 1)) if setup.get("symlog", None) is not None: p.ax.set_yscale("symlog", linthresh=setup["symlog"]) + if setup.get("xlog"): + p.ax.set_xscale("log") + if setup.get("ylog"): + p.ax.set_yscale("log") p.ax.set_xlim(xlims) p.ax.set_ylim(ylims) figname = f"{get_image_filename(setup['name'], limits_dict=setup['limits'])}.png" diff --git a/tests/regression_tests/test_hall_adiabatic_homo.py b/tests/regression_tests/test_hall_adiabatic_homo.py index 27235286..21df0bfd 100644 --- a/tests/regression_tests/test_hall_adiabatic_homo.py +++ b/tests/regression_tests/test_hall_adiabatic_homo.py @@ -17,6 +17,7 @@ }, "equilibrium_type": "adiabatic_homo", "hall_mhd": True, + "hall_substitution": True, "electron_fraction": 0.5, "cgs_units": True, "unit_density": 1.7e-14, diff --git a/tests/regression_tests/test_hall_harris_sheet.py b/tests/regression_tests/test_hall_harris_sheet.py index 65fc3815..586a4bda 100644 --- a/tests/regression_tests/test_hall_harris_sheet.py +++ b/tests/regression_tests/test_hall_harris_sheet.py @@ -19,6 +19,7 @@ "use_fixed_resistivity": True, "fixed_eta_value": 10 ** (-4.0), "hall_mhd": True, + "hall_substitution": True, "electron_fraction": 0.5, "cgs_units": True, "unit_density": 1.7e-14, diff --git a/tests/regression_tests/test_multi_elecinertia.py b/tests/regression_tests/test_multi_elecinertia.py new file mode 100644 index 00000000..50064dfd --- /dev/null +++ b/tests/regression_tests/test_multi_elecinertia.py @@ -0,0 +1,41 @@ +import numpy as np + +NB_RUNS = 24 +_kvals = np.logspace(-1, 4, NB_RUNS) + +elecinertia_setup = { + "name": "elecinertia_multi", + "config": { + "equilibrium_type": "adiabatic_homo", + "gridpoints": 51, + "number_of_runs": NB_RUNS, + "geometry": "Cartesian", + "x_start": 0, + "x_end": 1000, + "parameters": { + "k2": _kvals * np.sin(np.pi / 6), + "k3": _kvals * np.cos(np.pi / 6), + "cte_rho0": 1, + "cte_T0": 1, + "cte_B02": 0, + "cte_B03": 1, + }, + "hall_mhd": True, + "hall_substitution": True, + "electron_fraction": 0.5, + "elec_inertia": True, + "cgs_units": True, + "unit_density": 1.7e-14, + "unit_magneticfield": 10, + "unit_length": 7.534209349981049e-9, + "basename_datfile": "elecinertia", + "write_eigenfunctions": False, + "show_results": False, + "logging_level": 0, + }, + "xdata": _kvals, + "use_squared_omega": False, + "xlog": True, + "ylog": True, + "limits": {"xlims": (1e-1, 1e4), "ylims": (1e-2, 1e4)}, +}