Alternative Hall implementation and documentation update

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),

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
 

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

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

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
+!!
+!! - <tt>k2</tt> = 0
+!! - <tt>k3</tt> = 1
+!! - <tt>cte_rho0</tt> = 1
+!! - <tt>cte_T0</tt> = 1
+!! - <tt>cte_v02</tt> = 0
+!! - <tt>cte_v03</tt> = 1
+!! - <tt>viscosity</tt> = True
+!! - <tt>viscosity_value</tt> = 1e-3
+!!
+!! and can all be changed in the parfile. @endnote
 submodule (mod_equilibrium) smod_equil_couette_flow
   implicit none
 

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
+!!
+!! - <tt>k2</tt> = 0.12
+!! - <tt>k3</tt> = 0
+!! - <tt>cte_rho0</tt> = 1
+!! - <tt>cte_T0</tt> = 1
+!! - <tt>cte_B02</tt> = 1
+!! - <tt>cte_B03</tt> = 0 : guide field parameter.
+!! - <tt>alpha</tt> = 1 : used to set the width of the current sheet.
+!! - <tt>eq_bool</tt> = 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
 

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")