From c76ff953831a1343e48a82fe1e297426d78f05ea Mon Sep 17 00:00:00 2001
From: Nick Wogan <nicholaswogan@gmail.com>
Date: Sun, 3 Nov 2024 18:18:40 -0700
Subject: [PATCH] find_steady_state function

---
 photochem/cython/EvoAtmosphere.pyx            |  62 ++++++++++
 photochem/cython/EvoAtmosphere_pxd.pxd        |   4 +
 photochem/cython/PhotochemVars.pyx            |  37 +++++-
 photochem/cython/PhotochemVars_pxd.pxd        |  14 ++-
 photochem/cython/PhotochemWrk.pyx             |   7 ++
 photochem/cython/PhotochemWrk_pxd.pxd         |   2 +
 photochem/cython/_photochem.pyx               |   1 +
 photochem/fortran/EvoAtmosphere_wrapper.f90   |  41 +++++++
 photochem/fortran/PhotochemVars_wrapper.f90   |  64 ++++++++++
 photochem/fortran/PhotochemWrk_wrapper.f90    |  14 ++-
 src/evoatmosphere/photochem_evoatmosphere.f90 |  26 ++++
 .../photochem_evoatmosphere_integrate.f90     | 115 ++++++++++++++++++
 src/photochem_types.f90                       |  26 +++-
 tests/memtest_evo.f90                         |  27 ++++
 tests/testevo.f90                             |  26 +---
 15 files changed, 433 insertions(+), 33 deletions(-)

diff --git a/photochem/cython/EvoAtmosphere.pyx b/photochem/cython/EvoAtmosphere.pyx
index 6c9c4fe..462dffe 100644
--- a/photochem/cython/EvoAtmosphere.pyx
+++ b/photochem/cython/EvoAtmosphere.pyx
@@ -527,6 +527,68 @@ cdef class EvoAtmosphere:
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
 
