Hall by momentum equation substitution

.github/workflows/docs.yml

@@ -23,7 +23,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install ford graphviz
+        pip install ford==6.1.6 graphviz
         sudo apt-get install graphviz
         pip install -U sphinx
         pip install numpydoc sphinx-rtd-theme lazy-object-proxy==1.4 sphinx-autoapi

docs/general/parameter_file.md

@@ -76,7 +76,6 @@ This namelist includes all physics-related variables.
 | viscosity_value       | real    | constant value for the viscosity | 0 |
 | incompressible        | logical | if `.true.`, uses an incompressible approximation | `.false.` |
 | hall_mhd              | logical | inclusion of Hall effects | `.false.` |
-| elec_pressure         | logical | whether to include the electron pressure term in Ohm's law if `hall_mhd = .true.` | `.false.` |
 | elec_intertia         | logical | whether to include the electron intertia term in Ohm's law if `hall_mhd = .true.` | `.false.` |
 
 ## solvelist

post_processing/pylbo/automation/defaults.py

@@ -57,9 +57,10 @@
         ("viscosity_value", (int, np.integer, float)),
         ("hall_mhd", bool),
         ("hall_dropoff", bool),
-        ("elec_pressure", bool),
         ("elec_inertia", bool),
         ("inertia_dropoff", bool),
+        ("electron_fraction", (int, np.integer, float)),
+        ("incompressible", bool),
     ],
     "unitslist": [
         ("cgs_units", bool),

post_processing/pylbo/data_containers.py

@@ -400,7 +400,7 @@ def get_tube_speed(self, which_values=None):
         else:
             cA = self.get_alfven_speed()
             cs = self.get_sound_speed()
-            ct = cs * cA / np.sqrt(cs ** 2 + cA ** 2)
+            ct = cs * cA / np.sqrt(cs**2 + cA**2)
             return self._get_values(ct, which_values)
 
     def get_reynolds_nb(self, which_values=None):

post_processing/pylbo/utilities/toolbox.py

@@ -170,20 +170,20 @@ def solve_cubic_exact(a, b, c, d):
     q = c / a
     r = d / a
     Aterm = (
-        -2 * p ** 3
+        -2 * p**3
         + 9 * p * q
         - 27 * r
         + 3
         * np.sqrt(3)
         * np.sqrt(
-            -(p ** 2) * q ** 2
-            + 4 * q ** 3
-            + 4 * p ** 3 * r
+            -(p**2) * q**2
+            + 4 * q**3
+            + 4 * p**3 * r
             - 18 * p * q * r
-            + 27 * r ** 2
+            + 27 * r**2
         )
     ) ** (1 / 3) / (3 * 2 ** (1 / 3))
-    Bterm = (-(p ** 2) + 3 * q) / (9 * Aterm)
+    Bterm = (-(p**2) + 3 * q) / (9 * Aterm)
     cterm_min = (-1 - np.sqrt(3) * 1j) / 2
     cterm_pos = (-1 + np.sqrt(3) * 1j) / 2
     x1 = -p / 3 + Aterm - Bterm

post_processing/pylbo/visualisation/continua.py

@@ -35,7 +35,7 @@ def calculate_continua(ds):
 
     # Alfven and slow continuum (squared)
     alfven2 = (1 / rho) * ((k2 * B02 / ds.scale_factor) + k3 * B03) ** 2
-    slow2 = (gamma * p / (gamma * p + B0 ** 2)) * alfven2
+    slow2 = (gamma * p / (gamma * p + B0**2)) * alfven2
     # doppler shift equals dot product of k and v
     doppler = k2 * v02 / ds.scale_factor + k3 * v03
     slow_min = -np.sqrt(slow2)
@@ -50,7 +50,7 @@ def calculate_continua(ds):
     else:
         kpara = 0
     cs2 = gamma * p / rho  # sound speed
-    vA2 = B0 ** 2 / rho  # Alfvén speed
+    vA2 = B0**2 / rho  # Alfvén speed
     ci2 = p / rho  # isothermal sound speed
     dLdT = ds.equilibria["dLdT"]
     dLdrho = ds.equilibria["dLdrho"]
@@ -71,20 +71,20 @@ def calculate_continua(ds):
     elif slow_is_zero:
         thermal = 1j * (gamma - 1) * (rho * dLdrho - dLdT * (ci2 + vA2)) / (cs2 + vA2)
     else:
-        sigma_A2 = kpara ** 2 * vA2  # Alfvén frequency
+        sigma_A2 = kpara**2 * vA2  # Alfvén frequency
         sigma_c2 = cs2 * sigma_A2 / (vA2 + cs2)  # cusp frequency
         sigma_i2 = ci2 * sigma_A2 / (vA2 + ci2)  # isothermal cusp frequency
 
         # thermal and slow continuum are coupled in a third degree polynomial in omega,
         # coeffi means the coefficient corresponding to the term omega^i
         coeff3 = rho * (cs2 + vA2) * 1j / (gamma - 1)
         coeff2 = (
-            -(kappa_para * kpara ** 2 + rho * dLdT) * (ci2 + vA2) + rho ** 2 * dLdrho
+            -(kappa_para * kpara**2 + rho * dLdT) * (ci2 + vA2) + rho**2 * dLdrho
         )
         coeff1 = -rho * (cs2 + vA2) * sigma_c2 * 1j / (gamma - 1)
-        coeff0 = (kappa_para * kpara ** 2 + rho * dLdT) * (
+        coeff0 = (kappa_para * kpara**2 + rho * dLdT) * (
             ci2 + vA2
-        ) * sigma_i2 - rho ** 2 * dLdrho * sigma_A2
+        ) * sigma_i2 - rho**2 * dLdrho * sigma_A2
         # we have to solve this equation "gauss_gridpts" times.
         # the thermal continuum corresponds to the (only) purely imaginary solution,
         # modified slow continuum are other two solutions

post_processing/pylbo/visualisation/derived_eigenfunctions.py

@@ -78,7 +78,7 @@ def _get_title(self, ef_name):
         ef_name = ef_name.replace("curl", "\\nabla\\times")
         ef_name = ef_name.replace("para", "\\parallel")
         ef_name = ef_name.replace("perp", "\\perp")
-        name = fr"{part}(${ef_name}$)"
+        name = rf"{part}(${ef_name}$)"
         return name
 
     def _mark_points_without_data_written(self):

post_processing/pylbo/visualisation/profiles.py

@@ -177,7 +177,7 @@ def draw(self):
             + b02 * db02
             + b03 * db03
             + rho * g
-            - (dr_scale / r_scale) * (rho * v02 ** 2 - b02 ** 2)
+            - (dr_scale / r_scale) * (rho * v02**2 - b02**2)
         )
         equil_force[np.where(abs(equil_force) <= 1e-16)] = 0
         self.ax.plot(self.data.grid_gauss, equil_force, **self.kwargs)

