dedicated physics object + functions for resistivity

CMakeLists.txt

@@ -12,14 +12,40 @@ if (${Coverage})
     message(STATUS "Building with code coverage enabled.")
     message(STATUS "Default Fortran flags are disabled, optimisation disabled.")
     message(STATUS "====================")
-    set (CMAKE_Fortran_FLAGS "--coverage -o0 -g -cpp -ffree-line-length-none")
+    set(
+        CMAKE_Fortran_FLAGS
+        "--coverage -o0 -g -cpp -ffree-line-length-none -Wno-unused-dummy-argument"
+    )
 elseif(Debug)
-    set(CMAKE_Fortran_FLAGS "-fcheck=all -fbounds-check -Wall -ffree-line-length-none \
-                             -Wextra -Wconversion -pedantic -fbacktrace -cpp")
+    set(
+        CMAKE_Fortran_FLAGS
+        "-fcheck=all \
+        -fbounds-check \
+        -Wall \
+        -ffree-line-length-none \
+        -Wextra \
+        -Wconversion \
+        -pedantic \
+        -fbacktrace \
+        -cpp \
+        -Wno-unused-dummy-argument"
+    )
 else() # Release
-    set(CMAKE_Fortran_FLAGS "-fcheck=all -Wall -Wextra -ffree-line-length-none \
-                             -O3 -funroll-loops -ftree-vectorize -Wconversion \
-                             -pedantic -fbacktrace -cpp")
+    set(
+        CMAKE_Fortran_FLAGS
+        "-fcheck=all \
+        -Wall \
+        -Wextra \
+        -ffree-line-length-none \
+        -O3 \
+        -funroll-loops \
+        -ftree-vectorize \
+        -Wconversion \
+        -pedantic \
+        -fbacktrace \
+        -cpp \
+        -Wno-unused-dummy-argument"
+    )
 endif()
 
 if (DEFINED ENV{LEGOLASDIR})

src/physics/mod_physics.f08

@@ -0,0 +1,28 @@
+module mod_physics
+  use mod_resistivity, only: resistivity_t, new_resistivity
+  implicit none
+
+  private
+
+  type, public :: physics_t
+    type(resistivity_t) :: resistivity
+  contains
+    procedure, public :: delete
+  end type physics_t
+
+  public :: new_physics
+
+contains
+
+  function new_physics() result(physics)
+    type(physics_t) :: physics
+    physics%resistivity = new_resistivity()
+  end function new_physics
+
+
+  pure subroutine delete(this)
+    class(physics_t), intent(inout) :: this
+    call this%resistivity%delete()
+  end subroutine delete
+
+end module mod_physics

src/physics/mod_physics_utils.f08

@@ -0,0 +1,80 @@
+module mod_physics_utils
+  use mod_global_variables, only: dp
+  use mod_physical_constants, only: dpi
+  use mod_settings, only: settings_t
+  use mod_background, only: background_t
+  implicit none
+
+  private
+
+  interface
+    real(dp) function physics_i(x, settings, background)
+      use mod_global_variables, only: dp
+      use mod_settings, only: settings_t
+      use mod_background, only: background_t
+      real(dp), intent(in) :: x
+      type(settings_t), intent(in) :: settings
+      type(background_t), intent(in) :: background
+    end function physics_i
+  end interface
+
+  public :: physics_i
+  public :: get_dropoff
+  public :: get_dropoff_dr
+
+contains
+
+  real(dp) function get_dropoff(x, cte_value, settings, background)
+    real(dp), intent(in) :: x
+    real(dp), intent(in) :: cte_value
+    type(settings_t), intent(in) :: settings
+    type(background_t), intent(in) :: background
+    real(dp) :: edge_dist, width, sleft, sright, shift, stretch
+
+    edge_dist = settings%physics%dropoff_edge_dist
+    width = settings%physics%dropoff_width
+    sleft = settings%grid%get_grid_start() + 0.5_dp * width + edge_dist
+    sright = settings%grid%get_grid_end() - 0.5_dp * width - edge_dist
+    shift = cte_value * tanh(-dpi) / (tanh(-dpi) - tanh(dpi))
+    stretch = cte_value / (tanh(dpi) - tanh(-dpi))
+
+    if (sleft - 0.5_dp * width <= x .and. x <= sleft + 0.5_dp * width) then
+      get_dropoff = shift + stretch * tanh(2.0_dp * dpi * (x - sleft) / width)
+    else if (sleft + 0.5_dp * width < x .and. x < sright - 0.5_dp * width) then
+      get_dropoff = cte_value
+    else if (sright - 0.5_dp * width <= x .and. x <= sright + 0.5_dp * width) then
+      get_dropoff = shift + stretch * tanh(2.0_dp * dpi * (sright - x) / width)
+    else
+      get_dropoff = 0.0_dp
+    end if
+  end function get_dropoff
+
+
+  real(dp) function get_dropoff_dr(x, cte_value, settings, background)
+    real(dp), intent(in) :: x
+    real(dp), intent(in) :: cte_value
+    type(settings_t), intent(in) :: settings
+    type(background_t), intent(in) :: background
+    real(dp) :: edge_dist, width, sleft, sright, shift, stretch
+
+    edge_dist = settings%physics%dropoff_edge_dist
+    width = settings%physics%dropoff_width
+    sleft = settings%grid%get_grid_start() + 0.5_dp * width + edge_dist
+    sright = settings%grid%get_grid_end() - 0.5_dp * width - edge_dist
+    shift = cte_value * tanh(-dpi) / (tanh(-dpi) - tanh(dpi))
+    stretch = cte_value / (tanh(dpi) - tanh(-dpi))
+
+    if (sleft - 0.5_dp * width <= x .and. x <= sleft + 0.5_dp * width) then
+      get_dropoff_dr = (2.0_dp * dpi * stretch / width) / cosh( &
+        2.0_dp * dpi * (x - sleft) / width &
+      )**2
+    else if (sright - 0.5_dp * width <= x .and. x <= sright + 0.5_dp * width) then
+      get_dropoff_dr = (-2.0_dp * dpi * stretch / width) / cosh( &
+        2.0_dp * dpi * (sright - x) / width &
+      )**2
+    else
+      get_dropoff_dr = 0.0_dp
+    end if
+  end function get_dropoff_dr
+
+end module mod_physics_utils

