From 1ea672a1b0eb4f9001499b0c08277f0a6fbcae84 Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Fri, 24 Mar 2023 12:45:59 +0100 Subject: [PATCH] dedicated physics object + functions for resistivity --- CMakeLists.txt | 38 +++++++++-- src/physics/mod_physics.f08 | 28 ++++++++ src/physics/mod_physics_utils.f08 | 80 ++++++++++++++++++++++ src/physics/mod_resistivity.f08 | 107 +++++++++++++++++++++++++++++- 4 files changed, 246 insertions(+), 7 deletions(-) create mode 100644 src/physics/mod_physics.f08 create mode 100644 src/physics/mod_physics_utils.f08 diff --git a/CMakeLists.txt b/CMakeLists.txt index 2f0800be..d8da103c 100644 --- a/CMakeLists.txt +++ b/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}) diff --git a/src/physics/mod_physics.f08 b/src/physics/mod_physics.f08 new file mode 100644 index 00000000..0999facf --- /dev/null +++ b/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 diff --git a/src/physics/mod_physics_utils.f08 b/src/physics/mod_physics_utils.f08 new file mode 100644 index 00000000..02618b27 --- /dev/null +++ b/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 diff --git a/src/physics/mod_resistivity.f08 b/src/physics/mod_resistivity.f08 index adab4823..f04bb7d4 100644 --- a/src/physics/mod_resistivity.f08 +++ b/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