find_steady_state function

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.

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)
 

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)
 
-
-
+  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)

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)
+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)

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):

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)

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

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

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
 

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)

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')