convergence is passed upward. lower number of iters

equilibrate/cython/ChemEquiAnalysis_pxd.pxd

@@ -15,7 +15,7 @@ cdef extern void chemequianalysis_create_wrapper(void *ptr, char *thermofile,
 cdef extern void chemequianalysis_solve_wrapper(void *ptr, double *P, double *T,
                          cbool *molfracs_atoms_present, int *molfracs_atoms_dim, double *molfracs_atoms, 
                          cbool *molfracs_species_present, int *molfracs_species_dim, double *molfracs_species, 
-                         char *err)
+                         cbool *converged, char *err)
 
 # Getters and setters
 
@@ -52,5 +52,8 @@ cdef extern void chemequianalysis_molfracs_atoms_condensate_get(void *ptr, int *
 cdef extern void chemequianalysis_molfracs_species_condensate_get_size(void *ptr, int *dim1)
 cdef extern void chemequianalysis_molfracs_species_condensate_get(void *ptr, int *dim1, double *arr)
 
+cdef extern void chemequianalysis_verbose_get(void *ptr, cbool *val)
+cdef extern void chemequianalysis_verbose_set(void *ptr, cbool *val)
+
 cdef extern void chemequianalysis_mass_tol_get(void *ptr, double *val)
 cdef extern void chemequianalysis_mass_tol_set(void *ptr, double *val)

equilibrate/cython/_equilibrate.pyx

@@ -99,16 +99,19 @@ cdef class ChemEquiAnalysis:
       molfracs_species_ = molfracs_species
     cdef int molfracs_species_dim = molfracs_species_.shape[0]
 
+    cdef cbool converged
     cdef char err[ERR_LEN+1]
 
     cea_pxd.chemequianalysis_solve_wrapper(&self._ptr, &P, &T,
                          &molfracs_atoms_present, &molfracs_atoms_dim, <double *>molfracs_atoms_.data, 
                          &molfracs_species_present, &molfracs_species_dim, <double *>molfracs_species_.data, 
-                         err)
+                         &converged, err)
 
     if len(err.strip()) > 0:
       raise EquilibrateException(err.decode("utf-8").strip())
 
+    return converged
+
   property atoms_names:
     "List. Names of atoms"
     def __get__(self):
@@ -208,6 +211,15 @@ cdef class ChemEquiAnalysis:
       cea_pxd.chemequianalysis_molfracs_species_condensate_get(&self._ptr, &dim1, <double *>arr.data)
       return arr
 