+  def initialize_robust_stepper(self, ndarray[double, ndim=2] usol_start):
+    """Initializes a robust integration starting at `usol_start`.
+
+    Parameters
+    ----------
+    usol_start : ndarray[double,ndim=2]
+        Initial number densities (molecules/cm^3)
+    """   
+    cdef char err[ERR_LEN+1]
+    cdef int nq = self.dat.nq
+    cdef int nz = self.var.nz
+    cdef ndarray usol_start_ = np.asfortranarray(usol_start)
+    if usol_start_.shape[0] != nq or usol_start_.shape[1] != nz:
+      raise PhotoException("Input usol_start is the wrong size.")
+    ea_pxd.evoatmosphere_initialize_robust_stepper_wrapper(self._ptr, &nq, &nz,  <double *>usol_start_.data, err)
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
+    
+  def robust_step(self):
+    """Takes one internal robust integration step. Function `initialize_robust_stepper`
+    must have been called before this.
+
+    Returns
+    -------
+    give_up : bool
+         If True, then the algorithm thinks it is time to give up.
+    converged : bool
+        If True, then the integration has converged to a steady state.
+    """
+    cdef bool give_up, converged
+    cdef char err[ERR_LEN+1]
+    ea_pxd.evoatmosphere_robust_step_wrapper(self._ptr, &give_up, &converged, err)
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
+    return give_up, converged
+
+  def find_steady_state(self):
+    """Integrates using a robust stepper until a steady state has been achieved.
+
+    Returns
+    -------
+    converged : bool
+        If True, then the integration has converged to a steady state.
+    """
+
+    cdef bool converged, give_up
+    converged = False
+
+    self.initialize_robust_stepper(self.wrk.usol)
+
+    while True:
+      give_up, converged = self.robust_step()
+
+      if give_up:
+        converged = False
+        return converged
+      
+      if converged:
+        return converged
+
+      PyErr_CheckSignals()
+
   def production_and_loss(self, str species, ndarray[double, ndim=2] usol):
     """Computes the production and loss of input `species`.
     See ProductionLoss object in photochem_types.f90.
diff --git a/photochem/cython/EvoAtmosphere_pxd.pxd b/photochem/cython/EvoAtmosphere_pxd.pxd
index f55680a..452b2d3 100644
--- a/photochem/cython/EvoAtmosphere_pxd.pxd
+++ b/photochem/cython/EvoAtmosphere_pxd.pxd
@@ -57,6 +57,10 @@ cdef extern double evoatmosphere_step_wrapper(EvoAtmosphere *ptr, char *err)
 
 cdef extern void evoatmosphere_destroy_stepper_wrapper(EvoAtmosphere *ptr, char *err)
 
+cdef extern void evoatmosphere_initialize_robust_stepper_wrapper(EvoAtmosphere *ptr, int *nq, int *nz, double *usol_start, char *err)
+
+cdef extern void evoatmosphere_robust_step_wrapper(EvoAtmosphere *ptr, bool *give_up, bool *converged, char *err)
+
 cdef extern void evoatmosphere_production_and_loss_wrapper(EvoAtmosphere *ptr, char *species, int *nq, 
                                                         int *nz, double *usol, pl_pxd.ProductionLoss **pl_ptr, char *err)
 
diff --git a/photochem/cython/PhotochemVars.pyx b/photochem/cython/PhotochemVars.pyx
index ac953dd..9bc1e42 100644
--- a/photochem/cython/PhotochemVars.pyx
+++ b/photochem/cython/PhotochemVars.pyx
@@ -357,5 +357,38 @@ cdef class PhotochemVars:
     def __set__(self, bool val):
       var_pxd.photochemvars_upwind_molec_diff_set(self._ptr, &val)
 
-    
-  
\ No newline at end of file
+  property nerrors_before_giveup:
+    "int. Number of integration errors before giving up completely"
+    def __get__(self):
+      cdef int val
+      var_pxd.photochemvars_nerrors_before_giveup_get(self._ptr, &val)
+      return val
+    def __set__(self, int val):
+      var_pxd.photochemvars_nerrors_before_giveup_set(self._ptr, &val)
+
+  property nsteps_before_conv_check:
+    "int. Number of steps to take before checking for convergence"
+    def __get__(self):
+      cdef int val
+      var_pxd.photochemvars_nsteps_before_conv_check_get(self._ptr, &val)
+      return val
+    def __set__(self, int val):
+      var_pxd.photochemvars_nsteps_before_conv_check_set(self._ptr, &val)
+
+  property nsteps_before_reinit:
+    "int. Number of steps before reinitializing the integration"
+    def __get__(self):
+      cdef int val
+      var_pxd.photochemvars_nsteps_before_reinit_get(self._ptr, &val)
+      return val
+    def __set__(self, int val):
+      var_pxd.photochemvars_nsteps_before_reinit_set(self._ptr, &val)
+
+  property nsteps_before_giveup:
+    "int. Number of total steps to take before giving up."
+    def __get__(self):
+      cdef int val
+      var_pxd.photochemvars_nsteps_before_giveup_get(self._ptr, &val)
+      return val
+    def __set__(self, int val):
+      var_pxd.photochemvars_nsteps_before_giveup_set(self._ptr, &val)
\ No newline at end of file
diff --git a/photochem/cython/PhotochemVars_pxd.pxd b/photochem/cython/PhotochemVars_pxd.pxd
index 1fd176e..3a6755e 100644
--- a/photochem/cython/PhotochemVars_pxd.pxd
+++ b/photochem/cython/PhotochemVars_pxd.pxd
@@ -101,4 +101,16 @@ cdef extern void photochemvars_fast_arbitrary_rate_get(PhotochemVars *ptr, doubl
 cdef extern void photochemvars_fast_arbitrary_rate_set(PhotochemVars *ptr, double *val)
 
 cdef extern void photochemvars_upwind_molec_diff_get(PhotochemVars *ptr, bool *val)
-cdef extern void photochemvars_upwind_molec_diff_set(PhotochemVars *ptr, bool *val)
\ No newline at end of file
+cdef extern void photochemvars_upwind_molec_diff_set(PhotochemVars *ptr, bool *val)
+
+cdef extern void photochemvars_nerrors_before_giveup_get(PhotochemVars *ptr, int *val)
+cdef extern void photochemvars_nerrors_before_giveup_set(PhotochemVars *ptr, int *val)
+
+cdef extern void photochemvars_nsteps_before_conv_check_get(PhotochemVars *ptr, int *val)
+cdef extern void photochemvars_nsteps_before_conv_check_set(PhotochemVars *ptr, int *val)
+
+cdef extern void photochemvars_nsteps_before_reinit_get(PhotochemVars *ptr, int *val)
+cdef extern void photochemvars_nsteps_before_reinit_set(PhotochemVars *ptr, int *val)
+
+cdef extern void photochemvars_nsteps_before_giveup_get(PhotochemVars *ptr, int *val)
+cdef extern void photochemvars_nsteps_before_giveup_set(PhotochemVars *ptr, int *val)
\ No newline at end of file
diff --git a/photochem/cython/PhotochemWrk.pyx b/photochem/cython/PhotochemWrk.pyx
index 18dc3c6..8f1cf9b 100644
--- a/photochem/cython/PhotochemWrk.pyx
+++ b/photochem/cython/PhotochemWrk.pyx
@@ -23,6 +23,13 @@ cdef class PhotochemWrk:
   def __cinit__(self):
     self._ptr = NULL
 
+  property nsteps_total:
+    "int. Total number of steps in a robust integration."
+    def __get__(self):
+      cdef int val
+      wrk_pxd.photochemwrk_nsteps_total_get(self._ptr, &val)
+      return val
+
   property nsteps:
     "int. Number of integration steps excuted. Updated after every successful step."
     def __get__(self):
diff --git a/photochem/cython/PhotochemWrk_pxd.pxd b/photochem/cython/PhotochemWrk_pxd.pxd
index 75519cd..edf3e70 100644
--- a/photochem/cython/PhotochemWrk_pxd.pxd
+++ b/photochem/cython/PhotochemWrk_pxd.pxd
@@ -8,6 +8,8 @@ cdef extern from "PhotochemWrk.h":
   struct PhotochemWrkEvo:
     pass
 
+cdef extern void photochemwrk_nsteps_total_get(PhotochemWrk *ptr, int *val)
+
 cdef extern void photochemwrk_nsteps_get(PhotochemWrk *ptr, int *val)
 
 cdef extern void photochemwrk_t_history_get_size(PhotochemWrk *ptr, int *dim1)
diff --git a/photochem/cython/_photochem.pyx b/photochem/cython/_photochem.pyx
index db94ac5..b83dfd3 100644
--- a/photochem/cython/_photochem.pyx
+++ b/photochem/cython/_photochem.pyx
@@ -3,6 +3,7 @@ from libcpp cimport bool
 from libc.stdlib cimport malloc, free
 from libc.stdint cimport uintptr_t
 from cpython.object cimport PyObject_GenericSetAttr
+from cpython.exc cimport PyErr_CheckSignals
 import numpy as np
 import ctypes as ct
 import os
diff --git a/photochem/fortran/EvoAtmosphere_wrapper.f90 b/photochem/fortran/EvoAtmosphere_wrapper.f90
index 9e94440..851f900 100644
--- a/photochem/fortran/EvoAtmosphere_wrapper.f90
+++ b/photochem/fortran/EvoAtmosphere_wrapper.f90
@@ -478,6 +478,47 @@ subroutine evoatmosphere_destroy_stepper_wrapper(ptr, err) bind(c)
     endif
   end subroutine
 
+  subroutine evoatmosphere_initialize_robust_stepper_wrapper(ptr, nq, nz, usol_start, err) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: nq, nz
+    real(c_double), intent(in) :: usol_start(nq, nz)
+    character(len=c_char), intent(out) :: err(err_len+1)
+    
+    character(:), allocatable :: err_f
+  
+    type(EvoAtmosphere), pointer :: pc
+    call c_f_pointer(ptr, pc)
+    
+    call pc%initialize_robust_stepper(usol_start, err_f)
+    err(1) = c_null_char
+    if (allocated(err_f)) then
+      call copy_string_ftoc(err_f, err)
+    endif
+  end subroutine
+
+  subroutine evoatmosphere_robust_step_wrapper(ptr, give_up, converged, err) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    logical(c_bool), intent(out) :: give_up, converged
+    character(len=c_char), intent(out) :: err(err_len+1)
+    
+    logical :: give_up_f, converged_f
+    character(:), allocatable :: err_f
+    type(EvoAtmosphere), pointer :: pc
+
+    call c_f_pointer(ptr, pc)
+  
+    call pc%robust_step(give_up_f, converged_f, err_f)
+
+    give_up = give_up_f
+    converged = converged_f
+
+    err(1) = c_null_char
+    if (allocated(err_f)) then
+      call copy_string_ftoc(err_f, err)
+    endif
+
+  end subroutine
+
   subroutine evoatmosphere_production_and_loss_wrapper(ptr, species, nq, nz, usol, pl_ptr, err) bind(c)
     use photochem, only: ProductionLoss
     type(c_ptr), value, intent(in) :: ptr
diff --git a/photochem/fortran/PhotochemVars_wrapper.f90 b/photochem/fortran/PhotochemVars_wrapper.f90
index 5105c85..4cbddd7 100644
--- a/photochem/fortran/PhotochemVars_wrapper.f90
+++ b/photochem/fortran/PhotochemVars_wrapper.f90
@@ -503,4 +503,68 @@ subroutine photochemvars_upwind_molec_diff_set(ptr, val) bind(c)
     call c_f_pointer(ptr, var)
     var%upwind_molec_diff = val
   end subroutine
+
+  subroutine photochemvars_nerrors_before_giveup_get(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%nerrors_before_giveup
+  end subroutine
+  
+  subroutine photochemvars_nerrors_before_giveup_set(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%nerrors_before_giveup = val
+  end subroutine
+
+  subroutine photochemvars_nsteps_before_conv_check_get(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%nsteps_before_conv_check
+  end subroutine
+  
+  subroutine photochemvars_nsteps_before_conv_check_set(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%nsteps_before_conv_check = val
+  end subroutine
+
+  subroutine photochemvars_nsteps_before_reinit_get(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%nsteps_before_reinit
+  end subroutine
+  
+  subroutine photochemvars_nsteps_before_reinit_set(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%nsteps_before_reinit = val
+  end subroutine
+
+  subroutine photochemvars_nsteps_before_giveup_get(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%nsteps_before_giveup
+  end subroutine
+  
+  subroutine photochemvars_nsteps_before_giveup_set(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%nsteps_before_giveup = val
+  end subroutine
   
\ No newline at end of file
diff --git a/photochem/fortran/PhotochemWrk_wrapper.f90 b/photochem/fortran/PhotochemWrk_wrapper.f90
index d921009..e8bdf2d 100644
--- a/photochem/fortran/PhotochemWrk_wrapper.f90
+++ b/photochem/fortran/PhotochemWrk_wrapper.f90
@@ -3,12 +3,20 @@
   !!! getters and setters !!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+  subroutine photochemwrk_nsteps_total_get(ptr, val) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: val
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    val = wrk%nsteps_total
+  end subroutine
+
   subroutine photochemwrk_nsteps_get(ptr, val) bind(c)
     type(c_ptr), value, intent(in) :: ptr
     integer(c_int), intent(out) :: val
-    type(PhotochemWrk), pointer :: var
-    call c_f_pointer(ptr, var)
-    val = var%nsteps
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    val = wrk%nsteps
   end subroutine
 
   subroutine photochemwrk_t_history_get_size(ptr, dim1) bind(c)
diff --git a/src/evoatmosphere/photochem_evoatmosphere.f90 b/src/evoatmosphere/photochem_evoatmosphere.f90
index 51a5832..4bec504 100644
--- a/src/evoatmosphere/photochem_evoatmosphere.f90
+++ b/src/evoatmosphere/photochem_evoatmosphere.f90
@@ -54,6 +54,9 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
     procedure :: initialize_stepper
     procedure :: step
     procedure :: destroy_stepper
+    procedure :: initialize_robust_stepper
+    procedure :: robust_step
+    procedure :: find_steady_state
 
     !~~ photochem_evoatmosphere_utils.f90 ~~!
     procedure :: out2atmosphere_txt
@@ -202,6 +205,29 @@ module subroutine destroy_stepper(self, err)
       class(EvoAtmosphere), target, intent(inout) :: self
       character(:), allocatable, intent(out) :: err
     end subroutine
+
+    !> Initializes a robust integration starting at `usol_start`
+    module subroutine initialize_robust_stepper(self, usol_start, err)
+      class(EvoAtmosphere), target, intent(inout) :: self
+      real(dp), intent(in) :: usol_start(:,:) !! Initial number densities (molecules/cm^3)
+      character(:), allocatable, intent(out) :: err
+    end subroutine
+    
+    !> Takes one robust integration step. Function `initialize_robust_stepper`
+    !> must have been called before this.
+    module subroutine robust_step(self, give_up, converged, err)
+      class(EvoAtmosphere), target, intent(inout) :: self
+      logical, intent(out) :: give_up !! If .true., then the algorithm thinks it is time to give up.
+      logical, intent(out) :: converged !! If .true., then the integration has converged to a steady state.
+      character(:), allocatable, intent(out) :: err  
+    end subroutine
+    
+    !> Integrates using a robust stepper until a steady state has been achieved.
+    module function find_steady_state(self, err) result(converged)
+      class(EvoAtmosphere), target, intent(inout) :: self
+      character(:), allocatable, intent(out) :: err
+      logical :: converged !! If .true., then the integration has converged to a steady state.
+    end function
     
     module function RhsFn_evo(tn, sunvec_y, sunvec_f, user_data) &
                           result(ierr) bind(c, name='RhsFn_evo')
diff --git a/src/evoatmosphere/photochem_evoatmosphere_integrate.f90 b/src/evoatmosphere/photochem_evoatmosphere_integrate.f90
index 4c65d54..ebb52dc 100644
--- a/src/evoatmosphere/photochem_evoatmosphere_integrate.f90
+++ b/src/evoatmosphere/photochem_evoatmosphere_integrate.f90
@@ -1019,5 +1019,120 @@ module subroutine destroy_stepper(self, err)
     if (allocated(err)) return
     
   end subroutine
+
+  module subroutine initialize_robust_stepper(self, usol_start, err)
+    class(EvoAtmosphere), target, intent(inout) :: self
+    real(dp), intent(in) :: usol_start(:,:)
+    character(:), allocatable, intent(out) :: err
+
+    self%wrk%nsteps_total = 0
+    self%wrk%nerrors_total = 0
+    call self%initialize_stepper(usol_start, err)
+    if (allocated(err)) return
+
+  end subroutine
+
+  module subroutine robust_step(self, give_up, converged, err)
+    class(EvoAtmosphere), target, intent(inout) :: self
+    logical, intent(out) :: give_up
+    logical, intent(out) :: converged
+    character(:), allocatable, intent(out) :: err
+
+    real(dp) :: tn
+
+    type(PhotochemVars), pointer :: var
+    type(PhotochemWrkEvo), pointer :: wrk
+    
+    var => self%var
+    wrk => self%wrk
+
+    converged = .false.
+    give_up = .false.
+
+    if (wrk%nsteps_total < 0) then
+      err = "You must first initialize a robust stepper with 'initialize_robust_stepper'"
+      return
+    endif
+    if (var%nsteps_before_conv_check >= var%nsteps_before_reinit) then
+      err = "`nsteps_before_conv_check` should be < `nsteps_before_reinit`"
+      return
+    endif
+
+    tn = self%step(err)
+    if (.not.allocated(err)) then
+      ! If step worked, then we add it to counter
+      wrk%nsteps_total = wrk%nsteps_total + 1
+    else
+      ! There was an error
+      deallocate(err)
+      wrk%nerrors_total = wrk%nerrors_total + 1
+
+      ! If there are too many errors, then give up
+      if (wrk%nerrors_total > var%nerrors_before_giveup) then
+        give_up = .true.
+        return
+      endif
+
+      ! Trim negative numbers, and reinitialize
+      wrk%usol = max(wrk%usol, var%reinit_min_density)
+      call self%initialize_stepper(wrk%usol, err)
+      if (allocated(err)) return
+
+    endif
+
+    ! If we have reached the equilibrium time, then we have converged
+    if (tn > var%equilibrium_time) then
+      converged = .true.
+      return
+    endif
+
+    ! We allow convergence via other criteria, but only after
+    ! a minimum number of steps has been performed
+    if (self%wrk%nsteps > var%nsteps_before_conv_check) then
+      converged = self%check_for_convergence(err)
+      if (allocated(err)) return
+      if (converged) return
+    endif
+
+    ! Reinitialize integrator after some number of steps
+    if (self%wrk%nsteps > var%nsteps_before_reinit) then
+      wrk%usol = max(wrk%usol, var%reinit_min_density)
+      call self%initialize_stepper(wrk%usol, err)
+      if (allocated(err)) return
+    endif
+
+    ! Give up after a large number of steps
+    if (wrk%nsteps_total > var%nsteps_before_giveup) then
+      give_up = .true.
+      return
+    endif
+
+  end subroutine
+
+  module function find_steady_state(self, err) result(converged)
+    class(EvoAtmosphere), target, intent(inout) :: self
+    character(:), allocatable, intent(out) :: err
+    logical :: converged
+
+    logical :: give_up
+
+    converged = .false.
+
+    call self%initialize_robust_stepper(self%wrk%usol, err)
+    if (allocated(err)) return
+
+    do
+      call self%robust_step(give_up, converged, err)
+      if (allocated(err)) return
+
+      if (give_up) then
+        converged = .false.
+        return
+      endif
+
+      if (converged) return
+    enddo
+
+  end function
   
 end submodule
\ No newline at end of file
diff --git a/src/photochem_types.f90 b/src/photochem_types.f90
index aadbda5..4bd8538 100644
--- a/src/photochem_types.f90
+++ b/src/photochem_types.f90
@@ -456,12 +456,13 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
     ! other 
     !> number of times we initialize CVODE when it returns
     !> a potentially recoverable error. ONLY USED IN EVOATMOSPHERE (NOT ATMOSPHERE)
+    !> in the `evolve` method.
     integer :: max_error_reinit_attempts = 2 
     real(c_double) :: rtol = 1.0e-3_dp !! integration relative tolerance
     !> Integration absolute tolerance. If autodiff == .true., then the model
     !> works better when atol is smaller (e.g., atol = ~1.0e-18).
-    real(c_double) :: atol = 1.0e-25_dp 
-    integer :: mxsteps = 10000 !! max number of steps before integrator will give up.
+    real(c_double) :: atol = 1.0e-23_dp 
+    integer :: mxsteps = 100000 !! max number of steps before integrator will give up.
     !> seconds. atomsphere considered in equilibrium if integrations reaches this time.
     real(dp) :: equilibrium_time = 1.0e17_dp
     !> For convergence checking. Considers mixing ratio change between t_now and time 
@@ -490,6 +491,21 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
     !> diffusion terms instead of a centered scheme. This permits stability (at the cost 
     !> of accuracy) for atmospheres with strong molcular advection in the upper atmosphere.
     logical :: upwind_molec_diff = .false.
+
+    ! Settings for the `robust_stepper` and `find_steady_state` methods
+    !> Number of integration errors before giving up completely
+    integer :: nerrors_before_giveup = 10
+    !> Number of steps to take before checking for convergence
+    integer :: nsteps_before_conv_check = 300
+    !> Number of steps before reinitializing the integration
+    integer :: nsteps_before_reinit = 1000
+    !> Number of total steps to take before giving up.
+    integer :: nsteps_before_giveup = 100000
+    !> When the integrator reinitializes, this is the minimum
+    !> density allowed (molecules/cm^3)
+    real(dp) :: reinit_min_density = 1.0e-40_dp
+    ! End settings for `robust_stepper` and `find_steady_state`
+
   end type
 
   type :: SundialsData
@@ -518,6 +534,12 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
   type :: PhotochemWrk
     ! PhotochemWrk are work variables that change
     ! through the course of an integration
+
+    ! Total step counter for robust_step method
+    !> Total number of steps in a robust integration.
+    integer :: nsteps_total = -1
+    !> Total number of errors experienced in the robust integration.
+    integer :: nerrors_total = -1
     
     ! used in cvode
     integer(c_long) :: nsteps_previous = -10 !! For printing
diff --git a/tests/memtest_evo.f90 b/tests/memtest_evo.f90
index 9dbf675..38fd546 100644
--- a/tests/memtest_evo.f90
+++ b/tests/memtest_evo.f90
@@ -36,6 +36,7 @@ subroutine test_methods(filename)
     ! initialize_stepper : test_step
     call test_step(pcs)
     ! destroy_stepper : test_step
+    call test_robust_step(pcs)
     call test_out2atmosphere(pcs)
     call test_gas_fluxes(pcs) 
     call test_set_lower_bc(pcs)
@@ -358,6 +359,32 @@ subroutine test_step(pc)
     
   end subroutine
 
+  subroutine test_robust_step(pc)
+    type(EvoAtmosphere), intent(inout) :: pc
+    
+    character(:), allocatable :: err
+    logical :: give_up, converged
+    
+    call pc%initialize_robust_stepper(pc%var%usol_init, err)
+    if (allocated(err)) then
+      print*,trim(err)
+      stop 1
+    endif
+    
+    call pc%robust_step(give_up, converged, err)
+    if (allocated(err)) then
+      print*,trim(err)
+      stop 1
+    endif
+
+    call pc%robust_step(give_up, converged, err)
+    if (allocated(err)) then
+      print*,trim(err)
+      stop 1
+    endif
+    
+  end subroutine
+
   subroutine test_evolve(pc)
     use futils, only: linspace
     type(EvoAtmosphere), intent(inout) :: pc
diff --git a/tests/testevo.f90 b/tests/testevo.f90
index b2b464d..6e267d2 100644
--- a/tests/testevo.f90
+++ b/tests/testevo.f90
@@ -8,7 +8,6 @@ subroutine main_()
     character(:), allocatable :: err
     type(EvoAtmosphere) :: pc
     integer :: ind
-    real(dp) :: tn
     logical :: converged
 
     ! Print version
@@ -27,30 +26,7 @@ subroutine main_()
       stop 1
     endif
     
-    ! Initialize stepper
-    call pc%initialize_stepper(pc%var%usol_init, err)
-    if (allocated(err)) then
-      print*,trim(err)
-      stop 1
-    endif
-
-    ! Integrate to photochemical equilibrium
-    do
-      tn = pc%step(err)
-      if (allocated(err)) then
-        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
-
-    ! Destroy stepper
-    call pc%destroy_stepper(err)
+    converged = pc%find_steady_state(err)
     if (allocated(err)) then
       print*,trim(err)
       stop 1