From 998dec70f1e6bfa5f2b49bcd39188006070bfd21 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Tue, 30 Apr 2024 12:18:19 -0700 Subject: [PATCH] convergence is passed upward. lower number of iters --- equilibrate/cython/ChemEquiAnalysis_pxd.pxd | 5 +++- equilibrate/cython/_equilibrate.pyx | 14 +++++++++- equilibrate/fortran/equilibrate_c_api.f90 | 29 +++++++++++++++++---- src/equilibrate.f90 | 10 ++++--- src/equilibrate_cea.f90 | 18 ++++++++----- test/test_equilibrate.f90 | 10 ++++--- 6 files changed, 66 insertions(+), 20 deletions(-) diff --git a/equilibrate/cython/ChemEquiAnalysis_pxd.pxd b/equilibrate/cython/ChemEquiAnalysis_pxd.pxd index 4eadd9d..7cc53a1 100644 --- a/equilibrate/cython/ChemEquiAnalysis_pxd.pxd +++ b/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) \ No newline at end of file diff --git a/equilibrate/cython/_equilibrate.pyx b/equilibrate/cython/_equilibrate.pyx index 8af9e59..4a1044b 100644 --- a/equilibrate/cython/_equilibrate.pyx +++ b/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, molfracs_atoms_.data, &molfracs_species_present, &molfracs_species_dim, 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, 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. diff --git a/equilibrate/fortran/equilibrate_c_api.f90 b/equilibrate/fortran/equilibrate_c_api.f90 index 2f312c9..88678bc 100644 --- a/equilibrate/fortran/equilibrate_c_api.f90 +++ b/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 diff --git a/src/equilibrate.f90 b/src/equilibrate.f90 index d8597a1..a0ca48b 100644 --- a/src/equilibrate.f90 +++ b/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 \ No newline at end of file diff --git a/src/equilibrate_cea.f90 b/src/equilibrate_cea.f90 index 8dd1c1f..6b58d41 100755 --- a/src/equilibrate_cea.f90 +++ b/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 diff --git a/test/test_equilibrate.f90 b/test/test_equilibrate.f90 index 40a8a08..08b2cfb 100644 --- a/test/test_equilibrate.f90 +++ b/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