From 356ad04dcebada6f62ea42bfa22e5517bc02e8f6 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Wed, 19 Jul 2023 16:41:16 -0700 Subject: [PATCH] added simple ODE stepper for rc model --- src/CMakeLists.txt | 1 + src/clima_ode.f90 | 220 ++++++++++++++++++++++++++++++++++++++++++++ src/rc/clima_rc.f90 | 81 +++++++++++----- 3 files changed, 281 insertions(+), 21 deletions(-) create mode 100644 src/clima_ode.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 11710a9..ca7d5a0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,7 @@ set(BASE_SOURCES ${CMAKE_CURRENT_BINARY_DIR}/clima_version.f90 clima_const.f90 clima_useful.f90 + clima_ode.f90 clima_eqns.f90 clima_saturationdata.f90 clima_eqns_water.f90 diff --git a/src/clima_ode.f90 b/src/clima_ode.f90 new file mode 100644 index 0000000..50c2df4 --- /dev/null +++ b/src/clima_ode.f90 @@ -0,0 +1,220 @@ +module clima_ode + use clima_const, only: dp + implicit none + private + public :: ODEStepper + + !> A simple ODE stepper class using Heun's method (https://en.wikipedia.org/wiki/Heun%27s_method). + !> Adaptive stepsize selection follows pages 167 - 169 in + !> "Solving Ordinary Differential Equations I" by Hairer, Norsett and Wanner + type :: ODEStepper + + procedure(rhs_fcn), pointer :: fcn !! RHS function + integer :: n !! number of equations + + ! Current state of ODE integrator + real(dp) :: t !! current time + real(dp), allocatable :: u(:) !! current value of variables + real(dp) :: h !! current stepsize + + ! settings + real(dp) :: rtol = 1.0e-3_dp !! relative tolerance + real(dp) :: atol = 1.0e-6_dp !! absolute tolerance + real(dp) :: facmax = 3.0_dp + real(dp) :: facmin = 0.3_dp + real(dp) :: safety_factor = 0.8_dp + integer :: max_error_test_failures = 10 + + ! statistics for solver + integer :: nfevals !! number of function evalutions + integer :: naccept !! number of accepted steps + integer :: nreject !! number of steps that were rejected + + ! work space + real(dp), allocatable :: du(:) + real(dp), allocatable :: du1(:) + real(dp), allocatable :: u1(:) + real(dp), allocatable :: u_new(:) + real(dp) :: t_new + real(dp), allocatable :: sc(:) + real(dp) :: h_new + + contains + procedure :: initialize => ODEStepper_initialize + procedure :: initial_stepsize => ODEStepper_initial_stepsize + procedure :: step => ODEStepper_step + end type + + abstract interface + subroutine rhs_fcn(self, t, u, du, ierr) + import :: dp, ODEStepper + class(ODEStepper), intent(inout) :: self + real(dp), intent(in) :: t + real(dp), intent(in) :: u(:) + real(dp), intent(out) :: du(:) + integer, intent(out) :: ierr + end subroutine + end interface + +contains + + subroutine ODEStepper_initialize(self, fcn, t0, u0, ierr) + class(ODEStepper), intent(inout) :: self + procedure(rhs_fcn) :: fcn + real(dp), intent(in) :: t0 + real(dp), intent(in) :: u0(:) + integer, intent(out) :: ierr + + self%fcn => fcn + self%n = size(u0) + + self%t = t0 + if (allocated(self%u)) deallocate(self%u) + self%u = u0 + + self%nfevals = 0 + self%naccept = 0 + self%nreject = 0 + + if (allocated(self%du)) then + deallocate(self%du) + deallocate(self%du1) + deallocate(self%u1) + deallocate(self%u_new) + deallocate(self%sc) + endif + allocate(self%du(self%n)) + allocate(self%du1(self%n)) + allocate(self%u1(self%n)) + allocate(self%u_new(self%n)) + allocate(self%sc(self%n)) + + self%h = self%initial_stepsize(t0, u0, ierr) + if (ierr < 0) return + + end subroutine + + subroutine ODEStepper_step(self, ierr) + class(ODEStepper), intent(inout) :: self + integer, intent(out) :: ierr + + logical :: accept_step + real(dp) :: facmax, err + integer :: ierr1 + integer :: i, j + + facmax = self%facmax + + ierr = 0 + do j = 1,self%max_error_test_failures + + ! Try to perform a step of Heun's Method + self%t_new = self%t + self%h + call self%fcn(self%t, self%u, self%du, ierr1) + if (ierr1 < 0) then + ierr = -2 + return + endif + self%nfevals = self%nfevals + 1 + self%u1 = self%u + self%h*self%du ! first order approximation to u_new + call self%fcn(self%t_new, self%u1, self%du1, ierr1) + if (ierr1 < 0) then + ierr = -2 + return + endif + self%nfevals = self%nfevals + 1 + self%u_new = self%u + 0.5_dp*self%h*(self%du + self%du1) ! second-order approximation to u_new + + ! Error and timestep control. Methods generally follow Page 167 - 168 in + ! "Solving Ordinary Differential Equations I" by Hairer, Norsett and Wanner + do i = 1,self%n + self%sc(i) = self%atol + self%rtol*max(abs(self%u(i)), abs(self%u_new(i))) + enddo + + if (all(abs(self%u_new - self%u1) <= self%sc)) then + accept_step = .true. + facmax = self%facmax + else + accept_step = .false. + facmax = 1.0_dp + endif + + ! Compute new stepsize + err = sqrt((1.0_dp/real(self%n,dp))*sum(((self%u_new - self%u1)/self%sc)**2.0_dp)) + self%h_new = self%safety_factor*self%h*min(facmax, max(self%facmin, (1.0_dp/err)**(0.5_dp))) + self%h = self%h_new + + if (accept_step) then + self%t = self%t_new + self%u = self%u_new + self%naccept = self%naccept + 1 + return + else + self%nreject = self%nreject + 1 + endif + + enddo + + ! Will only get here if step failed for `self%max_error_test_failures` + ierr = -1 + return + + end subroutine + + !> Computes the initial stepsize. This follows 169 in + !> "Solving Ordinary Differential Equations I" by Hairer, Norsett and Wanner + function ODEStepper_initial_stepsize(self, t0, u0, ierr) result(h_init) + class(ODEStepper), intent(inout) :: self + real(dp), intent(in) :: t0 + real(dp), intent(in) :: u0(:) + integer, intent(out) :: ierr + real(dp) :: h_init + + integer :: i, ierr1 + real(dp) :: d0, d1, d2 + real(dp) :: h0, h1 + + ierr = 0 + + !~~ Step a) ~~! + call self%fcn(t0, u0, self%du, ierr1) + if (ierr1 < 0) then + ierr = -2 + return + endif + + do i = 1,self%n + self%sc(i) = self%atol + self%rtol*abs(u0(i)) + enddo + d0 = sqrt((1.0_dp/real(self%n,dp))*sum((u0/self%sc)**2.0_dp)) + d1 = sqrt((1.0_dp/real(self%n,dp))*sum((self%du/self%sc)**2.0_dp)) + + !~~ Step b) ~~! + h0 = 0.01_dp*(d0/d1) + if (d0 < 1.0e-5_dp .or. d1 < 1.0e-5_dp) then + h0 = 1.0e-6_dp + endif + + !~~ Step c) ~~! + self%u1 = u0 + h0*self%du + call self%fcn(t0+h0, self%u1, self%du1, ierr1) + if (ierr1 < 0) then + ierr = -2 + return + endif + + !~~ Step d) ~~! + d2 = sqrt((1.0_dp/real(self%n,dp))*sum(((self%du1 - self%du)/self%sc)**2.0_dp))/h0 + + !~~ Step e) ~~! + h1 = (0.01_dp/max(d1,d2))**(0.5_dp) + if (max(d1,d2) < 1.0e-15_dp) then + h1 = max(1.0e-6_dp, h0*1.0e-3_dp) + endif + + !~~ Step f) ~~! + h_init = min(100.0_dp*h0, h1) + + end function + +end module \ No newline at end of file diff --git a/src/rc/clima_rc.f90 b/src/rc/clima_rc.f90 index 0392a72..bd36b5f 100644 --- a/src/rc/clima_rc.f90 +++ b/src/rc/clima_rc.f90 @@ -2,6 +2,7 @@ module clima_rc use clima_types, only: Species use clima_radtran, only: Radtran use clima_const, only: dp, s_str_len + use clima_ode, only: ODEStepper implicit none private @@ -19,6 +20,10 @@ module clima_rc module procedure :: create_InitialConditions end interface + type, extends(ODEStepper) :: ODEStepper_custom + type(RadiativeConvectiveClimate), pointer :: c => null() + end type + type :: RadiativeConvectiveClimate ! Species @@ -72,6 +77,9 @@ module clima_rc real(dp), allocatable :: pdensities_r(:,:) !! (nz_r,np) real(dp), allocatable :: radii_r(:,:) !! (nz_r,np) + ! ODE stepper + type(ODEStepper_custom), allocatable :: ode + contains procedure :: make_profile => RadiativeConvectiveClimate_make_profile procedure :: TOA_fluxes => RadiativeConvectiveClimate_TOA_fluxes @@ -465,7 +473,7 @@ function create_InitialConditions(condensible_names, condensible_P, f_i, pdensit subroutine RadiativeConvectiveClimate_initialize_stepper(self, condensible_names, condensible_P, condensible_RH, & f_i, pdensities, radii, T_surf_guess, err) use clima_rc_adiabat, only: make_profile_rc - class(RadiativeConvectiveClimate), intent(inout) :: self + class(RadiativeConvectiveClimate), target, intent(inout) :: self character(*), intent(in) :: condensible_names(:) real(dp), intent(in) :: condensible_P(:) real(dp), intent(in) :: condensible_RH(:) @@ -476,37 +484,68 @@ subroutine RadiativeConvectiveClimate_initialize_stepper(self, condensible_names character(:), allocatable, intent(out) :: err real(dp) :: T_surf + real(dp), allocatable :: u(:) + integer :: ierr - ! Initial guess + ! Initial guess. This sets everything in the atmosphere, including `self%trop_ind` T_surf = self%surface_temperature(condensible_names, condensible_P, condensible_RH, & f_i, pdensities, radii, T_surf_guess, err) if (allocated(err)) return - block - real(dp), allocatable :: dT_dt(:) - real(dp), allocatable :: T(:) - allocate(dT_dt(size(self%T)+1)) - allocate(T(size(self%T)+1)) - T(1) = T_surf - T(2:) = self%T - call RadiativeConvectiveClimate_rhs(self, T, dT_dt, err) - if (allocated(err)) return - endblock - ! if (allocated(self%init)) deallocate(self%init) - ! allocate(self%init) - ! self%init = InitialConditions(condensible_names, condensible_P, f_i, pdensities, T_init, radii) + ! Initialize ODE stepper. + allocate(u(self%nz+1)) + u(1) = T_surf + u(2:) = self%T + if (allocated(self%ode)) deallocate(self%ode) + allocate(self%ode) + self%ode%c => self ! associate pointer to self + call self%ode%initialize(ODEStepper_rhs, 0.0_dp, u, ierr) ! initialize at t = 0 + if (ierr < 0) then + err = "ODEStepper failed to initialize" + return + endif + + end subroutine + + subroutine ODEStepper_rhs(self, t, u, du, ierr) + class(ODEStepper), intent(inout) :: self + real(dp), intent(in) :: t + real(dp), intent(in) :: u(:) + real(dp), intent(out) :: du(:) + integer, intent(out) :: ierr + + character(:), allocatable :: err + + ierr = 0 + select type (self) + class is (ODEStepper_custom) + call RadiativeConvectiveClimate_rhs(self%c, u, du, err) + if (allocated(err)) then + ierr = -1 + return + endif + end select end subroutine - ! function RadiativeConvectiveClimate_step(self) result(tn) - ! class(RadiativeConvectiveClimate), intent(inout) :: self - ! real(dp) :: tn + function RadiativeConvectiveClimate_step(self, err) result(tn) + class(RadiativeConvectiveClimate), intent(inout) :: self + character(:), allocatable, intent(out) :: err + real(dp) :: tn - ! ! Do one step from Radiative dT/dt. + integer :: ierr - ! ! Adjust tropopause lapse rate. + ! Do one step from Radiative dT/dt. + call self%ode%step(ierr) + if (ierr < 0) then + err = 'ODE integration step failed' + return + endif + + ! Convective adjustment - ! end function + + end function subroutine check_inputs(nz, ng, np, f_i, pdensities, radii, err) integer, intent(in) :: nz, ng, np