post_processing/pylbo/visualisation/spectra.py

@@ -329,7 +329,7 @@ def _get_ydata(self):
         """
         ydata = np.empty(len(self.dataseries), dtype=object)
         for i, ds in enumerate(self.dataseries):
-            ydata[i] = ds.eigenvalues ** self._w_pow
+            ydata[i] = ds.eigenvalues**self._w_pow
             if self.use_real_parts:
                 ydata[i] = ydata[i].real
             else:

src/boundaries/smod_essential_boundaries.f08

@@ -1,4 +1,7 @@
 submodule (mod_boundary_manager) smod_essential_boundaries
+  use mod_global_variables, only: dp_LIMIT, boundary_type, geometry, coaxial
+  use mod_equilibrium_params, only: k2, k3
+
   implicit none
 
 contains
@@ -11,11 +14,23 @@
     ! on the left side we have a zero in a quadratic basis function (number 2), which
     ! zeroes out the odd rows/columns. We explicitly handle this by introducing an
     ! element on the diagonal for these indices as well. The odd numbered indices for
-    ! the quadratic elements correspond to rho, v2, v2, T, a1 = 1, 5, 6, 9, 11.
+    ! the quadratic elements correspond to rho, v2, v3, T, a1 = 1, 5, 7, 9, 11.
     call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[1, 5, 7, 9, 11])
-    ! wall or regularity conditions: v1, a2 and a3 have to be zero. These are cubic
+    ! wall or regularity conditions: v1 and (k3 * a2 - k2 * a3) have to be zero. These are cubic
     ! elements, so we force non-zero basis functions (odd rows & columns) to zero.
-    call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[3, 13, 15])
+    ! for wall we force a2 = a3 = 0 regardless of k2 and k3, for wall_weak we
+    ! only force a2 resp. a3 to zero if k3 resp. k2 is nonzero
+    if (boundary_type == 'wall') then
+      call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[3, 13, 15])
+    else if (boundary_type == 'wall_weak') then
+      call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[3])
+      if (abs(k2) > dp_LIMIT) then
+        call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[15])
+      end if
+      if (abs(k3) > dp_LIMIT) then
+        call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[13])
+      end if
+    end if
     ! if T boundary conditions are needed, set even row/colum (quadratic)
     if (apply_T_bounds) then
       call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[10])
@@ -33,8 +48,18 @@
 
     diagonal_factor = get_diagonal_factor(matrix)
 
-    ! for a fixed wall v1, a2 and a3 should be zero, same reasoning as for left side.
-    call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[19, 29, 31])
+    ! for a fixed wall v1 and (k3 * a2 - k2 * a3) should be zero, same reasoning as for left side.
+    if (boundary_type == 'wall') then
+      call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[19, 29, 31])
+    else if (boundary_type == 'wall_weak') then
+      call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[19])
+      if (abs(k2) > dp_LIMIT) then
+        call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[31])
+      end if
+      if (abs(k3) > dp_LIMIT) then
+        call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[29])
+      end if
+    end if
     ! T condition
     if (apply_T_bounds) then
       call set_row_col_to_value(quadblock, val=diagonal_factor, idxs=[26])
@@ -84,4 +109,4 @@ function get_diagonal_factor(matrix) result(diagonal_factor)
   end function get_diagonal_factor
 
 
-end submodule smod_essential_boundaries
+end submodule smod_essential_boundaries

src/boundaries/smod_natural_bounds_hall.f08

@@ -3,77 +3,6 @@
 
 contains
 
-  module procedure add_natural_hall_terms
-    use mod_global_variables, only: hall_mhd, elec_pressure
-    use mod_equilibrium, only: hall_field
-
-    real(dp)  :: eps
-    real(dp)  :: rho, T0, B01, B02, B03
-    real(dp)  :: eta_H
-
-    if (.not. hall_mhd) then
-      return
-    end if
-
-    eps = eps_grid(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)
-
-    ! ==================== Quadratic * Quadratic ====================
-    call reset_factor_positions(new_size=1)
-    ! H(6, 6)
-    factors(1) = -eta_H * eps * (k3 * B02 - k2 * B03 / eps) / rho
-    positions(1, :) = [6, 6]
-    call subblock(quadblock, factors, positions, weight, h_quad, h_quad)
-    if (elec_pressure) then
-      call reset_factor_positions(new_size=2)
-      ! H(6, 1)
-      factors(1) = -eta_H * T0 / rho
-      positions(1, :) = [6, 1]
-      ! H(6, 5)
-      factors(2) = -eta_H
-      positions(2, :) = [6, 5]
-      call subblock(quadblock, factors, positions, weight, h_quad, h_quad)
-    end if
-
-    ! ==================== 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 procedure add_natural_hall_terms
-
   module procedure add_natural_hall_Bterms
     use mod_global_variables, only: hall_mhd, elec_inertia
     use mod_equilibrium, only: hall_field
@@ -112,4 +41,55 @@
 
   end procedure add_natural_hall_Bterms
 
+  module procedure add_natural_hall_terms
+    use mod_global_variables, only: hall_mhd, viscosity, viscosity_value
+    use mod_equilibrium, only: hall_field
+
+    real(dp)  :: eps, deps
+    real(dp)  :: rho
+    real(dp)  :: eta_H, mu
+
+    if (.not. (hall_mhd .and. viscosity)) then
+      return
+    end if
+
+    eps = eps_grid(grid_idx)
+    deps = d_eps_grid_dr(grid_idx)
+    rho = rho_field % rho0(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)
+
+  end procedure add_natural_hall_terms
+
 end submodule

src/boundaries/smod_natural_bounds_viscosity.f08

@@ -4,7 +4,8 @@
 contains
 
   module procedure add_natural_viscosity_terms
-    use mod_global_variables, only: viscosity, viscous_heating, viscosity_value
+    use mod_global_variables, only: viscosity, viscous_heating, viscosity_value, &
+                                    gamma_1, incompressible
     use mod_equilibrium, only: v_field
 
     real(dp)  :: eps, deps
@@ -62,8 +63,8 @@
     positions(1, :) = [4, 4]
     ! Sigma(5, 4)
     factors(2) = (0.0d0, 0.0d0)
-    if (viscous_heating) then
-      factors(2) = 2.0d0 * ic * mu * dv03
+    if (viscous_heating .and. (.not. incompressible)) then
+      factors(2) = 2.0d0 * ic * gamma_1 * mu * dv03
     end if
     positions(2, :) = [5, 4]
     call subblock(quadblock, factors, positions, weight, h_quad, h_quad)
@@ -72,8 +73,8 @@
     call reset_factor_positions(new_size=1)
     ! Sigma(5, 2)
     factors(1) = (0.0d0, 0.0d0)
-    if (viscous_heating) then
-      factors(1) = 2.0d0 * mu * dv01
+    if (viscous_heating .and. (.not. incompressible)) then
+      factors(1) = 2.0d0 * gamma_1 * mu * dv01
     end if
     positions(1, :) = [5, 2]
     call subblock(quadblock, factors, positions, weight, h_quad, h_cubic)

src/dataIO/datfile_changes.txt

@@ -38,3 +38,6 @@ added this in a backwards compatible way to the pylbo reader
 1) added eigenfunction subset info to datfile header
 2) added ef_written_flags and ef_written_idxs to datfile, before writing eigenfunctions
 Changes needed an update of the pylbo reader
+
+=== version 1.2.0 --> 1.2.1 ===
+added electron_fraction to param_list after g

src/dataIO/mod_input.f08

@@ -48,7 +48,7 @@ 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_pressure, elec_inertia, inertia_dropoff
+        hall_mhd, 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,8 +223,8 @@ subroutine print_console_info()
 
     call log_message("hall mhd           : " // str(hall_mhd))
     if (hall_mhd) then
-      call log_message("    electron pressure : " // str(elec_pressure))
-      call log_message("    electron intertia : " // str(elec_inertia))
+      call log_message("    electron fraction : " // str(electron_fraction))
+      call log_message("    electron inertia  : " // str(elec_inertia))
     end if
 
     call log_message("            << Solver settings >>")