python wrapper for convergence checking

photochem/cython/Atmosphere.pyx

@@ -73,7 +73,16 @@ cdef class Atmosphere:
       wrk = PhotochemWrk(alloc = False)
       wrk._ptr = self._wrk_ptr
       return wrk
-
+
+  def check_for_convergence(self):
+    "Determines if integration has converged to photochemical steady-state."
+    cdef bool converged
+    cdef char err[ERR_LEN+1]
+    a_pxd.atmosphere_check_for_convergence_wrapper(&self._ptr, &converged, err)
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
+    return converged
+
   def photochemical_equilibrium(self):
     "Integrates to photochemical equilibrium starting from `self.var.usol_init`"  
     cdef bool success

photochem/cython/Atmosphere_pxd.pxd

@@ -15,6 +15,7 @@ cdef extern void atmosphere_create_wrapper(void *ptr, char *mechanism_file,
                                         char *atmosphere_txt, char *data_dir, void *dat_ptr, void *var_ptr,
                                         void *wrk_ptr, char *err);
 
+cdef extern void atmosphere_check_for_convergence_wrapper(void *ptr, bool *converged, char* err)
 cdef extern void atmosphere_photochemical_equilibrium_wrapper(void *ptr, bool *success, char* err)
 
 cdef extern void atmosphere_out2atmosphere_txt_wrapper(void *ptr, char *filename, bool *overwrite, bool *clip, char *err)                         

photochem/cython/EvoAtmosphere.pyx

@@ -393,6 +393,15 @@ cdef class EvoAtmosphere:
       raise PhotoException(err.decode("utf-8").strip())
     return success
 