+  property verbose:
+    "bool. Determines amount of printing."
+    def __get__(self):
+      cdef cbool val
+      cea_pxd.chemequianalysis_verbose_get(&self._ptr, &val)
+      return val
+    def __set__(self, cbool val):
+      cea_pxd.chemequianalysis_verbose_set(&self._ptr, &val)
+
   property mass_tol:
     """float. Degree to which mass will be balanced. Gordon & McBride's default 
     is 1.0e-6, but it seems like 1.0e-2 is OK.

equilibrate/fortran/equilibrate_c_api.f90

@@ -94,7 +94,7 @@ subroutine chemequianalysis_create_wrapper(ptr, thermofile, &
   subroutine chemequianalysis_solve_wrapper(ptr, P, T, &
                                             molfracs_atoms_present, molfracs_atoms_dim, molfracs_atoms, &
                                             molfracs_species_present, molfracs_species_dim, molfracs_species, &
-                                            err) bind(c)
+                                            converged, err) bind(c)
     use equilibrate, only: ChemEquiAnalysis
     type(c_ptr), intent(in) :: ptr
     real(c_double), intent(in) :: P
@@ -105,6 +105,7 @@ subroutine chemequianalysis_solve_wrapper(ptr, P, T, &
     logical(c_bool), intent(in) :: molfracs_species_present
     integer(c_int), intent(in) :: molfracs_species_dim
     real(c_double), intent(in) :: molfracs_species(molfracs_species_dim)
+    logical(c_bool), intent(out) :: converged
     character(kind=c_char), intent(out) :: err(err_len+1)
 
     character(:), allocatable :: err_f
@@ -113,13 +114,13 @@ subroutine chemequianalysis_solve_wrapper(ptr, P, T, &
     call c_f_pointer(ptr, cea)
 
     if (molfracs_atoms_present .and. molfracs_species_present) then
-      call cea%solve(P, T, molfracs_atoms=molfracs_atoms, molfracs_species=molfracs_species, err=err_f)
+      converged = cea%solve(P, T, molfracs_atoms=molfracs_atoms, molfracs_species=molfracs_species, err=err_f)
     elseif (molfracs_atoms_present .and. .not.molfracs_species_present) then
-      call cea%solve(P, T, molfracs_atoms=molfracs_atoms, err=err_f)
+      converged = cea%solve(P, T, molfracs_atoms=molfracs_atoms, err=err_f)
     elseif (.not.molfracs_atoms_present .and. molfracs_species_present) then
-      call cea%solve(P, T, molfracs_species=molfracs_species, err=err_f)
+      converged = cea%solve(P, T, molfracs_species=molfracs_species, err=err_f)
     elseif (.not.molfracs_atoms_present .and. .not.molfracs_species_present) then
-      call cea%solve(P, T, err=err_f)
+      converged = cea%solve(P, T, err=err_f)
     endif
 
     err(1) = c_null_char
@@ -380,6 +381,24 @@ subroutine chemequianalysis_molfracs_species_condensate_get(ptr, dim1, arr) bind
     arr = cea%molfracs_species_condensate
   end subroutine
 
+  subroutine chemequianalysis_verbose_get(ptr, val) bind(c)
+    use equilibrate, only: ChemEquiAnalysis
+    type(c_ptr), intent(in) :: ptr
+    logical(c_bool), intent(out) :: val
+    type(ChemEquiAnalysis), pointer :: cea
+    call c_f_pointer(ptr, cea)
+    val = cea%verbose
+  end subroutine
+
+  subroutine chemequianalysis_verbose_set(ptr, val) bind(c)
+    use equilibrate, only: ChemEquiAnalysis
+    type(c_ptr), intent(in) :: ptr
+    logical(c_bool), intent(in) :: val
+    type(ChemEquiAnalysis), pointer :: cea
+    call c_f_pointer(ptr, cea)
+    cea%verbose = val
+  end subroutine
+
   subroutine chemequianalysis_mass_tol_get(ptr, val) bind(c)
     use equilibrate, only: ChemEquiAnalysis
     type(c_ptr), intent(in) :: ptr

src/equilibrate.f90

@@ -161,7 +161,7 @@ function create_ChemEquiAnalysis(thermopath, atoms, species, err) result(cea)
   !> If successful, then the equilibrium composition will be stored in a number
   !> of attributes (e.g., self%molfracs_species). If unsuccessful, then the variable
   !> `err` is allocated with an error message.
-  subroutine solve(self, P, T, molfracs_atoms, molfracs_species, err)
+  function solve(self, P, T, molfracs_atoms, molfracs_species, err) result(converged)
     class(ChemEquiAnalysis), intent(inout) :: self
     real(dp), intent(in) :: P !! Pressure in bars
     real(dp), intent(in) :: T !! Temperature in Kelvin
@@ -170,10 +170,13 @@ subroutine solve(self, P, T, molfracs_atoms, molfracs_species, err)
     !> Species mole fractions in the same order and length as self%species_names.
     real(dp), optional, intent(in) :: molfracs_species(:)
     character(:), allocatable, intent(out) :: err
+    logical :: converged
 
     real(dp), allocatable :: molfracs_atoms_(:)
     integer :: i, j, jj
 
+    converged = .false.
+
     if (present(molfracs_atoms) .and. present(molfracs_species)) then
       err = 'Both "molfracs_atoms" and "molfracs_species" are inputs, but only one is allowed.'
       return
@@ -226,7 +229,8 @@ subroutine solve(self, P, T, molfracs_atoms, molfracs_species, err)
     if (self%dat%error) then
       err = trim(self%dat%err_msg)
       return
-    endif 
+    endif
+    converged = self%dat%converged
 
     ! Compute mole fractions of gases and condensates in the solution
     j = 1
@@ -265,6 +269,6 @@ subroutine solve(self, P, T, molfracs_atoms, molfracs_species, err)
     self%molfracs_atoms_condensate = self%molfracs_atoms_condensate/max(sum(self%molfracs_atoms_condensate),tiny(1.0_dp))
     self%molfracs_atoms_gas = self%molfracs_atoms_gas/max(sum(self%molfracs_atoms_gas),tiny(1.0_dp))
 
-  end subroutine
+  end function
 
 end module

src/equilibrate_cea.f90

@@ -12,6 +12,7 @@ module equilibrate_cea
                                        !! Gordon & McBride's default is 1.0e-6, but
                                        !! this seems too strict.
       logical :: verbose2 !! actual verbosity variable
+      logical :: converged !! for passing convergence
 
       ! Run variables
       integer     :: iter_max
@@ -705,6 +706,7 @@ subroutine solve(self,mode,verbo,verbose2,N_atoms_in,N_reactants_in,molfracs_ato
          RETURN
       end if
 
+      self%converged = .false.
       self%verbose = .FALSE.
       self%verbose_cond = (verbo == 'vy')
       self%verbose2 = verbose2
@@ -884,7 +886,7 @@ recursive subroutine ec_COMP_EQU_CHEM(self, N_atoms_use, N_reac, molfracs_atoms,
       slowed = .FALSE.
       call ec_INIT_ALL_VALS(self,N_atoms_use,self%N_reactants,n,n_spec,pi_atom)
 
-      self%iter_max = 50000 + self%N_reactants/2
+      self%iter_max = 50 + self%N_reactants/2
       current_solids_number = 0
 
       MMW = 0d0
@@ -935,10 +937,11 @@ recursive subroutine ec_COMP_EQU_CHEM(self, N_atoms_use, N_reac, molfracs_atoms,
          END IF
       END DO
 
+      self%converged = converged
       IF (.NOT. converged) THEN
-         self%error = .true.
-         self%err_msg = 'One or more convergence criteria not satisfied (gas only).'
-         return
+        if (self%verbose2) then
+          print*,'One or more convergence criteria not satisfied (gas only).'
+        endif
       END IF
 
       converged = .FALSE.
@@ -1057,6 +1060,7 @@ recursive subroutine ec_COMP_EQU_CHEM(self, N_atoms_use, N_reac, molfracs_atoms,
 
                END DO
 
+               self%converged = converged
                IF (.NOT. converged) THEN
                   IF (self%quick) THEN
                      self%quick = .FALSE.
@@ -1072,9 +1076,9 @@ recursive subroutine ec_COMP_EQU_CHEM(self, N_atoms_use, N_reac, molfracs_atoms,
                      slowed = .TRUE.
                      EXIT
                   ELSE
-                     self%error = .true.
-                     self%err_msg = 'One or more convergence criteria not satisfied (gas & condensates).'
-                     return
+                    if (self%verbose2) then
+                       print*,'One or more convergence criteria not satisfied (gas & condensates).'
+                    endif
                   END IF
                END IF
 

test/test_equilibrate.f90

@@ -12,6 +12,7 @@ subroutine test()
     character(s_str_len), allocatable :: species(:)
     character(s_str_len), allocatable :: atoms(:)
     real(dp), allocatable :: X(:), correct_answer(:)
+    logical :: converged
     integer :: i
 
     species = [ &
@@ -262,17 +263,19 @@ subroutine test()
       stop 1
     endif
 
-    call cea%solve(1.0_dp, 1000.0_dp, molfracs_atoms=X, err=err)
+    converged = cea%solve(1.0_dp, 1000.0_dp, molfracs_atoms=X, err=err)
     if (allocated(err)) then
       print*,err
       stop 1
     endif
+    if (.not.converged) stop 1
 
-    call cea2%solve(1.0_dp, 1000.0_dp, molfracs_atoms=X, err=err)
+    converged = cea2%solve(1.0_dp, 1000.0_dp, molfracs_atoms=X, err=err)
     if (allocated(err)) then
       print*,err
       stop 1
     endif
+    if (.not.converged) stop 1
 
     do i = 1,size(cea%molfracs_species)
       if (.not.is_close(cea%molfracs_species(i),correct_answer(i),tol=1.0e-4_dp) .and. cea%molfracs_species(i) > 1.0e-50_dp) then
@@ -290,11 +293,12 @@ subroutine test()
       endif
     enddo
 
-    call cea%solve(1.0_dp, 1000.0_dp, molfracs_species=cea%molfracs_species, err=err)
+    converged = cea%solve(1.0_dp, 1000.0_dp, molfracs_species=cea%molfracs_species, err=err)
     if (allocated(err)) then
       print*,err
       stop 1
     endif
+    if (.not.converged) stop 1
 
     do i = 1,size(cea%molfracs_species)
       if (.not.is_close(cea%molfracs_species(i),correct_answer(i),tol=1.0e-4_dp) .and. cea%molfracs_species(i) > 1.0e-50_dp) then