diff --git a/src/atmosphere/photochem_atmosphere.f90 b/src/atmosphere/photochem_atmosphere.f90 index c9ba882..3a8db43 100644 --- a/src/atmosphere/photochem_atmosphere.f90 +++ b/src/atmosphere/photochem_atmosphere.f90 @@ -18,6 +18,7 @@ module photochem_atmosphere contains ! photochem_atmosphere_rhs.f90 + procedure, private :: prep_atm_background_gas procedure :: prep_atmosphere => prep_all_background_gas procedure :: right_hand_side_chem procedure :: production_and_loss @@ -26,6 +27,7 @@ module photochem_atmosphere ! photochem_atmosphere_integrate.f90 procedure :: evolve + procedure :: check_for_convergence procedure :: photochemical_equilibrium procedure :: initialize_stepper procedure :: step @@ -67,6 +69,13 @@ module function create_Atmosphere(mechanism_file, settings_file, flux_file, atmo !~~ photochem_atmosphere_rhs.f90 ~~! + module subroutine prep_atm_background_gas(self, usol_in, usol, molecules_per_particle) + class(Atmosphere), target, intent(inout) :: self + real(dp), intent(in) :: usol_in(:,:) + real(dp), intent(out) :: usol(:,:) + real(dp), intent(out) :: molecules_per_particle(:,:) + end subroutine + !> Given usol_in, the mixing ratios of each species in the atmosphere, !> this subroutine calculates reaction rates, photolysis rates, etc. !> and puts this information into self%wrk. self%wrk contains all the @@ -137,6 +146,13 @@ module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, er logical :: success !! If True, then integration was successful. character(:), allocatable, intent(out) :: err end function + + !> Determines if integration has converged to photochemical equilibrium or not. + module function check_for_convergence(self, err) result(converged) + class(Atmosphere), target, intent(inout) :: self + character(:), allocatable, intent(out) :: err + logical :: converged + end function !> Integrates to photochemical equilibrium starting from self%var%usol_init module subroutine photochemical_equilibrium(self, success, err) diff --git a/src/atmosphere/photochem_atmosphere_integrate.f90 b/src/atmosphere/photochem_atmosphere_integrate.f90 index 50b0cf6..8e465a3 100644 --- a/src/atmosphere/photochem_atmosphere_integrate.f90 +++ b/src/atmosphere/photochem_atmosphere_integrate.f90 @@ -363,7 +363,7 @@ module subroutine photochemical_equilibrium(self, success, err) type(PhotochemWrk), pointer :: wrk real(dp) :: tn - integer :: nsteps + logical :: converged dat => self%dat var => self%var @@ -373,17 +373,23 @@ module subroutine photochemical_equilibrium(self, success, err) if (allocated(err)) return success = .true. - tn = 0 - nsteps = 0 - do while(tn < var%equilibrium_time) + do tn = self%step(err) if (allocated(err)) then ! If the step fails then exit success = .false. exit endif - nsteps = nsteps + 1 - if (nsteps > var%mxsteps) then + + ! Check for convergence + converged = self%check_for_convergence(err) + if (allocated(err)) then + success = .false. + exit + endif + if (converged) exit + + if (wrk%nsteps > var%mxsteps) then ! If we did more steps then max steps, then exit success = .false. exit @@ -398,6 +404,71 @@ module subroutine photochemical_equilibrium(self, success, err) end subroutine + module function check_for_convergence(self, err) result(converged) + use, intrinsic :: iso_c_binding + class(Atmosphere), target, intent(inout) :: self + character(:), allocatable, intent(out) :: err + logical :: converged + + integer :: i,j,ind + type(PhotochemData), pointer :: dat + type(PhotochemVars), pointer :: var + type(PhotochemWrk), pointer :: wrk + + dat => self%dat + var => self%var + wrk => self%wrk + + converged = .false. + if (.not.c_associated(wrk%sun%cvode_mem)) then + err = "You must first initialize the stepper with 'initialize_stepper'" + return + endif + + ! If we reach equilibrium time, then converged + if (wrk%tn > var%equilibrium_time) then + converged = .true. + return + endif + + ! Now consider step history. + + ! Can't do analysis on step 0 + if (wrk%nsteps == 0) return + + ! Find index in history closest to time of interest. Note that this will only + ! consider a limited step history. We can not save all history. + ind = minloc(abs(wrk%t_history - var%conv_hist_factor*wrk%t_history(1)),1) + + ! Can't be current time, so we will check the previous step if needed. + if (ind == 1) ind = 2 + + ! Compute difference between current mixing ratios, and mixing ratios + ! at our index of interest. + do j = 1,var%nz + do i = 1,dat%nq + if (wrk%mix_history(i,j,1) > var%conv_min_mix) then + wrk%dmix(i,j) = abs(wrk%mix_history(i,j,1) - wrk%mix_history(i,j,ind)) + else + ! Ignore small mixing ratios + wrk%dmix(i,j) = 0.0_dp + endif + enddo + enddo + + ! Maximum normalized change + wrk%longdy = maxval(abs(wrk%dmix/wrk%mix_history(:,:,1))) + ! Also consider that change over time + wrk%longdydt = wrk%longdy/(wrk%t_history(1) - wrk%t_history(ind)) + + ! Check for convergence + if (wrk%longdy < var%conv_longdy .and. wrk%longdydt < var%conv_longdydt) then + converged = .true. + return + endif + + end function + module subroutine initialize_stepper(self, usol_start, err) use, intrinsic :: iso_c_binding use fcvode_mod, only: CV_BDF, CV_NORMAL, CV_ONE_STEP, FCVodeInit, FCVodeSStolerances, & @@ -465,6 +536,17 @@ module subroutine initialize_stepper(self, usol_start, err) wrk%sun%yvec(i) = var%lower_fix_mr(i) endif enddo + + ! Load initial conditions into history vector + wrk%nsteps = 0 + wrk%t_history = -1.0_dp + wrk%t_history(1) = 0.0_dp + wrk%mix_history = -1.0_dp + block + real(c_double), pointer :: usol_tmp(:,:) + usol_tmp(1:dat%nq,1:var%nz) => wrk%sun%yvec(1:var%neqs) + wrk%mix_history(:,:,1) = usol_tmp + endblock ! create SUNDIALS N_Vector wrk%sun%sunvec_y => FN_VMake_Serial(neqs_long, wrk%sun%yvec) @@ -549,27 +631,57 @@ module subroutine initialize_stepper(self, usol_start, err) end subroutine module function step(self, err) result(tn) - use iso_c_binding, only: c_null_ptr, c_int, c_double, c_associated - use fcvode_mod, only: CV_ONE_STEP, FCVode + use iso_c_binding, only: c_null_ptr, c_int, c_double, c_associated, c_long + use fcvode_mod, only: CV_ONE_STEP, FCVode, FCVodeGetNumSteps class(Atmosphere), target, intent(inout) :: self character(:), allocatable, intent(out) :: err real(dp) :: tn integer(c_int) :: ierr + integer(c_long) :: nsteps_(1) + integer :: i, k real(c_double), parameter :: dum = 0.0_dp real(c_double) :: tcur(1) + type(PhotochemData), pointer :: dat + type(PhotochemVars), pointer :: var + type(PhotochemWrk), pointer :: wrk + + dat => self%dat + var => self%var + wrk => self%wrk - if (.not.c_associated(self%wrk%sun%cvode_mem)) then + if (.not.c_associated(wrk%sun%cvode_mem)) then err = "You must first initialize the stepper with 'initialize_stepper'" return endif - ierr = FCVode(self%wrk%sun%cvode_mem, dum, self%wrk%sun%sunvec_y, tcur, CV_ONE_STEP) + ierr = FCVode(wrk%sun%cvode_mem, dum, wrk%sun%sunvec_y, tcur, CV_ONE_STEP) if (ierr /= 0) then err = "CVODE step failed" return endif tn = tcur(1) + + ! Update nsteps + ierr = FCVodeGetNumSteps(wrk%sun%cvode_mem, nsteps_) + wrk%nsteps = nsteps_(1) + + ! Move over t and mix history + k = min(wrk%nsteps+1,size(wrk%t_history)) + do i = k,2,-1 + wrk%t_history(i) = wrk%t_history(i-1) + wrk%mix_history(:,:,i) = wrk%mix_history(:,:,i-1) + enddo + + ! Save current t and mix + wrk%t_history(1) = tn + block + real(c_double), pointer :: usol_tmp(:,:) + usol_tmp(1:dat%nq,1:var%nz) => wrk%sun%yvec(1:var%neqs) + call self%prep_atm_background_gas(usol_tmp, wrk%usol, wrk%molecules_per_particle) + wrk%mix_history(:,:,1) = wrk%usol + endblock + end function module subroutine destroy_stepper(self, err) diff --git a/src/atmosphere/photochem_atmosphere_rhs.f90 b/src/atmosphere/photochem_atmosphere_rhs.f90 index 30f5faf..6d73042 100644 --- a/src/atmosphere/photochem_atmosphere_rhs.f90 +++ b/src/atmosphere/photochem_atmosphere_rhs.f90 @@ -385,7 +385,7 @@ subroutine diffusion_coefficients(dat, var, den, mubar, & end subroutine - subroutine prep_atm_background_gas(self, usol_in, usol, molecules_per_particle) + module subroutine prep_atm_background_gas(self, usol_in, usol, molecules_per_particle) use photochem_common, only: molec_per_particle use photochem_const, only: small_real use photochem_enum, only: MixingRatioBC diff --git a/src/evoatmosphere/photochem_evoatmosphere.f90 b/src/evoatmosphere/photochem_evoatmosphere.f90 index ef9a42c..4106ef4 100644 --- a/src/evoatmosphere/photochem_evoatmosphere.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere.f90 @@ -50,6 +50,7 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo) !!! photochem_evoatmosphere_integrate.f90 !!! procedure :: evolve + procedure :: check_for_convergence procedure :: initialize_stepper procedure :: step procedure :: destroy_stepper @@ -166,6 +167,12 @@ module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, re ! Evolve atmosphere through time, and saves output in a binary Fortran file. end function + module function check_for_convergence(self, err) result(converged) + class(EvoAtmosphere), target, intent(inout) :: self + character(:), allocatable, intent(out) :: err + logical :: converged + end function + !> Initializes an integration starting at `usol_start` module subroutine initialize_stepper(self, usol_start, err) class(EvoAtmosphere), target, intent(inout) :: self diff --git a/src/evoatmosphere/photochem_evoatmosphere_integrate.f90 b/src/evoatmosphere/photochem_evoatmosphere_integrate.f90 index 19f738a..32e1beb 100644 --- a/src/evoatmosphere/photochem_evoatmosphere_integrate.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere_integrate.f90 @@ -686,6 +686,71 @@ subroutine read_end_of_evo_file(self, filename, t_eval, restart_index, tcur, top end subroutine + module function check_for_convergence(self, err) result(converged) + use, intrinsic :: iso_c_binding + class(EvoAtmosphere), target, intent(inout) :: self + character(:), allocatable, intent(out) :: err + logical :: converged + + integer :: i,j,ind + type(PhotochemData), pointer :: dat + type(PhotochemVars), pointer :: var + type(PhotochemWrkEvo), pointer :: wrk + + dat => self%dat + var => self%var + wrk => self%wrk + + converged = .false. + if (.not.c_associated(wrk%sun%cvode_mem)) then + err = "You must first initialize the stepper with 'initialize_stepper'" + return + endif + + ! If we reach equilibrium time, then converged + if (wrk%tn > var%equilibrium_time) then + converged = .true. + return + endif + + ! Now consider step history. + + ! Can't do analysis on step 0 + if (wrk%nsteps == 0) return + + ! Find index in history closest to time of interest. Note that this will only + ! consider a limited step history. We can not save all history. + ind = minloc(abs(wrk%t_history - var%conv_hist_factor*wrk%t_history(1)),1) + + ! Can't be current time, so we will check the previous step if needed. + if (ind == 1) ind = 2 + + ! Compute difference between current mixing ratios, and mixing ratios + ! at our index of interest. + do j = 1,var%nz + do i = 1,dat%nq + if (wrk%mix_history(i,j,1) > var%conv_min_mix) then + wrk%dmix(i,j) = abs(wrk%mix_history(i,j,1) - wrk%mix_history(i,j,ind)) + else + ! Ignore small mixing ratios + wrk%dmix(i,j) = 0.0_dp + endif + enddo + enddo + + ! Maximum normalized change + wrk%longdy = maxval(abs(wrk%dmix/wrk%mix_history(:,:,1))) + ! Also consider that change over time + wrk%longdydt = wrk%longdy/(wrk%t_history(1) - wrk%t_history(ind)) + + ! Check for convergence + if (wrk%longdy < var%conv_longdy .and. wrk%longdydt < var%conv_longdydt) then + converged = .true. + return + endif + + end function + module subroutine initialize_stepper(self, usol_start, err) use, intrinsic :: iso_c_binding use fcvode_mod, only: CV_BDF, CV_NORMAL, CV_ONE_STEP, FCVodeInit, FCVodeSVtolerances, & @@ -776,6 +841,13 @@ module subroutine initialize_stepper(self, usol_start, err) enddo enddo + ! Load initial conditions into history vector + wrk%nsteps = 0 + wrk%t_history = -1.0_dp + wrk%t_history(1) = 0.0_dp + wrk%mix_history = -1.0_dp + wrk%mix_history(:,:,1) = wrk%mix ! set by self%set_trop_ind + wrk%sun%abstol_nvec => FN_VMake_Serial(neqs_long, wrk%sun%abstol) if (.not. associated(wrk%sun%abstol_nvec)) then err = "CVODE setup error." @@ -865,15 +937,24 @@ module subroutine initialize_stepper(self, usol_start, err) end subroutine module function step(self, err) result(tn) - use iso_c_binding, only: c_null_ptr, c_int, c_double, c_associated - use fcvode_mod, only: CV_ONE_STEP, FCVode + use iso_c_binding, only: c_null_ptr, c_int, c_double, c_associated, c_long + use fcvode_mod, only: CV_ONE_STEP, FCVode, FCVodeGetNumSteps class(EvoAtmosphere), target, intent(inout) :: self character(:), allocatable, intent(out) :: err real(dp) :: tn integer(c_int) :: ierr + integer(c_long) :: nsteps_(1) + integer :: i, k real(c_double), parameter :: dum = 0.0_dp real(c_double) :: tcur(1) + type(PhotochemData), pointer :: dat + type(PhotochemVars), pointer :: var + type(PhotochemWrkEvo), pointer :: wrk + + dat => self%dat + var => self%var + wrk => self%wrk if (.not.c_associated(self%wrk%sun%cvode_mem)) then err = "You must first initialize the stepper with 'initialize_stepper'" @@ -886,6 +967,30 @@ module function step(self, err) result(tn) return endif tn = tcur(1) + + ! Update nsteps + ierr = FCVodeGetNumSteps(wrk%sun%cvode_mem, nsteps_) + wrk%nsteps = nsteps_(1) + + ! Move over t and mix history + k = min(wrk%nsteps+1,size(wrk%t_history)) + do i = k,2,-1 + wrk%t_history(i) = wrk%t_history(i-1) + wrk%mix_history(:,:,i) = wrk%mix_history(:,:,i-1) + enddo + + ! Save current t and mix + wrk%t_history(1) = tn + block + real(c_double), pointer :: usol_tmp(:,:) + usol_tmp(1:dat%nq,1:var%nz) => wrk%sun%yvec + call self%prep_atm_evo_gas(usol_tmp, wrk%usol, & + wrk%molecules_per_particle, wrk%pressure, wrk%density, wrk%mix, wrk%mubar, & + wrk%pressure_hydro, wrk%density_hydro, err) + if (allocated(err)) return + wrk%mix_history(:,:,1) = wrk%mix + endblock + end function module subroutine destroy_stepper(self, err) diff --git a/src/photochem_const.f90 b/src/photochem_const.f90 index 80366e7..90ee404 100644 --- a/src/photochem_const.f90 +++ b/src/photochem_const.f90 @@ -22,4 +22,7 @@ module photochem_const real(dp), parameter :: small_real = tiny(1.0_dp)**0.25_dp real(dp), parameter :: ln_small_real = log(small_real) + !> number of steps to save during integration, for steady-state checking + integer, parameter :: nsteps_save = 500 + end module diff --git a/src/photochem_types.f90 b/src/photochem_types.f90 index 087fb7d..d80c8eb 100644 --- a/src/photochem_types.f90 +++ b/src/photochem_types.f90 @@ -449,7 +449,17 @@ subroutine time_dependent_rate_fcn(tn, nz, rate) real(c_double) :: atol = 1.d-27 !! integration absolute tolerance integer :: mxsteps = 10000 !! max number of steps before integrator will give up. !> seconds. atomsphere considered in equilibrium if integrations reaches this time. - real(dp) :: equilibrium_time = 1.e17_dp + real(dp) :: equilibrium_time = 1.0e17_dp + !> For convergence checking. Considers mixing ratio change between t_now and time + !> t = t_now*conv_hist_factor to see if atmosphere is changing. + real(dp) :: conv_hist_factor = 0.5_dp + !> Minimum mixing ratio considered in convergence checking. + real(dp) :: conv_min_mix = 1.0e-20 + !> Threshold normalized change in mixing ratios for converchecking check. + real(dp) :: conv_longdy = 0.01_dp + !> Threshold normalized change in mixing ratios per time change for + !> convergence checking. + real(dp) :: conv_longdydt = 1.0e-4 real(c_double) :: initial_dt = 1.0e-6_dp !! intial timestep size (seconds) integer(c_int) :: max_err_test_failures = 15 !! CVODE max error test failures integer(c_int) :: max_order = 5 !! CVODE max order for BDF method. @@ -487,8 +497,26 @@ subroutine time_dependent_rate_fcn(tn, nz, rate) ! through the course of an integration ! used in cvode - integer(c_long) :: nsteps_previous = -10 + integer(c_long) :: nsteps_previous = -10 !! For printing type(SundialsData) :: sun !! CVODE data + ! All for determining convergence + integer :: nsteps = 0 !! Number of integration steps excuted. Updated + !! after every successful step. + !> History of times at previous integration steps. Index 1 is current, + !> while index 2, 3, 4 are previous steps. Updated after every successful step. + real(dp), allocatable :: t_history(:) + !> History of mixing ratios at previous integration steps. Index 1 is + !> current, while index 2, 3, 4 are previous steps. . Updated after + !> every successful step. + real(dp), allocatable :: mix_history(:,:,:) + !> Change in mixing ratio over some number of integrations steps. Updated + !> during convergence checking. + real(dp), allocatable :: dmix(:,:) + !> Normalized change in mixing ratios over some number of integrations steps + real(dp) :: longdy = 0.0_dp + !> Normalized change in mixing ratios divided by change in time + !> over some number of integrations steps. + real(dp) :: longdydt = 0.0_dp !> The current time (seconds). The is updated with each call to the !> right hand side, and jacobian. It is only important if @@ -566,10 +594,14 @@ subroutine init_PhotochemWrkEvo(self, nsp, np, nq, nz, nrT, kj, nw) end subroutine subroutine init_PhotochemWrk(self, nsp, np, nq, nz, nrT, kj, nw) + use photochem_const, only: nsteps_save class(PhotochemWrk), intent(inout) :: self integer, intent(in) :: nsp, np, nq, nz, nrT, kj, nw if (allocated(self%usol)) then + deallocate(self%t_history) + deallocate(self%mix_history) + deallocate(self%dmix) deallocate(self%usol) deallocate(self%mubar) deallocate(self%pressure) @@ -599,6 +631,9 @@ subroutine init_PhotochemWrk(self, nsp, np, nq, nz, nrT, kj, nw) deallocate(self%rainout_rates) endif + allocate(self%t_history(nsteps_save)) + allocate(self%mix_history(nq,nz,nsteps_save)) + allocate(self%dmix(nq,nz)) allocate(self%usol(nq,nz)) allocate(self%mubar(nz)) allocate(self%pressure(nz)) diff --git a/tests/testevo.f90 b/tests/testevo.f90 index 50c0e1e..9ba02c3 100644 --- a/tests/testevo.f90 +++ b/tests/testevo.f90 @@ -60,6 +60,7 @@ subroutine test_stepper() type(EvoAtmosphere) :: pc integer :: i real(dp) :: tn + logical :: converged pc = EvoAtmosphere("../photochem/data/reaction_mechanisms/zahnle_earth.yaml", & "../tests/testevo_settings2.yaml", & @@ -84,6 +85,12 @@ subroutine test_stepper() print*,trim(err) stop 1 endif + converged = pc%check_for_convergence(err) + if (allocated(err)) then + print*,trim(err) + stop 1 + endif + if (converged) exit enddo call pc%destroy_stepper(err)