+  def check_for_convergence(self):
+    "Determines if integration has converged to photochemical steady-state."  
+    cdef bool converged
+    cdef char err[ERR_LEN+1]
+    ea_pxd.evoatmosphere_check_for_convergence_wrapper(&self._ptr, &converged, err)
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
+    return converged
+
   def initialize_stepper(self, ndarray[double, ndim=2] usol_start):
     """Initializes an integration starting at `usol_start`.
 

photochem/cython/EvoAtmosphere_pxd.pxd

@@ -34,6 +34,8 @@ cdef extern void evoatmosphere_evolve_wrapper(void *ptr, char *filename,
                 double *tstart, int *nq, int *nz, double *usol, 
                 int *nt, double *t_eval, bool *overwrite, bool *restart_from_file, bool *success, char *err)
 
+cdef extern void evoatmosphere_check_for_convergence_wrapper(void *ptr, bool *converged, char* err)
+
 cdef extern void evoatmosphere_initialize_stepper_wrapper(void *ptr, int *nq, int *nz, double *usol_start, char *err)
 
 cdef extern double evoatmosphere_step_wrapper(void *ptr, char *err)

photochem/cython/PhotochemVars.pyx

@@ -230,6 +230,46 @@ cdef class PhotochemVars:
       return val
     def __set__(self, double val):
       var_pxd.photochemvars_equilibrium_time_set(&self._ptr, &val)
+
+  property conv_hist_factor:
+    """double. For convergence checking. Considers mixing ratio change between t_now and time 
+    t = t_now*conv_hist_factor to see if atmosphere is changing.
+    """
+    def __get__(self):
+      cdef double val
+      var_pxd.photochemvars_conv_hist_factor_get(&self._ptr, &val)
+      return val
+    def __set__(self, double val):
+      var_pxd.photochemvars_conv_hist_factor_set(&self._ptr, &val)
+
+  property conv_min_mix:
+    "double. Minimum mixing ratio considered in convergence checking."
+    def __get__(self):
+      cdef double val
+      var_pxd.photochemvars_conv_min_mix_get(&self._ptr, &val)
+      return val
+    def __set__(self, double val):
+      var_pxd.photochemvars_conv_min_mix_set(&self._ptr, &val)
+
+  property conv_longdy:
+    "double. Threshold normalized change in mixing ratios for converchecking check."
+    def __get__(self):
+      cdef double val
+      var_pxd.photochemvars_conv_longdy_get(&self._ptr, &val)
+      return val
+    def __set__(self, double val):
+      var_pxd.photochemvars_conv_longdy_set(&self._ptr, &val)
+
+  property conv_longdydt:
+    """double. Threshold normalized change in mixing ratios per time change for 
+    convergence checking.
+    """
+    def __get__(self):
+      cdef double val
+      var_pxd.photochemvars_conv_longdydt_get(&self._ptr, &val)
+      return val
+    def __set__(self, double val):
+      var_pxd.photochemvars_conv_longdydt_set(&self._ptr, &val)
 
   property verbose:
     "int. 0 == no printing. 1 == some printing. 2 == bunch of printing."

photochem/cython/PhotochemVars_pxd.pxd

@@ -65,6 +65,18 @@ cdef extern void photochemvars_mxsteps_set(void *ptr, int *val)
 cdef extern void photochemvars_equilibrium_time_get(void *ptr, double *val)
 cdef extern void photochemvars_equilibrium_time_set(void *ptr, double *val)
 
+cdef extern void photochemvars_conv_hist_factor_get(void *ptr, double *val)
+cdef extern void photochemvars_conv_hist_factor_set(void *ptr, double *val)
+
+cdef extern void photochemvars_conv_min_mix_get(void *ptr, double *val)
+cdef extern void photochemvars_conv_min_mix_set(void *ptr, double *val)
+
+cdef extern void photochemvars_conv_longdy_get(void *ptr, double *val)
+cdef extern void photochemvars_conv_longdy_set(void *ptr, double *val)
+
+cdef extern void photochemvars_conv_longdydt_get(void *ptr, double *val)
+cdef extern void photochemvars_conv_longdydt_set(void *ptr, double *val)
+
 cdef extern void photochemvars_verbose_get(void *ptr, int *val)
 cdef extern void photochemvars_verbose_set(void *ptr, int *val)
 

photochem/cython/PhotochemWrk.pyx

@@ -33,6 +33,51 @@ cdef class PhotochemWrk:
       wrk_pxd.deallocate_photochemwrk(&self._ptr)
       self._ptr = NULL
 
+  property nsteps:
+    "int. Number of integration steps excuted. Updated after every successful step."
+    def __get__(self):
+      cdef int val
+      wrk_pxd.photochemwrk_nsteps_get(&self._ptr, &val)
+      return val
+
+  property t_history:
+    """ndarray[double,dim=1], shape (500). History of times at previous integration steps. 
+    Index 1 is current, while index 2, 3, 4 are previous steps. Updated after every successful step.
+    """
+    def __get__(self):
+      cdef int dim1
+      wrk_pxd.photochemwrk_t_history_get_size(&self._ptr, &dim1)
+      cdef ndarray arr = np.empty((dim1), np.double, order="F")
+      wrk_pxd.photochemwrk_t_history_get(&self._ptr, &dim1, <double *>arr.data)
+      return arr
+
+  property mix_history:
+    """ndarray[double,dim=3], shape (nq,nz,500). 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.
+    """
+    def __get__(self):
+      cdef int dim1,dim2,dim3
+      wrk_pxd.photochemwrk_mix_history_get_size(&self._ptr, &dim1, &dim2, &dim3)
+      cdef ndarray arr = np.empty((dim1,dim2,dim3), np.double, order="F")
+      wrk_pxd.photochemwrk_mix_history_get(&self._ptr, &dim1, &dim2, &dim3, <double *>arr.data)
+      return arr
+
+  property longdy:
+    "double. Normalized change in mixing ratios over some number of integrations steps."
+    def __get__(self):
+      cdef double val
+      wrk_pxd.photochemwrk_longdy_get(&self._ptr, &val)
+      return val
+
+  property longdydt:
+    """double. Normalized change in mixing ratios divided by change in time over some 
+    number of integrations steps.
+    """
+    def __get__(self):
+      cdef double val
+      wrk_pxd.photochemwrk_longdydt_get(&self._ptr, &val)
+      return val
+
   property tn:
     "double. The current time of integration."
     def __get__(self):

photochem/cython/PhotochemWrk_pxd.pxd

@@ -5,6 +5,18 @@ cdef extern from "<stdbool.h>":
 cdef extern void allocate_photochemwrk(void *ptr)
 cdef extern void deallocate_photochemwrk(void *ptr)
 
+cdef extern void photochemwrk_nsteps_get(void *ptr, int *val)
+
+cdef extern void photochemwrk_t_history_get_size(void *ptr, int *dim1)
+cdef extern void photochemwrk_t_history_get(void *ptr, int *dim1, double *arr)
+
+cdef extern void photochemwrk_mix_history_get_size(void *ptr, int *dim1, int *dim2, int *dim3)
+cdef extern void photochemwrk_mix_history_get(void *ptr, int *dim1, int *dim2, int *dim3, double *arr)
+
+cdef extern void photochemwrk_longdy_get(void *ptr, double *val)
+
+cdef extern void photochemwrk_longdydt_get(void *ptr, double *val)
+
 cdef extern void photochemwrk_tn_get(void *ptr, double *val)
 cdef extern void photochemwrk_tn_set(void *ptr, double *val)
 

photochem/fortran/Atmosphere_wrapper.f90

@@ -102,6 +102,26 @@ subroutine atmosphere_photochemical_equilibrium_wrapper(ptr, success, err) bind(
       call copy_string_ftoc(err_f, err)
     endif
   end subroutine
+
+  subroutine atmosphere_check_for_convergence_wrapper(ptr, converged, err) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    logical(c_bool), intent(out) :: converged
+    character(len=c_char), intent(out) :: err(err_len+1)
+
+    character(:), allocatable :: err_f
+    logical :: converged_f
+
+    type(Atmosphere), pointer :: pc
+    call c_f_pointer(ptr, pc)
+
+    converged_f = pc%check_for_convergence(err_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 atmosphere_out2atmosphere_txt_wrapper(ptr, filename, overwrite, clip, err) bind(c)
     type(c_ptr), intent(in) :: ptr

photochem/fortran/EvoAtmosphere_wrapper.f90

@@ -340,6 +340,26 @@ subroutine evoatmosphere_evolve_wrapper(ptr, filename, tstart, nq, nz, usol, nt,
     endif
   end subroutine
 
+  subroutine evoatmosphere_check_for_convergence_wrapper(ptr, converged, err) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    logical(c_bool), intent(out) :: converged
+    character(len=c_char), intent(out) :: err(err_len+1)
+
+    character(:), allocatable :: err_f
+    logical :: converged_f
+
+    type(EvoAtmosphere), pointer :: pc
+    call c_f_pointer(ptr, pc)
+
+    converged_f = pc%check_for_convergence(err_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_initialize_stepper_wrapper(ptr, nq, nz, usol_start, err) bind(c)
     type(c_ptr), intent(in) :: ptr
     integer(c_int), intent(in) :: nq, nz

photochem/fortran/PhotochemVars_wrapper.f90

@@ -321,6 +321,70 @@ subroutine photochemvars_equilibrium_time_set(ptr, val) bind(c)
     var%equilibrium_time = val
   end subroutine
 
+  subroutine photochemvars_conv_hist_factor_get(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%conv_hist_factor
+  end subroutine
+
+  subroutine photochemvars_conv_hist_factor_set(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%conv_hist_factor = val
+  end subroutine
+
+  subroutine photochemvars_conv_min_mix_get(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%conv_min_mix
+  end subroutine
+
+  subroutine photochemvars_conv_min_mix_set(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%conv_min_mix = val
+  end subroutine
+
+  subroutine photochemvars_conv_longdy_get(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%conv_longdy
+  end subroutine
+
+  subroutine photochemvars_conv_longdy_set(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%conv_longdy = val
+  end subroutine
+
+  subroutine photochemvars_conv_longdydt_get(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(out) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%conv_longdydt
+  end subroutine
+
+  subroutine photochemvars_conv_longdydt_set(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(in) :: val
+    type(PhotochemVars), pointer :: var
+    call c_f_pointer(ptr, var)
+    var%conv_longdydt = val
+  end subroutine
+
   subroutine photochemvars_verbose_get(ptr, val) bind(c)
     type(c_ptr), intent(in) :: ptr
     integer(c_int), intent(out) :: val

photochem/fortran/PhotochemWrk_wrapper.f90

@@ -21,6 +21,66 @@ subroutine deallocate_photochemwrk(ptr) bind(c)
   !!! getters and setters !!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+  subroutine photochemwrk_nsteps_get(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(out) :: val
+    type(PhotochemWrk), pointer :: var
+    call c_f_pointer(ptr, var)
+    val = var%nsteps
+  end subroutine
+
+  subroutine photochemwrk_t_history_get_size(ptr, dim1) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(out) :: dim1
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    dim1 = size(wrk%t_history,1)
+  end subroutine
+
+  subroutine photochemwrk_t_history_get(ptr, dim1, arr) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(in) :: dim1
+    real(c_double), intent(out) :: arr(dim1)
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    arr = wrk%t_history
+  end subroutine
+
+  subroutine photochemwrk_mix_history_get_size(ptr, dim1, dim2, dim3) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(out) :: dim1, dim2, dim3
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    dim1 = size(wrk%mix_history,1)
+    dim2 = size(wrk%mix_history,2)
+    dim3 = size(wrk%mix_history,3)
+  end subroutine
+
+  subroutine photochemwrk_mix_history_get(ptr, dim1, dim2, dim3, arr) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(in) :: dim1, dim2, dim3
+    real(c_double), intent(out) :: arr(dim1, dim2, dim3)
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    arr = wrk%mix_history
+  end subroutine
+
+  subroutine photochemwrk_longdy_get(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(out) :: val
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    val = wrk%longdy
+  end subroutine
+
+  subroutine photochemwrk_longdydt_get(ptr, val) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    real(c_double), intent(out) :: val
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    val = wrk%longdydt
+  end subroutine
+
   subroutine photochemwrk_tn_get(ptr, val) bind(c)
     type(c_ptr), intent(in) :: ptr
     real(c_double), intent(out) :: val

src/atmosphere/photochem_atmosphere.f90

@@ -147,7 +147,7 @@ module function evolve(self, filename, tstart, usol_start, t_eval, overwrite, er
       character(:), allocatable, intent(out) :: err
     end function
 
-    !> Determines if integration has converged to photochemical equilibrium or not.
+    !> Determines if integration has converged to photochemical steady-state.
     module function check_for_convergence(self, err) result(converged)
       class(Atmosphere), target, intent(inout) :: self
       character(:), allocatable, intent(out) :: err

src/evoatmosphere/photochem_evoatmosphere.f90

@@ -167,6 +167,7 @@ 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
 
+    !> Determines if integration has converged to photochemical steady-state.
     module function check_for_convergence(self, err) result(converged)
       class(EvoAtmosphere), target, intent(inout) :: self
       character(:), allocatable, intent(out) :: err