src/physics/mod_resistivity.f08

@@ -3,16 +3,26 @@
 !! and sets the resistivity values based on the equilibrium configuration.
 module mod_resistivity
   use mod_global_variables, only: dp
-  use mod_physical_constants, only: dpi, Z_ion, coulomb_log
+  use mod_physical_constants, only: dpi, Z_ion, coulomb_log, ec_cgs, me_cgs, kB_cgs
   use mod_settings, only: settings_t
   use mod_background, only: background_t
   use mod_function_utils, only: from_function
+  use mod_physics_utils, only: physics_i, get_dropoff, get_dropoff_dr
   use mod_grid, only: grid_gauss
   implicit none
 
   private
 
+  type, public :: resistivity_t
+    procedure(physics_i), pointer, nopass :: eta
+    procedure(physics_i), pointer, nopass :: detadT
+    procedure(physics_i), pointer, nopass :: detadr
+  contains
+    procedure, public :: delete
+  end type resistivity_t
+
   public :: set_resistivity_values
+  public :: new_resistivity
 
 contains
 
@@ -179,4 +189,99 @@ subroutine get_constants(ec, me, kB)
     kB = kB_cgs
   end subroutine get_constants
 
+
+  function new_resistivity() result(resistivity)
+    type(resistivity_t) :: resistivity
+
+    resistivity%eta => default_eta
+    resistivity%detadT => default_detadT
+    resistivity%detadr => default_detadr
+  end function new_resistivity
+
+  !> Default profile for the resistivity $\eta$. Returns either the Spitzer resistivity
+  !! based on the equilibrium temperature profile, or a fixed resistivity value if
+  !! that was specified through the settings. Values are normalised on return.
+  real(dp) function default_eta(x, settings, background)
+    real(dp), intent(in) :: x
+    type(settings_t), intent(in) :: settings
+    type(background_t), intent(in) :: background
+    real(dp) :: unit_temperature, unit_resistivity
+    real(dp) :: T0
+
+    if (settings%physics%resistivity%has_fixed_resistivity()) then
+      default_eta = settings%physics%resistivity%get_fixed_resistivity()
+      if (settings%physics%resistivity%use_dropoff) then
+        default_eta = get_dropoff(x, default_eta, settings, background)
+      end if
+      return
+    end if
+
+    unit_temperature = settings%units%get_unit_temperature()
+    unit_resistivity = settings%units%get_unit_resistivity()
+    T0 = background%temperature%T0(x) * unit_temperature
+    default_eta = ( &
+      (4.0_dp / 3.0_dp) &
+      * sqrt(2.0_dp * dpi) &
+      * Z_ion &
+      * ec_cgs**2 &
+      * sqrt(me_cgs) &
+      * coulomb_log &
+      / (kB_cgs * T0)**(3.0_dp / 2.0_dp) &
+    ) / unit_resistivity
+  end function default_eta
+
+
+  !> Default profile for the derivative of the resistivity $\eta$ with respect to
+  !! the temperature $T$. Returns normalised values, or zero if a fixed resistivity
+  !! was specified in the settings.
+  real(dp) function default_detadT(x, settings, background)
+    real(dp), intent(in) :: x
+    type(settings_t), intent(in) :: settings
+    type(background_t), intent(in) :: background
+    real(dp) :: unit_temperature, unit_detadT
+    real(dp) :: T0
+
+    if (settings%physics%resistivity%has_fixed_resistivity()) then
+      default_detadT = 0.0_dp
+      return
+    end if
+    unit_temperature = settings%units%get_unit_temperature()
+    unit_detadT = settings%units%get_unit_resistivity() / unit_temperature
+    T0 = background%temperature%T0(x) * unit_temperature
+    default_detadT = ( &
+      -2.0_dp * sqrt(2.0_dp * dpi) * Z_ion * ec_cgs**2 * sqrt(me_cgs) * coulomb_log &
+      / (kB_cgs**(3.0_dp / 2.0_dp) * T0**(5.0_dp / 2.0_dp)) &
+    ) / unit_detadT
+  end function default_detadT
+
+
+  !> Default profile for the derivative of the resistivity $\eta$ with respect to
+  !! position. Returns zero unless a dropoff profile was chosen.
+  real(dp) function default_detadr(x, settings, background)
+    real(dp), intent(in) :: x
+    type(settings_t), intent(in) :: settings
+    type(background_t), intent(in) :: background
+
+    default_detadr = 0.0_dp
+    if (settings%physics%resistivity%has_fixed_resistivity()) then
+      if (settings%physics%resistivity%use_dropoff) then
+        default_detadr = get_dropoff_dr( &
+          x, &
+          settings%physics%resistivity%get_fixed_resistivity(), &
+          settings, &
+          background &
+        )
+      end if
+      return
+    end if
+  end function default_detadr
+
+
+  pure subroutine delete(this)
+    class(resistivity_t), intent(inout) :: this
+    nullify(this%eta)
+    nullify(this%detadT)
+    nullify(this%detadr)
+  end subroutine delete
+
 end module mod_resistivity