From ef69d541c8c89587f7b4778a19eda14b922866d7 Mon Sep 17 00:00:00 2001 From: rs028 Date: Fri, 19 Dec 2025 11:49:50 +0000 Subject: [PATCH 01/11] Upgrade openlibm and cvode --- tools/install/install_cvode.sh | 3 ++- tools/install/install_openlibm.sh | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/tools/install/install_cvode.sh b/tools/install/install_cvode.sh index 3a9cec86..090ca546 100755 --- a/tools/install/install_cvode.sh +++ b/tools/install/install_cvode.sh @@ -26,7 +26,7 @@ # ./install_cvode.sh ~/path/to/dependencies/directory /path/to/fortran/compiler # ----------------------------------------------------------------------------- -SUNDIALS_VERSION="2.7.0" +SUNDIALS_VERSION="7.5.0" # path to dependencies directory if [ -z "$1" ] ; then @@ -89,6 +89,7 @@ cmake -DCMAKE_INSTALL_PREFIX="$DEP_DIR/cvode" \ -DFCMIX_ENABLE:BOOL=ON \ -DEXAMPLES_ENABLE:BOOL=OFF \ -DCMAKE_MACOSX_RPATH:BOOL=ON \ + -DBUILD_FORTRAN_MODULE_INTERFACE:BOOL=ON \ .. if [ $? -ne 0 ] ; then printf "\n[cvode] cmake --> FAIL\n" diff --git a/tools/install/install_openlibm.sh b/tools/install/install_openlibm.sh index af4f7332..5405a56d 100755 --- a/tools/install/install_openlibm.sh +++ b/tools/install/install_openlibm.sh @@ -21,7 +21,7 @@ # ./install_openlibm.sh ~/path/to/dependencies/directory # ----------------------------------------------------------------------------- -OPENLIBM_VERSION="0.8.6" +OPENLIBM_VERSION="0.8.7" # path to dependencies directory if [ -z "$1" ] ; then From 29aec00d568890750e0e3352fc1e5ebae9d6a681 Mon Sep 17 00:00:00 2001 From: rs028 Date: Sat, 27 Dec 2025 15:26:23 +0000 Subject: [PATCH 02/11] Adapt to use sundials 7.5 --- src/atchem2.f90 | 216 ++++++++++++++++++++++++++++------------ src/outputFunctions.f90 | 121 +++++++++++++++------- 2 files changed, 233 insertions(+), 104 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 49692166..3b32d490 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -1,14 +1,14 @@ ! ----------------------------------------------------------------------------- ! -! Copyright (c) 2009 - 2012 Chris Martin, Kasia Boronska, Jenny Young, +! Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young, ! Peter Jimack, Mike Pilling ! -! Copyright (c) 2017 - 2018 Sam Cox, Roberto Sommariva +! Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, Neil Butcher ! ! This file is part of the AtChem2 software package. ! -! This file is covered by the MIT license which can be found in the file -! LICENSE.md at the top level of the AtChem2 distribution. +! This file is licensed under the MIT license, which can be found in the file +! `LICENSE` at the top level of the AtChem2 distribution. ! ! ----------------------------------------------------------------------------- @@ -18,6 +18,17 @@ ! ! ******************************************************************** ! +module cvode_rhs_mod + use, intrinsic :: iso_c_binding + implicit none + + type, bind(C) :: UserData + integer(c_int) :: ipar(10) + real(c_double) :: rpar(1) + end type UserData + +end module cvode_rhs_mod + PROGRAM ATCHEM2 use, intrinsic :: iso_fortran_env, only : stderr => error_unit @@ -41,6 +52,17 @@ PROGRAM ATCHEM2 use output_functions_mod use constraint_functions_mod, only : addConstrainedSpeciesToProbSpec, removeConstrainedSpeciesFromProbSpec use solver_functions_mod, only : jfy, proc + use cvode_rhs_mod, only : UserData + + ! Sundials module + use fsundials_core_mod + use fnvector_serial_mod + use fcvode_mod + use fsunlinsol_spgmr_mod + use fsunlinsol_dense_mod + use fsunmatrix_dense_mod + use fsunmatrix_band_mod + use fsunlinsol_band_mod implicit none ! interface to linux API @@ -77,10 +99,10 @@ function dlclose( handle ) bind ( c, name="dlclose" ) ! Declarations for solver parameters integer(kind=QI) :: ier integer :: meth, itmeth, iatol, itask, currentNumTimestep - integer(kind=NPI) :: iout(21), ipar(10) + integer(kind=NPI) :: ipar(10) integer(kind=NPI) :: neq real(kind=DP) :: t, tout - real(kind=DP) :: rout(6), rpar(1) + real(kind=DP) :: rpar(1) ! Walltime variables integer(kind=QI) :: runStart, runEnd, runTime, clockRate @@ -110,6 +132,17 @@ function dlclose( handle ) bind ( c, name="dlclose" ) integer(c_int), parameter :: rtld_lazy=1 ! value extracted from the C header file integer(c_int), parameter :: rtld_now=2 ! value extracted from the C header file + !sundials declarations + + type(c_ptr) :: ctx ! SUNDIALS context for the simulation + type(N_Vector), pointer :: sunvec_u ! sundials vector + type(SUNLinearSolver), pointer :: sunls ! sundials linear solver + type(SUNMatrix), pointer :: sunmat_A ! sundials matrix (empty) + type(c_ptr) :: cvode_mem ! CVODE memory + real(c_double) :: t_arr(1) + + type(UserData), target :: udata + ! ***************************************************************** ! Explicit declaration of FCVFUN() interface, which is a ! user-supplied function to CVODE. @@ -136,6 +169,19 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) integer(kind=NPI) :: i end subroutine FCVFUN + integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) + use, intrinsic :: iso_c_binding + use fsundials_core_mod + use fnvector_serial_mod + + implicit none + + real(c_double), value, intent(in) :: t + type(N_Vector), intent(inout) :: y + type(N_Vector), intent(inout) :: ydot + type(c_ptr), value, intent(in) :: user_data + end function rhs_fn + end interface ! ***************************************************************** @@ -145,9 +191,7 @@ end subroutine FCVFUN call SYSTEM_CLOCK( runStart ) ! Initialise some variables used by CVODE functions to invalid values - iout(:) = -1_NPI ipar(:) = -1_NPI - rout(:) = -1.0_DP rpar(:) = -1.0_DP write (*, '(A)') 'AtChem2 v1.3-dev' @@ -269,10 +313,8 @@ end subroutine FCVFUN t = modelStartTime call calcCurrentDateParameters( t ) tout = timestepSize + t - ! Parameters for FCVMALLOC(). (Comments from cvode guide) meth - ! specifies the basic integration: 1 for Adams (nonstiff) or 2 for - ! BDF stiff) - meth = 2 + ! Parameters for FCVodeCreate(). Adams (nonstiff) or BDF (stiff) + meth = CV_BDF ! itmeth specifies the nonlinear iteration method: 1 for functional ! iteration or 2 for Newton iteration. itmeth = 2 @@ -352,66 +394,82 @@ end subroutine FCVFUN ! CONFIGURE SOLVER ! ***************************************************************** + ! create the SUNDIALS context + ier = FSUNContext_Create(SUN_COMM_NULL, ctx) ipar(1) = neq ipar(2) = numReac - call FNVINITS( 1, neq, ier ) - if ( ier /= 0 ) then - write (stderr, 20) ier - 20 format (///' SUNDIALS_ERROR: FNVINITS() returned ier = ', I5) - stop + ! create SUNDIALS N_Vector + sunvec_u => FN_VMake_Serial(neq, z, ctx) + if (.not. associated(sunvec_u)) then + print *, 'ERROR: sunvec = NULL' + stop 1 end if + + write (*, '(A30, 1P e15.3) ') ' t0 = ', t write (*,*) - call FCVMALLOC( t, z, meth, itmeth, iatol, rtol, atol, & - iout, rout, ipar, rpar, ier ) - if ( ier /= 0 ) then - write (stderr, 30) ier - 30 format (///' SUNDIALS_ERROR: FCVMALLOC() returned ier = ', I5) - stop + + ! create and initialize CVode memory + cvode_mem = FCVodeCreate(meth, ctx) + if (.not. c_associated(cvode_mem)) print *, 'ERROR: cvode_mem = NULL' + + ier = FCVodeInit(cvode_mem, c_funloc(rhs_fn), t, sunvec_u) + if (ier /= 0) then + print *, 'Error in FCVodeInit, ierr = ', ier, '; halting' + stop 1 + end if + + ier = FCVodeSStolerances(cvode_mem, rtol, atol) + if (ier /= 0) then + print *, 'Error in FCVodeSStolerances, ierr = ', ier, '; halting' + stop 1 + end if + + ier = FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) + if (ier /= 0) then + print *, 'Error in FCVodeSetMaxNumSteps, ierr = ', ier, '; halting' + stop 1 end if - call FCVSETIIN( 'MAX_NSTEPS', maxNumInternalSteps, ier ) - write (*, '(A, I0)') ' setting maxnumsteps ier = ', ier - - call FCVSETRIN( 'MAX_STEP', maxStep, ier ) + ier = FCVodeSetMaxStep(cvode_mem, real(maxStep, kind=C_DOUBLE)) write (*, '(A, I0)') ' setting maxstep ier = ', ier write (*,*) + udata%ipar = ipar + udata%rpar = rpar + + ier = FCVodeSetUserData(cvode_mem, c_loc(udata)) + if (ier /= 0) stop 'User data setup failed' + ! SELECT SOLVER TYPE ACCORDING TO FILE INPUT ! SPGMR SOLVER if ( solverType == 1 ) then - call FCVSPGMR( 0, 1, lookBack, deltaMain, ier ) - ! SPGMR SOLVER WITH BANDED PRECONDITIONER + sunls => FSUNLinSol_SPGMR(sunvec_u, SUN_PREC_NONE, int(lookBack, kind=C_INT), ctx) + ! SPGMR SOLVER WITH BANDED PRECONDITIONER else if ( solverType == 2 ) then - call FCVSPGMR( 1, 1, lookBack, deltaMain, ier ) - call FCVBPINIT( neq, preconBandUpper, preconBandLower, ier ) - if ( ier /= 0 ) then - write (stderr,*) 'SUNDIALS_ERROR: preconditioner returned ier = ', ier ; - call FCVFREE() - stop - end if - ! DENSE SOLVER + sunmat_A => FSUNBandMatrix(neq, preconBandUpper, preconBandLower, ctx) + sunls => FSUNLinSol_Band(sunvec_u, sunmat_A, ctx) + ! DENSE SOLVER else if ( solverType == 3 ) then - call FCVDENSE( neq, ier ) - ! UNEXPECTED SOLVER TYPE + ! Create dense SUNMatrix for use in linear solves + sunmat_A => FSUNDenseMatrix(neq, neq, ctx) + sunls => FSUNLinSol_Dense(sunvec_u, sunmat_A, ctx) + ! UNEXPECTED SOLVER TYPE else write (stderr,*) 'Error with solverType input, input = ', solverType write (stderr,*) 'Available options are 1, 2, 3.' stop end if - ! ERROR HANDLING - if ( ier /= 0 ) then - write (stderr,*) ' SUNDIALS_ERROR: SOLVER returned ier = ', ier - call FCVFREE() - stop - end if + + ! Attach the matrix and linear solver + ier = FCVodeSetLinearSolver(cvode_mem, sunls, sunmat_A); + ! ERROR HANDLING if ( ier /= 0 ) then - write (stderr, 40) ier - 40 format (///' SUNDIALS_ERROR: FCVDENSE() returned ier = ', I5) - call FCVFREE() + write (stderr,*) ' SUNDIALS_ERROR: FCVodeSetLinearSolver returned ier = ', ier + call FCVodeFree(cvode_mem) stop end if @@ -440,10 +498,12 @@ end subroutine FCVFUN end if ! Get concentrations for unconstrained species - call FCVODE( tout, t, z, itask, ier ) + t_arr(1) =t + ier = FCVode(cvode_mem, tout, sunvec_u, t_arr, itask) if ( ier /= 0 ) then write (*, '(A, I0)') ' ier POST FCVODE()= ', ier end if + t = t_arr(1) flush(6) time = nint( t ) @@ -470,9 +530,6 @@ end subroutine FCVFUN call outputreactionRates( time ) end if - ! Output CVODE solver parameters and timestep sizes - call outputSolverParameters( t, rout(3), rout(2), iout, solverType ) - ! Output envVar values ro2 = ro2sum( speciesConcs ) call outputEnvVar( t ) @@ -480,9 +537,8 @@ end subroutine FCVFUN ! Error handling if ( ier < 0 ) then fmt = "(///' SUNDIALS_ERROR: FCVODE() returned ier = ', I5, /, 'Linear Solver returned ier = ', I5) " - write (stderr, fmt) ier, iout (15) ! free memory - call FCVFREE() + call FCVodeFree(cvode_mem) stop end if @@ -495,21 +551,13 @@ end subroutine FCVFUN ! Output final model concentrations, in a usable format for model ! restart call outputFinalModelState( getSpeciesList(), speciesConcs ) - write (*,*) - + write (*, '(A)') '------------------' write (*, '(A)') ' Final statistics' write (*, '(A)') '------------------' + call PrintFinalStats(cvode_mem) + write (*,*) - ! Final on-screen output - fmt = "(' No. steps = ', I0, ' No. f-s = ', I0, " // & - "' No. J-s = ', I0, ' No. LU-s = ', I0/" // & - "' No. nonlinear iterations = ', I0/" // & - "' No. nonlinear convergence failures = ', I0/" // & - "' No. error test failures = ', I0/) " - - write (*, fmt) iout (3), iout (4), iout (17), iout (8), & - iout (7), iout (6), iout (5) call SYSTEM_CLOCK( runEnd, clockRate ) runTime = ( runEnd - runStart ) / clockRate @@ -521,7 +569,7 @@ end subroutine FCVFUN ! ***************************************************************** ! deallocate CVODE internal data - call FCVFREE() + call FCVodeFree(cvode_mem) deallocate (speciesConcs, z) deallocate (reacDetailedRatesSpecies, prodDetailedRatesSpecies) deallocate (detailedRatesSpeciesName, speciesOfInterest) @@ -626,3 +674,39 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) end subroutine FCVFUN ! ******************************************************************** ! +integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) + use, intrinsic :: iso_c_binding + use types_mod + use fsundials_core_mod + use fnvector_serial_mod + use cvode_rhs_mod, only : UserData + use types_mod + + implicit none + + real(c_double), value, intent(in) :: t + type(N_Vector), intent(inout) :: y + type(N_Vector), intent(inout) :: ydot + type(c_ptr), value, intent(in) :: user_data + + ! Local pointers to vector data + real(c_double), pointer :: ydata(:) + real(c_double), pointer :: ydotdata(:) + + type(UserData), pointer :: ud + integer(kind=NPI) :: ier + + integer(kind=NPI) :: ipar_f(10) + integer :: i + call c_f_pointer(user_data, ud) + ipar_f = [(int(ud%ipar(i), kind=NPI), i=1, size(ud%ipar))] + + ! Get raw data arrays from N_Vector + ydata => FN_VGetArrayPointer(y) + ydotdata => FN_VGetArrayPointer(ydot) + + + call FCVFUN( t, ydata, ydotdata, ipar_f, ud%rpar, ier ) + + rhs_fn = ier ! 0 = success +end function rhs_fn diff --git a/src/outputFunctions.f90 b/src/outputFunctions.f90 index b9693086..e809b089 100644 --- a/src/outputFunctions.f90 +++ b/src/outputFunctions.f90 @@ -1,14 +1,14 @@ ! ----------------------------------------------------------------------------- ! -! Copyright (c) 2009 - 2012 Chris Martin, Kasia Boronska, Jenny Young, +! Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young, ! Peter Jimack, Mike Pilling ! -! Copyright (c) 2017 Sam Cox, Roberto Sommariva +! Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, Neil Butcher ! ! This file is part of the AtChem2 software package. ! -! This file is covered by the MIT license which can be found in the file -! LICENSE.md at the top level of the AtChem2 distribution. +! This file is licensed under the MIT license, which can be found in the file +! `LICENSE` at the top level of the AtChem2 distribution. ! ! ----------------------------------------------------------------------------- @@ -63,46 +63,91 @@ subroutine outputEnvVar( t ) return end subroutine outputEnvVar - ! ----------------------------------------------------------------- - ! Write parameters output by CVODE solver to file. - subroutine outputSolverParameters( t, prev, this, array, solver_type ) - use, intrinsic :: iso_fortran_env, only : stderr => error_unit - use types_mod + subroutine PrintFinalStats(cvode_mem) - real(kind=DP), intent(in) :: t, prev, this - integer(kind=NPI), intent(in) :: array(:) - integer(kind=SI), intent(in) :: solver_type - integer(kind=SI) :: i - logical :: first_time = .true. + !======= Inclusions =========== + use iso_c_binding + use fcvode_mod - if ( ( solver_type == 1 ) .or. ( solver_type == 2 ) ) then - ! CVSPILS type solver - if ( first_time .eqv. .true. ) then - write (57, '(A9, 2A17, 20A9) ') 't', 'currentStepSize', 'previousStepSize', 'LENRW', 'LENIW', 'NST', 'NFE', & - 'NETF', 'NCFN', 'NNI', 'NSETUPS', 'QU', 'QCUR', 'NOR', 'LENRWLS', 'LENIWLS', & - 'LS_FLAG', 'NFELS', 'NJTV', 'NPE', 'NPS', 'NLI', 'NCFL' - first_time = .false. - end if - write (57, '(1P e9.2, 2 (ES17.8E3), 20I9) ') t, prev, this, (array(i), i = 1, 11), (array(i), i = 13, 21) - - else if ( solver_type == 3 ) then - ! CVDLS type solver - if ( first_time .eqv. .true. ) then - write (57, '(A9, 2A17, 16A9) ') 't', 'currentStepSize', 'previousStepSize', 'LENRW', 'LENIW', 'NST', 'NFE', & - 'NETF', 'NCFN', 'NNI', 'NSETUPS', 'QU', 'QCUR', 'NOR', 'LENRWLS', 'LENIWLS', & - 'LS_FLAG', 'NFELS', 'NJE' - first_time = .false. - end if - write (57, '(1P e9.2, 2 (ES17.8E3), 16I9) ') t, prev, this, (array(i), i = 1, 11), (array(i), i = 13, 17) + !======= Declarations ========= + implicit none + + type(c_ptr), intent(in) :: cvode_mem ! solver memory structure - else - write (stderr,*) 'outputSolverParameters(): Error with solver_type = ', solver_type - write (stderr,*) 'Available options are 1, 2, 3.' - stop + integer(c_int) :: retval ! error flag + + integer(c_long) :: nsteps(1) ! num steps + integer(c_long) :: nfe(1) ! num function evals + integer(c_long) :: netfails(1) ! num error test fails + integer(c_long) :: nniters(1) ! nonlinear solver iterations + integer(c_long) :: nncfails(1) ! nonlinear solver fails + integer(c_long) :: njacevals(1) ! number of Jacobian evaluations + integer(c_long) :: nluevals(1) ! number of LU evals + integer(c_long) :: ngevals(1) ! number of root evals + + !======= Internals ============ + + retval = FCVodeGetNumSteps(cvode_mem, nsteps) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumSteps, retval = ', retval, '; halting' + stop 1 end if + retval = FCVodeGetNumRhsEvals(cvode_mem, nfe) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumRhsEvals, retval = ', retval, '; halting' + stop 1 + end if + + retval = FCVodeGetNumLinSolvSetups(cvode_mem, nluevals) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumLinSolvSetups, retval = ', retval, '; halting' + stop 1 + end if + + retval = FCVodeGetNumErrTestFails(cvode_mem, netfails) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumErrTestFails, retval = ', retval, '; halting' + stop 1 + end if + + retval = FCVodeGetNumNonlinSolvIters(cvode_mem, nniters) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumNonlinSolvIters, retval = ', retval, '; halting' + stop 1 + end if + + retval = FCVodeGetNumNonlinSolvConvFails(cvode_mem, nncfails) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumNonlinSolvConvFails, retval = ', retval, '; halting' + stop 1 + end if + + retval = FCVodeGetNumJacEvals(cvode_mem, njacevals) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumJacEvals, retval = ', retval, '; halting' + stop 1 + end if + + retval = FCVodeGetNumGEvals(cvode_mem, ngevals) + if (retval /= 0) then + print *, 'Error in FCVodeGetNumGEvals, retval = ', retval, '; halting' + stop 1 + end if + + print '(4x,A,i9)', 'Total internal steps taken =', nsteps + print '(4x,A,i9)', 'Total rhs function calls =', nfe + print '(4x,A,i9)', 'Total Jacobian function calls =', njacevals + print '(4x,A,i9)', 'Total root function calls =', ngevals + print '(4x,A,i9)', 'Total LU function calls =', nluevals + print '(4x,A,i9)', 'Num error test failures =', netfails + print '(4x,A,i9)', 'Num nonlinear solver iters =', nniters + print '(4x,A,i9)', 'Num nonlinear solver fails =', nncfails + print *, ' ' + return - end subroutine outputSolverParameters + + end subroutine PrintFinalStats ! ----------------------------------------------------------------- ! Write parameters used in calculation of photolysis rates to file. From 44fc3c02b3bc23b120fd7c259ef22a0100a77fa2 Mon Sep 17 00:00:00 2001 From: rs028 Date: Sat, 27 Dec 2025 16:50:56 +0000 Subject: [PATCH 03/11] Update credits and copyright --- doc/latex/Credits.tex | 7 ++++--- tests/makefile.tests | 8 ++++---- tools/install/Makefile.skel | 3 ++- tools/install/install_cvode.sh | 6 +++--- tools/install/install_openlibm.sh | 6 +++--- 5 files changed, 16 insertions(+), 14 deletions(-) diff --git a/doc/latex/Credits.tex b/doc/latex/Credits.tex index 5f96a989..3df31e10 100644 --- a/doc/latex/Credits.tex +++ b/doc/latex/Credits.tex @@ -1,11 +1,11 @@ % ----------------------------------------------------------------------------- % -% Copyright (c) 2017 Sam Cox, Roberto Sommariva +% Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva % % This file is part of the AtChem2 software package. % -% This file is covered by the MIT license which can be found in the file -% LICENSE.md at the top level of the AtChem2 distribution. +% This file is licensed under the MIT license, which can be found in the file +% `LICENSE` at the top level of the AtChem2 distribution. % % ----------------------------------------------------------------------------- @@ -27,6 +27,7 @@ \section{Credits} \label{sec:credits} \begin{itemize} \item James Allsopp +\item Neil Butcher \item Will Drysdale \item Maarten Fabr{\'e} \item Alfred Mayhew diff --git a/tests/makefile.tests b/tests/makefile.tests index 088f0b8e..7db9e066 100644 --- a/tests/makefile.tests +++ b/tests/makefile.tests @@ -1,14 +1,14 @@ # ----------------------------------------------------------------------------- # -# Copyright (c) 2009 - 2012 Chris Martin, Kasia Boronska, Jenny Young, +# Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young, # Peter Jimack, Mike Pilling # -# Copyright (c) 2017 Sam Cox, Roberto Sommariva +# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, James Allsopp # # This file is part of the AtChem2 software package. # -# This file is covered by the MIT license which can be found in the file -# LICENSE.md at the top level of the AtChem2 distribution. +# This file is licensed under the MIT license, which can be found in the file +# `LICENSE` at the top level of the AtChem2 distribution. # # ----------------------------------------------------------------------------- diff --git a/tools/install/Makefile.skel b/tools/install/Makefile.skel index 0a145eac..53d83687 100644 --- a/tools/install/Makefile.skel +++ b/tools/install/Makefile.skel @@ -3,7 +3,8 @@ # Copyright (c) 2009-2012 Chris Martin, Kasia Boronska, Jenny Young, # Peter Jimack, Mike Pilling # -# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva +# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, James Allsopp, +# Neil Butcher # # This file is part of the AtChem2 software package. # diff --git a/tools/install/install_cvode.sh b/tools/install/install_cvode.sh index 090ca546..777cbd13 100755 --- a/tools/install/install_cvode.sh +++ b/tools/install/install_cvode.sh @@ -1,12 +1,12 @@ #!/bin/sh # ----------------------------------------------------------------------------- # -# Copyright (c) 2017 Sam Cox, Roberto Sommariva +# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva, Neil Butcher # # This file is part of the AtChem2 software package. # -# This file is covered by the MIT license which can be found in the file -# LICENSE.md at the top level of the AtChem2 distribution. +# This file is licensed under the MIT license, which can be found in the file +# `LICENSE` at the top level of the AtChem2 distribution. # # ----------------------------------------------------------------------------- diff --git a/tools/install/install_openlibm.sh b/tools/install/install_openlibm.sh index 5405a56d..b75bce54 100755 --- a/tools/install/install_openlibm.sh +++ b/tools/install/install_openlibm.sh @@ -1,12 +1,12 @@ #!/bin/sh # ----------------------------------------------------------------------------- # -# Copyright (c) 2017 Sam Cox, Roberto Sommariva +# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva # # This file is part of the AtChem2 software package. # -# This file is covered by the MIT license which can be found in the file -# LICENSE.md at the top level of the AtChem2 distribution. +# This file is licensed under the MIT license, which can be found in the file +# `LICENSE` at the top level of the AtChem2 distribution. # # ----------------------------------------------------------------------------- From ace0e3594a8c263b3e62044514d159a793d9deef Mon Sep 17 00:00:00 2001 From: rs028 Date: Mon, 29 Dec 2025 11:16:49 +0000 Subject: [PATCH 04/11] Modify Makefile for sundials upgrade --- tools/install/Makefile.skel | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/tools/install/Makefile.skel b/tools/install/Makefile.skel index 53d83687..5a02145a 100644 --- a/tools/install/Makefile.skel +++ b/tools/install/Makefile.skel @@ -23,6 +23,7 @@ FORTC = "gnu" # Set the dependencies paths. Use full paths, not relative paths. # For example: $(HOME)/path/to/dependencies/directory/cvode/lib CVODELIBDIR = cvode/lib +CVODEOBJDIR = cvode/fortran OPENLIBMDIR = openlibm FRUITDIR = fruit_3.4.3 @@ -84,8 +85,17 @@ endif # set the CVODE and openlibm compilation flags LDFLAGS = -L$(CVODELIBDIR) -L$(OPENLIBMDIR) \ -Wl,$(RPATH_OPTION),/usr/lib/:$(CVODELIBDIR):$(OPENLIBMDIR) \ - -lopenlibm -lsundials_fcvode -lsundials_cvode -lsundials_fnvecserial \ - -lsundials_nvecserial -ldl + -lopenlibm \ + -lsundials_fcvode_mod \ + -lsundials_fnvecserial_mod \ + -lsundials_fcore_mod \ + -lsundials_fsunlinsolspgmr_mod \ + -lsundials_fsunmatrixband_mod \ + -lsundials_fsunlinsoldense_mod \ + -lsundials_fsunmatrixdense_mod \ + -lsundials_cvode \ + -lsundials_nvecserial \ + -ldl # ---------------------------------------------------- # EXECUTABLE setup @@ -109,7 +119,7 @@ SRCS = $(CORE_SRCS) $(SRC)/atchem2.f90 # the executable is rebuilt every time one of the fortran files (in $SRCS) # is modified $(AOUT): $(SRCS) - $(FORT_COMP) -o $(AOUT) -J$(OBJ) -I$(OBJ) $(SRCS) $(FFLAGS) $(LDFLAGS) + $(FORT_COMP) -o $(AOUT) -J$(OBJ) -I$(OBJ) -I$(CVODEOBJDIR) $(SRCS) $(FFLAGS) $(LDFLAGS) # secondary makefile for the Testsuite include tests/makefile.tests From 564cf754a4c33ef41083d664b2b461267313990e Mon Sep 17 00:00:00 2001 From: rs028 Date: Wed, 7 Jan 2026 14:58:02 +0000 Subject: [PATCH 05/11] Fix fortraan style --- src/atchem2.f90 | 57 ++++++++++++++++++++--------------------- src/outputFunctions.f90 | 38 +++++++++++++-------------- 2 files changed, 47 insertions(+), 48 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 3b32d490..2d390831 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -23,8 +23,8 @@ module cvode_rhs_mod implicit none type, bind(C) :: UserData - integer(c_int) :: ipar(10) - real(c_double) :: rpar(1) + integer(c_int) :: ipar(10) + real(c_double) :: rpar(1) end type UserData end module cvode_rhs_mod @@ -54,7 +54,7 @@ PROGRAM ATCHEM2 use solver_functions_mod, only : jfy, proc use cvode_rhs_mod, only : UserData - ! Sundials module + ! SUNDIALS module use fsundials_core_mod use fnvector_serial_mod use fcvode_mod @@ -132,13 +132,12 @@ function dlclose( handle ) bind ( c, name="dlclose" ) integer(c_int), parameter :: rtld_lazy=1 ! value extracted from the C header file integer(c_int), parameter :: rtld_now=2 ! value extracted from the C header file - !sundials declarations - - type(c_ptr) :: ctx ! SUNDIALS context for the simulation - type(N_Vector), pointer :: sunvec_u ! sundials vector - type(SUNLinearSolver), pointer :: sunls ! sundials linear solver - type(SUNMatrix), pointer :: sunmat_A ! sundials matrix (empty) - type(c_ptr) :: cvode_mem ! CVODE memory + ! SUNDIALS declarations + type(c_ptr) :: ctx ! SUNDIALS context for the simulation + type(N_Vector), pointer :: sunvec_u ! sundials vector + type(SUNLinearSolver), pointer :: sunls ! sundials linear solver + type(SUNMatrix), pointer :: sunmat_A ! sundials matrix (empty) + type(c_ptr) :: cvode_mem ! CVODE memory real(c_double) :: t_arr(1) type(UserData), target :: udata @@ -177,9 +176,9 @@ integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) implicit none real(c_double), value, intent(in) :: t - type(N_Vector), intent(inout) :: y - type(N_Vector), intent(inout) :: ydot - type(c_ptr), value, intent(in) :: user_data + type(N_Vector), intent(inout) :: y + type(N_Vector), intent(inout) :: ydot + type(c_ptr), value, intent(in) :: user_data end function rhs_fn end interface @@ -410,23 +409,23 @@ end function rhs_fn write (*, '(A30, 1P e15.3) ') ' t0 = ', t write (*,*) - + ! create and initialize CVode memory cvode_mem = FCVodeCreate(meth, ctx) if (.not. c_associated(cvode_mem)) print *, 'ERROR: cvode_mem = NULL' - + ier = FCVodeInit(cvode_mem, c_funloc(rhs_fn), t, sunvec_u) if (ier /= 0) then print *, 'Error in FCVodeInit, ierr = ', ier, '; halting' stop 1 end if - + ier = FCVodeSStolerances(cvode_mem, rtol, atol) if (ier /= 0) then print *, 'Error in FCVodeSStolerances, ierr = ', ier, '; halting' stop 1 end if - + ier = FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) if (ier /= 0) then print *, 'Error in FCVodeSetMaxNumSteps, ierr = ', ier, '; halting' @@ -463,13 +462,13 @@ end function rhs_fn stop end if - + ! Attach the matrix and linear solver - ier = FCVodeSetLinearSolver(cvode_mem, sunls, sunmat_A); + ier = FCVodeSetLinearSolver(cvode_mem, sunls, sunmat_A); ! ERROR HANDLING if ( ier /= 0 ) then write (stderr,*) ' SUNDIALS_ERROR: FCVodeSetLinearSolver returned ier = ', ier - call FCVodeFree(cvode_mem) + call FCVodeFree( cvode_mem ) stop end if @@ -538,7 +537,7 @@ end function rhs_fn if ( ier < 0 ) then fmt = "(///' SUNDIALS_ERROR: FCVODE() returned ier = ', I5, /, 'Linear Solver returned ier = ', I5) " ! free memory - call FCVodeFree(cvode_mem) + call FCVodeFree( cvode_mem ) stop end if @@ -551,11 +550,11 @@ end function rhs_fn ! Output final model concentrations, in a usable format for model ! restart call outputFinalModelState( getSpeciesList(), speciesConcs ) - + write (*, '(A)') '------------------' write (*, '(A)') ' Final statistics' write (*, '(A)') '------------------' - call PrintFinalStats(cvode_mem) + call PrintFinalStats( cvode_mem ) write (*,*) @@ -569,7 +568,7 @@ end function rhs_fn ! ***************************************************************** ! deallocate CVODE internal data - call FCVodeFree(cvode_mem) + call FCVodeFree( cvode_mem ) deallocate (speciesConcs, z) deallocate (reacDetailedRatesSpecies, prodDetailedRatesSpecies) deallocate (detailedRatesSpeciesName, speciesOfInterest) @@ -685,9 +684,9 @@ integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) implicit none real(c_double), value, intent(in) :: t - type(N_Vector), intent(inout) :: y - type(N_Vector), intent(inout) :: ydot - type(c_ptr), value, intent(in) :: user_data + type(N_Vector), intent(inout) :: y + type(N_Vector), intent(inout) :: ydot + type(c_ptr), value, intent(in) :: user_data ! Local pointers to vector data real(c_double), pointer :: ydata(:) @@ -697,8 +696,8 @@ integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) integer(kind=NPI) :: ier integer(kind=NPI) :: ipar_f(10) - integer :: i - call c_f_pointer(user_data, ud) + integer :: i + call c_f_pointer( user_data, ud ) ipar_f = [(int(ud%ipar(i), kind=NPI), i=1, size(ud%ipar))] ! Get raw data arrays from N_Vector diff --git a/src/outputFunctions.f90 b/src/outputFunctions.f90 index e809b089..ca5d5b29 100644 --- a/src/outputFunctions.f90 +++ b/src/outputFunctions.f90 @@ -63,7 +63,7 @@ subroutine outputEnvVar( t ) return end subroutine outputEnvVar - subroutine PrintFinalStats(cvode_mem) + subroutine PrintFinalStats( cvode_mem ) !======= Inclusions =========== use iso_c_binding @@ -74,16 +74,16 @@ subroutine PrintFinalStats(cvode_mem) type(c_ptr), intent(in) :: cvode_mem ! solver memory structure - integer(c_int) :: retval ! error flag + integer(c_int) :: retval ! error flag - integer(c_long) :: nsteps(1) ! num steps - integer(c_long) :: nfe(1) ! num function evals - integer(c_long) :: netfails(1) ! num error test fails - integer(c_long) :: nniters(1) ! nonlinear solver iterations - integer(c_long) :: nncfails(1) ! nonlinear solver fails - integer(c_long) :: njacevals(1) ! number of Jacobian evaluations - integer(c_long) :: nluevals(1) ! number of LU evals - integer(c_long) :: ngevals(1) ! number of root evals + integer(c_long) :: nsteps(1) ! num steps + integer(c_long) :: nfe(1) ! num function evals + integer(c_long) :: netfails(1) ! num error test fails + integer(c_long) :: nniters(1) ! nonlinear solver iterations + integer(c_long) :: nncfails(1) ! nonlinear solver fails + integer(c_long) :: njacevals(1) ! number of Jacobian evaluations + integer(c_long) :: nluevals(1) ! number of LU evals + integer(c_long) :: ngevals(1) ! number of root evals !======= Internals ============ @@ -134,15 +134,15 @@ subroutine PrintFinalStats(cvode_mem) print *, 'Error in FCVodeGetNumGEvals, retval = ', retval, '; halting' stop 1 end if - - print '(4x,A,i9)', 'Total internal steps taken =', nsteps - print '(4x,A,i9)', 'Total rhs function calls =', nfe - print '(4x,A,i9)', 'Total Jacobian function calls =', njacevals - print '(4x,A,i9)', 'Total root function calls =', ngevals - print '(4x,A,i9)', 'Total LU function calls =', nluevals - print '(4x,A,i9)', 'Num error test failures =', netfails - print '(4x,A,i9)', 'Num nonlinear solver iters =', nniters - print '(4x,A,i9)', 'Num nonlinear solver fails =', nncfails + + print '(4x, A, i9)', 'Total internal steps taken =', nsteps + print '(4x, A, i9)', 'Total rhs function calls =', nfe + print '(4x, A, i9)', 'Total Jacobian function calls =', njacevals + print '(4x, A, i9)', 'Total root function calls =', ngevals + print '(4x, A, i9)', 'Total LU function calls =', nluevals + print '(4x, A, i9)', 'Num error test failures =', netfails + print '(4x, A, i9)', 'Num nonlinear solver iters =', nniters + print '(4x, A, i9)', 'Num nonlinear solver fails =', nncfails print *, ' ' return From d6ba05df12e153f2ff9749bed2a77fd5bdfe2932 Mon Sep 17 00:00:00 2001 From: rs028 Date: Wed, 7 Jan 2026 15:22:48 +0000 Subject: [PATCH 06/11] URL for numdiff --- tools/install/install_numdiff.sh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tools/install/install_numdiff.sh b/tools/install/install_numdiff.sh index 315501b8..517f75be 100755 --- a/tools/install/install_numdiff.sh +++ b/tools/install/install_numdiff.sh @@ -1,12 +1,12 @@ #!/bin/sh # ----------------------------------------------------------------------------- # -# Copyright (c) 2017 Sam Cox, Roberto Sommariva +# Copyright (c) 2017-2025 Sam Cox, Roberto Sommariva # # This file is part of the AtChem2 software package. # -# This file is covered by the MIT license which can be found in the file -# LICENSE.md at the top level of the AtChem2 distribution. +# This file is licensed under the MIT license, which can be found in the file +# `LICENSE` at the top level of the AtChem2 distribution. # # ----------------------------------------------------------------------------- @@ -42,7 +42,7 @@ fi # download archive NUMDIFF_DIR="numdiff-${NUMDIFF_VERSION}" NUMDIFF_ARCHIVE="${NUMDIFF_DIR}.tar.gz" -wget "https://savannah.nongnu.org/download/numdiff/${NUMDIFF_ARCHIVE}" +wget "https://download.savannah.gnu.org/releases/numdiff/${NUMDIFF_ARCHIVE}" if [ $? -ne 0 ] ; then printf "\n[numdiff] wget --> FAIL\n" exit 1 From 7b6641daddccac671fb7be3a03119d2bcfd10ff9 Mon Sep 17 00:00:00 2001 From: rs028 Date: Wed, 7 Jan 2026 15:27:26 +0000 Subject: [PATCH 07/11] Temporary deactivation of indent test --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1b436b75..8c1f5d39 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -141,7 +141,7 @@ jobs: - name: Indent and Style tests run: | - make indenttest + #make indenttest #temporaarily deactivated make styletest - name: Unit and Model tests From 7e6e1629ce26734efad092a7676cd390d47d0410 Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Wed, 14 Jan 2026 10:21:35 +0000 Subject: [PATCH 08/11] indent test fix --- src/atchem2.f90 | 1288 +++++++++++++++++++++++------------------------ 1 file changed, 644 insertions(+), 644 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 2d390831..33bd8aed 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -23,592 +23,592 @@ module cvode_rhs_mod implicit none type, bind(C) :: UserData - integer(c_int) :: ipar(10) - real(c_double) :: rpar(1) - end type UserData + integer(c_int) :: ipar(10) + real(c_double) :: rpar(1) +end type UserData end module cvode_rhs_mod PROGRAM ATCHEM2 - use, intrinsic :: iso_fortran_env, only : stderr => error_unit - use iso_c_binding - use types_mod - use species_mod - use constraints_mod - use interpolation_method_mod - use reaction_structure_mod - use photolysis_rates_mod, only : photoX, photoY, photoNumberOfPoints, PR_type - use reaction_rates_mod - use env_vars_mod - use date_mod, only : calcInitialDateParameters, calcCurrentDateParameters - use directories_mod, only : output_dir, configuration_dir, shared_library - use storage_mod, only : maxSpecLength, maxPhotoRateNameLength, maxFilepathLength - use solver_params_mod - use model_params_mod - use argparse_mod - use input_functions_mod - use config_functions_mod - use output_functions_mod - use constraint_functions_mod, only : addConstrainedSpeciesToProbSpec, removeConstrainedSpeciesFromProbSpec - use solver_functions_mod, only : jfy, proc - use cvode_rhs_mod, only : UserData - - ! SUNDIALS module +use, intrinsic :: iso_fortran_env, only : stderr => error_unit +use iso_c_binding +use types_mod +use species_mod +use constraints_mod +use interpolation_method_mod +use reaction_structure_mod +use photolysis_rates_mod, only : photoX, photoY, photoNumberOfPoints, PR_type +use reaction_rates_mod +use env_vars_mod +use date_mod, only : calcInitialDateParameters, calcCurrentDateParameters +use directories_mod, only : output_dir, configuration_dir, shared_library +use storage_mod, only : maxSpecLength, maxPhotoRateNameLength, maxFilepathLength +use solver_params_mod +use model_params_mod +use argparse_mod +use input_functions_mod +use config_functions_mod +use output_functions_mod +use constraint_functions_mod, only : addConstrainedSpeciesToProbSpec, removeConstrainedSpeciesFromProbSpec +use solver_functions_mod, only : jfy, proc +use cvode_rhs_mod, only : UserData + +! SUNDIALS module +use fsundials_core_mod +use fnvector_serial_mod +use fcvode_mod +use fsunlinsol_spgmr_mod +use fsunlinsol_dense_mod +use fsunmatrix_dense_mod +use fsunmatrix_band_mod +use fsunlinsol_band_mod +implicit none + +! interface to linux API +interface + + function dlopen( filename, mode ) bind ( c, name="dlopen" ) + ! void *dlopen(const char *filename, int mode); + use iso_c_binding + type(c_ptr) :: dlopen + character(c_char), intent(in) :: filename(*) + integer(c_int), value :: mode + end function + + function dlsym( handle, name ) bind ( c, name="dlsym" ) + ! void *dlsym(void *handle, const char *name); + use iso_c_binding + type(c_funptr) :: dlsym + type(c_ptr), value :: handle + character(c_char), intent(in) :: name(*) + end function + + function dlclose( handle ) bind ( c, name="dlclose" ) + ! int dlclose(void *handle); + use iso_c_binding + integer(c_int) :: dlclose + type(c_ptr), value :: handle + end function +end interface + +! ***************************************************************** +! DECLARATIONS +! ***************************************************************** + +! Declarations for solver parameters +integer(kind=QI) :: ier +integer :: meth, itmeth, iatol, itask, currentNumTimestep +integer(kind=NPI) :: ipar(10) +integer(kind=NPI) :: neq +real(kind=DP) :: t, tout +real(kind=DP) :: rpar(1) + +! Walltime variables +integer(kind=QI) :: runStart, runEnd, runTime, clockRate +! Number of species and reactions +integer(kind=NPI) :: numSpec, numReac + +! Declarations for detailed rates output +type(reaction_frequency_pair) :: invalid_reaction_frequency_pair +type(reaction_frequency_pair), allocatable :: reacDetailedRatesSpecies(:,:), prodDetailedRatesSpecies(:,:) +integer(kind=NPI), allocatable :: reacDetailedRatesSpeciesLengths(:), prodDetailedRatesSpeciesLengths(:) +integer(kind=NPI), allocatable :: detailedRatesSpecies(:) +character(len=maxSpecLength), allocatable :: detailedRatesSpeciesName(:) +! Declarations for concentration outputs +character(len=maxSpecLength), allocatable :: speciesOfInterest(:) +! simulation output time variables +integer(kind=QI) :: time, elapsed +! species concentrations with and without constrained species +real(kind=DP), allocatable :: speciesConcs(:) +real(kind=DP), allocatable :: z(:) + +character(len=400) :: fmt + +type(c_ptr) :: handle +character(len=maxFilepathLength) :: library +type(c_funptr) :: proc_addr +integer :: closure +integer(c_int), parameter :: rtld_lazy=1 ! value extracted from the C header file +integer(c_int), parameter :: rtld_now=2 ! value extracted from the C header file + +! SUNDIALS declarations +type(c_ptr) :: ctx ! SUNDIALS context for the simulation +type(N_Vector), pointer :: sunvec_u ! sundials vector +type(SUNLinearSolver), pointer :: sunls ! sundials linear solver +type(SUNMatrix), pointer :: sunmat_A ! sundials matrix (empty) +type(c_ptr) :: cvode_mem ! CVODE memory +real(c_double) :: t_arr(1) + +type(UserData), target :: udata + +! ***************************************************************** +! Explicit declaration of FCVFUN() interface, which is a +! user-supplied function to CVODE. +interface + + subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) + use types_mod + use species_mod + use constraints_mod + use reaction_structure_mod + use interpolation_functions_mod, only : getVariableConstrainedSpeciesConcentrationAtT + use constraint_functions_mod + + ! Fortran routine for right-hand side function. + real(kind=DP), intent(in) :: t, y(*) + real(kind=DP), intent(out) :: ydot(*) + integer(kind=NPI), intent(in) :: ipar(*) + real(kind=DP), intent(in) :: rpar(*) + integer(kind=NPI), intent(out) :: ier + + integer(kind=NPI) :: nConSpec, np, numReac + real(kind=DP) :: concAtT, dummy + real(kind=DP), allocatable :: dy(:), z(:) + integer(kind=NPI) :: i + end subroutine FCVFUN + + integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) + use, intrinsic :: iso_c_binding use fsundials_core_mod use fnvector_serial_mod - use fcvode_mod - use fsunlinsol_spgmr_mod - use fsunlinsol_dense_mod - use fsunmatrix_dense_mod - use fsunmatrix_band_mod - use fsunlinsol_band_mod - implicit none - ! interface to linux API - interface - - function dlopen( filename, mode ) bind ( c, name="dlopen" ) - ! void *dlopen(const char *filename, int mode); - use iso_c_binding - type(c_ptr) :: dlopen - character(c_char), intent(in) :: filename(*) - integer(c_int), value :: mode - end function - - function dlsym( handle, name ) bind ( c, name="dlsym" ) - ! void *dlsym(void *handle, const char *name); - use iso_c_binding - type(c_funptr) :: dlsym - type(c_ptr), value :: handle - character(c_char), intent(in) :: name(*) - end function - - function dlclose( handle ) bind ( c, name="dlclose" ) - ! int dlclose(void *handle); - use iso_c_binding - integer(c_int) :: dlclose - type(c_ptr), value :: handle - end function - end interface - - ! ***************************************************************** - ! DECLARATIONS - ! ***************************************************************** - - ! Declarations for solver parameters - integer(kind=QI) :: ier - integer :: meth, itmeth, iatol, itask, currentNumTimestep - integer(kind=NPI) :: ipar(10) - integer(kind=NPI) :: neq - real(kind=DP) :: t, tout - real(kind=DP) :: rpar(1) - - ! Walltime variables - integer(kind=QI) :: runStart, runEnd, runTime, clockRate - ! Number of species and reactions - integer(kind=NPI) :: numSpec, numReac - - ! Declarations for detailed rates output - type(reaction_frequency_pair) :: invalid_reaction_frequency_pair - type(reaction_frequency_pair), allocatable :: reacDetailedRatesSpecies(:,:), prodDetailedRatesSpecies(:,:) - integer(kind=NPI), allocatable :: reacDetailedRatesSpeciesLengths(:), prodDetailedRatesSpeciesLengths(:) - integer(kind=NPI), allocatable :: detailedRatesSpecies(:) - character(len=maxSpecLength), allocatable :: detailedRatesSpeciesName(:) - ! Declarations for concentration outputs - character(len=maxSpecLength), allocatable :: speciesOfInterest(:) - ! simulation output time variables - integer(kind=QI) :: time, elapsed - ! species concentrations with and without constrained species - real(kind=DP), allocatable :: speciesConcs(:) - real(kind=DP), allocatable :: z(:) - - character(len=400) :: fmt - - type(c_ptr) :: handle - character(len=maxFilepathLength) :: library - type(c_funptr) :: proc_addr - integer :: closure - integer(c_int), parameter :: rtld_lazy=1 ! value extracted from the C header file - integer(c_int), parameter :: rtld_now=2 ! value extracted from the C header file - - ! SUNDIALS declarations - type(c_ptr) :: ctx ! SUNDIALS context for the simulation - type(N_Vector), pointer :: sunvec_u ! sundials vector - type(SUNLinearSolver), pointer :: sunls ! sundials linear solver - type(SUNMatrix), pointer :: sunmat_A ! sundials matrix (empty) - type(c_ptr) :: cvode_mem ! CVODE memory - real(c_double) :: t_arr(1) - - type(UserData), target :: udata - - ! ***************************************************************** - ! Explicit declaration of FCVFUN() interface, which is a - ! user-supplied function to CVODE. - interface - - subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) - use types_mod - use species_mod - use constraints_mod - use reaction_structure_mod - use interpolation_functions_mod, only : getVariableConstrainedSpeciesConcentrationAtT - use constraint_functions_mod - - ! Fortran routine for right-hand side function. - real(kind=DP), intent(in) :: t, y(*) - real(kind=DP), intent(out) :: ydot(*) - integer(kind=NPI), intent(in) :: ipar(*) - real(kind=DP), intent(in) :: rpar(*) - integer(kind=NPI), intent(out) :: ier - - integer(kind=NPI) :: nConSpec, np, numReac - real(kind=DP) :: concAtT, dummy - real(kind=DP), allocatable :: dy(:), z(:) - integer(kind=NPI) :: i - end subroutine FCVFUN + implicit none - integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) - use, intrinsic :: iso_c_binding - use fsundials_core_mod - use fnvector_serial_mod - - implicit none - - real(c_double), value, intent(in) :: t - type(N_Vector), intent(inout) :: y - type(N_Vector), intent(inout) :: ydot - type(c_ptr), value, intent(in) :: user_data - end function rhs_fn - - end interface - - ! ***************************************************************** - ! MODEL SETUP AND CONFIGURATION - ! ***************************************************************** - - call SYSTEM_CLOCK( runStart ) - - ! Initialise some variables used by CVODE functions to invalid values - ipar(:) = -1_NPI - rpar(:) = -1.0_DP - - write (*, '(A)') 'AtChem2 v1.3-dev' - write (*,*) - write (*, '(A)') '-------------' - write (*, '(A)') ' Directories' - write (*, '(A)') '-------------' - call get_and_set_directories_from_command_arguments() - write (*,*) - - ! Open files for output - open (unit=50, file=trim( output_dir ) // "/speciesConcentrations.output") - open (unit=51, file=trim( output_dir ) // "/errors.output") - open (unit=52, file=trim( output_dir ) // "/environmentVariables.output") - open (unit=53, file=trim( output_dir ) // "/finalModelState.output") - open (unit=55, file=trim( output_dir ) // "/jacobian.output") - open (unit=56, file=trim( output_dir ) // "/lossRates.output") - open (unit=57, file=trim( output_dir ) // "/mainSolverParameters.output") - open (unit=58, file=trim( output_dir ) // "/photolysisRates.output") - open (unit=59, file=trim( output_dir ) // "/photolysisRatesParameters.output") - open (unit=60, file=trim( output_dir ) // "/productionRates.output") - flush(6) - - library = trim( shared_library ) - handle = dlopen(trim( library )//c_null_char, RTLD_LAZY) - if (.not. c_associated(handle)) then - write(*, '(2A)') 'Unable to load DLL ', trim( library ) - stop - end if + real(c_double), value, intent(in) :: t + type(N_Vector), intent(inout) :: y + type(N_Vector), intent(inout) :: ydot + type(c_ptr), value, intent(in) :: user_data +end function rhs_fn - proc_addr=dlsym( handle, "update_p"//c_null_char ) - if ( .not. c_associated(proc_addr) ) then - write(*,*) 'Unable to load the procedure update_p' - stop - end if - call c_f_procpointer( proc_addr, proc ) - - write (*, '(A)') '-----------------------' - write (*, '(A)') ' Species and reactions' - write (*, '(A)') '-----------------------' - call readNumberOfSpeciesAndReactions() - write (*,*) - numSpec = getNumberOfSpecies() - numReac = getNumberOfReactions() - - ! Set array sizes = number of species - allocate (speciesConcs(numSpec), z(numSpec)) - ! Set array size = number of reactions - allocate (reactionRates(numReac)) - - ! Read in reactions - call readReactions() - write (*,*) - neq = numSpec - - ! Read species names and numbers - call setSpeciesList( readSpecies() ) - write (*,*) - - ! Set initial concentrations for all species for which this info is - ! provided, 0.0 for any unspecified. - call readAndSetInitialConcentrations( speciesConcs ) - write (*,*) - - write (*, '(A)') '----------------------------------------' - write (*, '(A)') ' Species requiring detailed rate output' - write (*, '(A)') '----------------------------------------' - - ! Read in species requiring detailed reats output, and set up variables to - ! hold these - write (*, '(A)') ' Reading which species require detailed rate output...' - call readDetailedRatesSpeciesNames( trim( configuration_dir ) // '/outputRates.config', detailedRatesSpeciesName ) - write (*, '(A)') ' Finished reading which species require detailed rate output.' - - allocate (detailedRatesSpecies(size( detailedRatesSpeciesName ))) - allocate (reacDetailedRatesSpecies(size( detailedRatesSpeciesName ), size( clhs, 2 ))) - allocate (prodDetailedRatesSpecies(size( detailedRatesSpeciesName ), size( crhs, 2 ))) - invalid_reaction_frequency_pair = reaction_frequency_pair(-1_NPI, 0_NPI, 0.0_DP) - reacDetailedRatesSpecies(:,:) = invalid_reaction_frequency_pair - prodDetailedRatesSpecies(:,:) = invalid_reaction_frequency_pair - allocate (reacDetailedRatesSpeciesLengths(size( detailedRatesSpeciesName ))) - allocate (prodDetailedRatesSpeciesLengths(size( detailedRatesSpeciesName ))) - - ! Fill detailedRatesSpecies with a list of the numbers of the - ! species requiring detailed rates output, with numbers from their ordering in - ! speciesList - call matchNameToNumber( getSpeciesList(), detailedRatesSpeciesName, detailedRatesSpecies ) - - ! reac/prodDetailedRatesSpecies will each eventually hold one row per species - ! requiring detailed rate output, with the first element being the number of that - ! species, and the remaining elements being the numbers of the - ! reactions in which that species appears as a reactant/product respectively - ! - ! Fill the remaining elements of each row of reac/prodDetailedRatesSpecies with the - ! numbers of the reactions in which that species appears as a reactant/product respectively - call findReactionsWithProductOrReactant( detailedRatesSpecies, clhs, clcoeff, reacDetailedRatesSpecies, & +end interface + +! ***************************************************************** +! MODEL SETUP AND CONFIGURATION +! ***************************************************************** + +call SYSTEM_CLOCK( runStart ) + +! Initialise some variables used by CVODE functions to invalid values +ipar(:) = -1_NPI +rpar(:) = -1.0_DP + +write (*, '(A)') 'AtChem2 v1.3-dev' +write (*,*) +write (*, '(A)') '-------------' +write (*, '(A)') ' Directories' +write (*, '(A)') '-------------' +call get_and_set_directories_from_command_arguments() +write (*,*) + +! Open files for output +open (unit=50, file=trim( output_dir ) // "/speciesConcentrations.output") +open (unit=51, file=trim( output_dir ) // "/errors.output") +open (unit=52, file=trim( output_dir ) // "/environmentVariables.output") +open (unit=53, file=trim( output_dir ) // "/finalModelState.output") +open (unit=55, file=trim( output_dir ) // "/jacobian.output") +open (unit=56, file=trim( output_dir ) // "/lossRates.output") +open (unit=57, file=trim( output_dir ) // "/mainSolverParameters.output") +open (unit=58, file=trim( output_dir ) // "/photolysisRates.output") +open (unit=59, file=trim( output_dir ) // "/photolysisRatesParameters.output") +open (unit=60, file=trim( output_dir ) // "/productionRates.output") +flush(6) + +library = trim( shared_library ) +handle = dlopen(trim( library )//c_null_char, RTLD_LAZY) +if (.not. c_associated(handle)) then +write(*, '(2A)') 'Unable to load DLL ', trim( library ) +stop +end if + +proc_addr=dlsym( handle, "update_p"//c_null_char ) +if ( .not. c_associated(proc_addr) ) then +write(*,*) 'Unable to load the procedure update_p' +stop +end if +call c_f_procpointer( proc_addr, proc ) + +write (*, '(A)') '-----------------------' +write (*, '(A)') ' Species and reactions' +write (*, '(A)') '-----------------------' +call readNumberOfSpeciesAndReactions() +write (*,*) +numSpec = getNumberOfSpecies() +numReac = getNumberOfReactions() + +! Set array sizes = number of species +allocate (speciesConcs(numSpec), z(numSpec)) +! Set array size = number of reactions +allocate (reactionRates(numReac)) + +! Read in reactions +call readReactions() +write (*,*) +neq = numSpec + +! Read species names and numbers +call setSpeciesList( readSpecies() ) +write (*,*) + +! Set initial concentrations for all species for which this info is +! provided, 0.0 for any unspecified. +call readAndSetInitialConcentrations( speciesConcs ) +write (*,*) + +write (*, '(A)') '----------------------------------------' +write (*, '(A)') ' Species requiring detailed rate output' +write (*, '(A)') '----------------------------------------' + +! Read in species requiring detailed reats output, and set up variables to +! hold these +write (*, '(A)') ' Reading which species require detailed rate output...' +call readDetailedRatesSpeciesNames( trim( configuration_dir ) // '/outputRates.config', detailedRatesSpeciesName ) +write (*, '(A)') ' Finished reading which species require detailed rate output.' + +allocate (detailedRatesSpecies(size( detailedRatesSpeciesName ))) +allocate (reacDetailedRatesSpecies(size( detailedRatesSpeciesName ), size( clhs, 2 ))) +allocate (prodDetailedRatesSpecies(size( detailedRatesSpeciesName ), size( crhs, 2 ))) +invalid_reaction_frequency_pair = reaction_frequency_pair(-1_NPI, 0_NPI, 0.0_DP) +reacDetailedRatesSpecies(:,:) = invalid_reaction_frequency_pair +prodDetailedRatesSpecies(:,:) = invalid_reaction_frequency_pair +allocate (reacDetailedRatesSpeciesLengths(size( detailedRatesSpeciesName ))) +allocate (prodDetailedRatesSpeciesLengths(size( detailedRatesSpeciesName ))) + +! Fill detailedRatesSpecies with a list of the numbers of the +! species requiring detailed rates output, with numbers from their ordering in +! speciesList +call matchNameToNumber( getSpeciesList(), detailedRatesSpeciesName, detailedRatesSpecies ) + +! reac/prodDetailedRatesSpecies will each eventually hold one row per species +! requiring detailed rate output, with the first element being the number of that +! species, and the remaining elements being the numbers of the +! reactions in which that species appears as a reactant/product respectively +! +! Fill the remaining elements of each row of reac/prodDetailedRatesSpecies with the +! numbers of the reactions in which that species appears as a reactant/product respectively +call findReactionsWithProductOrReactant( detailedRatesSpecies, clhs, clcoeff, reacDetailedRatesSpecies, & reacDetailedRatesSpeciesLengths ) - call findReactionsWithProductOrReactant( detailedRatesSpecies, crhs, crcoeff, prodDetailedRatesSpecies, & +call findReactionsWithProductOrReactant( detailedRatesSpecies, crhs, crcoeff, prodDetailedRatesSpecies, & prodDetailedRatesSpeciesLengths ) - write (*, '(A, I0)') ' Species requiring detailed rate output (number of species found): ', size( detailedRatesSpeciesName ) - write (*,*) - - ! Read in the numbers of RO2 species from mechanism.ro2 - call readRO2species() - - ! Read in and set solver parameters - call set_solver_parameters( getParametersFromFile( trim( configuration_dir ) // "/solver.parameters" ) ) - write (*,*) - - ! Read in and set model parameters - call set_model_parameters( getParametersFromFile( trim( configuration_dir ) // "/model.parameters" ) ) - write (*,*) - - ! Set the parameters of MODULE date_mod to their value based on - ! startDay, startMonth, startYear - call calcInitialDateParameters() - - ! Hard coded solver parameters - t = modelStartTime - call calcCurrentDateParameters( t ) - tout = timestepSize + t - ! Parameters for FCVodeCreate(). Adams (nonstiff) or BDF (stiff) - meth = CV_BDF - ! itmeth specifies the nonlinear iteration method: 1 for functional - ! iteration or 2 for Newton iteration. - itmeth = 2 - ! IATOL specifies the type for absolute tolerance ATOL: 1 for scalar - ! or 2 for array. If IATOL= 3, the arguments RTOL and ATOL are - ! ignored and the user is expected to subsequently call FCVEWTSET() - ! and provide the function FCVEWT(). - iatol = 1 - - ! Parameter for FCVODE(). Comment from cvode guide: ITASK is a task - ! indicator and should be set to 1 for normal mode (overshoot TOUT - ! and interpolate), or to 2 for one-step mode (return after each - ! internal step taken) - itask = 1 - - ! currentNumTimestep counts the number of iterative steps. Set to - ! zero. Calculation will terminate when - ! currentNumTimestep>=maxNumTimesteps. - currentNumTimestep = 0 - - ! fill speciesOfInterest with the names of species to output to - ! concentration.output - write (*, '(A)') '---------------------' - write (*, '(A)') ' Species of Interest' - write (*, '(A)') '---------------------' - speciesOfInterest = readSpeciesOfInterest() - write (*,*) - - flush(stderr) - - ! ***************************************************************** - ! SET PHOTOLYSIS RATES - ! ***************************************************************** - - write (*, '(A)') '------------' - write (*, '(A)') ' Photolysis' - write (*, '(A)') '------------' - call readPhotoRates() - write (*,*) - - ! Read in environment variables (FIXED, CONSTRAINED, CALC or - ! NOTUSED, see environmentVariables.config) - write (*, '(A)') '-----------------------' - write (*, '(A)') ' Environment variables' - write (*, '(A)') '-----------------------' - call readEnvVar() - write (*,*) - - ! ***************************************************************** - ! SET CONSTRAINTS - ! ***************************************************************** - - write (*, '(A)') '-------------' - write (*, '(A)') ' Constraints' - write (*, '(A)') '-------------' - call readSpeciesConstraints( t, speciesConcs ) - write (*,*) - - call outputSpeciesOfInterest( t, speciesOfInterest, speciesConcs ) - - ! This outputs z, which is speciesConcs with all the constrained - ! species removed. - call removeConstrainedSpeciesFromProbSpec( speciesConcs, getConstrainedSpecies(), z ) - - ! ADJUST PROBLEM SPECIFICATION TO GIVE NUMBER OF SPECIES TO BE - ! SOLVED FOR (N - C = M) - neq = numSpec - getNumberOfConstrainedSpecies() - write (*, '(A)') '---------------' - write (*, '(A)') ' Problem stats' - write (*, '(A)') '---------------' - write (*, '(A30, I0) ') ' neq = ', neq - write (*, '(A30, I0) ') ' numberOfConstrainedSpecies = ', getNumberOfConstrainedSpecies() - - flush(stderr) - - ! ***************************************************************** - ! CONFIGURE SOLVER - ! ***************************************************************** - - ! create the SUNDIALS context - ier = FSUNContext_Create(SUN_COMM_NULL, ctx) - ipar(1) = neq - ipar(2) = numReac - - ! create SUNDIALS N_Vector - sunvec_u => FN_VMake_Serial(neq, z, ctx) - if (.not. associated(sunvec_u)) then - print *, 'ERROR: sunvec = NULL' - stop 1 +write (*, '(A, I0)') ' Species requiring detailed rate output (number of species found): ', size( detailedRatesSpeciesName ) +write (*,*) + +! Read in the numbers of RO2 species from mechanism.ro2 +call readRO2species() + +! Read in and set solver parameters +call set_solver_parameters( getParametersFromFile( trim( configuration_dir ) // "/solver.parameters" ) ) +write (*,*) + +! Read in and set model parameters +call set_model_parameters( getParametersFromFile( trim( configuration_dir ) // "/model.parameters" ) ) +write (*,*) + +! Set the parameters of MODULE date_mod to their value based on +! startDay, startMonth, startYear +call calcInitialDateParameters() + +! Hard coded solver parameters +t = modelStartTime +call calcCurrentDateParameters( t ) +tout = timestepSize + t +! Parameters for FCVodeCreate(). Adams (nonstiff) or BDF (stiff) +meth = CV_BDF +! itmeth specifies the nonlinear iteration method: 1 for functional +! iteration or 2 for Newton iteration. +itmeth = 2 +! IATOL specifies the type for absolute tolerance ATOL: 1 for scalar +! or 2 for array. If IATOL= 3, the arguments RTOL and ATOL are +! ignored and the user is expected to subsequently call FCVEWTSET() +! and provide the function FCVEWT(). +iatol = 1 + +! Parameter for FCVODE(). Comment from cvode guide: ITASK is a task +! indicator and should be set to 1 for normal mode (overshoot TOUT +! and interpolate), or to 2 for one-step mode (return after each +! internal step taken) +itask = 1 + +! currentNumTimestep counts the number of iterative steps. Set to +! zero. Calculation will terminate when +! currentNumTimestep>=maxNumTimesteps. +currentNumTimestep = 0 + +! fill speciesOfInterest with the names of species to output to +! concentration.output +write (*, '(A)') '---------------------' +write (*, '(A)') ' Species of Interest' +write (*, '(A)') '---------------------' +speciesOfInterest = readSpeciesOfInterest() +write (*,*) + +flush(stderr) + +! ***************************************************************** +! SET PHOTOLYSIS RATES +! ***************************************************************** + +write (*, '(A)') '------------' +write (*, '(A)') ' Photolysis' +write (*, '(A)') '------------' +call readPhotoRates() +write (*,*) + +! Read in environment variables (FIXED, CONSTRAINED, CALC or +! NOTUSED, see environmentVariables.config) +write (*, '(A)') '-----------------------' +write (*, '(A)') ' Environment variables' +write (*, '(A)') '-----------------------' +call readEnvVar() +write (*,*) + +! ***************************************************************** +! SET CONSTRAINTS +! ***************************************************************** + +write (*, '(A)') '-------------' +write (*, '(A)') ' Constraints' +write (*, '(A)') '-------------' +call readSpeciesConstraints( t, speciesConcs ) +write (*,*) + +call outputSpeciesOfInterest( t, speciesOfInterest, speciesConcs ) + +! This outputs z, which is speciesConcs with all the constrained +! species removed. +call removeConstrainedSpeciesFromProbSpec( speciesConcs, getConstrainedSpecies(), z ) + +! ADJUST PROBLEM SPECIFICATION TO GIVE NUMBER OF SPECIES TO BE +! SOLVED FOR (N - C = M) +neq = numSpec - getNumberOfConstrainedSpecies() +write (*, '(A)') '---------------' +write (*, '(A)') ' Problem stats' +write (*, '(A)') '---------------' +write (*, '(A30, I0) ') ' neq = ', neq +write (*, '(A30, I0) ') ' numberOfConstrainedSpecies = ', getNumberOfConstrainedSpecies() + +flush(stderr) + +! ***************************************************************** +! CONFIGURE SOLVER +! ***************************************************************** + +! create the SUNDIALS context +ier = FSUNContext_Create(SUN_COMM_NULL, ctx) +ipar(1) = neq +ipar(2) = numReac + +! create SUNDIALS N_Vector +sunvec_u => FN_VMake_Serial(neq, z, ctx) +if (.not. associated(sunvec_u)) then +print *, 'ERROR: sunvec = NULL' +stop 1 +end if + + + +write (*, '(A30, 1P e15.3) ') ' t0 = ', t +write (*,*) + +! create and initialize CVode memory +cvode_mem = FCVodeCreate(meth, ctx) +if (.not. c_associated(cvode_mem)) print *, 'ERROR: cvode_mem = NULL' + +ier = FCVodeInit(cvode_mem, c_funloc(rhs_fn), t, sunvec_u) +if (ier /= 0) then +print *, 'Error in FCVodeInit, ierr = ', ier, '; halting' +stop 1 +end if + +ier = FCVodeSStolerances(cvode_mem, rtol, atol) +if (ier /= 0) then +print *, 'Error in FCVodeSStolerances, ierr = ', ier, '; halting' +stop 1 +end if + +ier = FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) +if (ier /= 0) then +print *, 'Error in FCVodeSetMaxNumSteps, ierr = ', ier, '; halting' +stop 1 +end if + +ier = FCVodeSetMaxStep(cvode_mem, real(maxStep, kind=C_DOUBLE)) +write (*, '(A, I0)') ' setting maxstep ier = ', ier +write (*,*) + +udata%ipar = ipar +udata%rpar = rpar + +ier = FCVodeSetUserData(cvode_mem, c_loc(udata)) +if (ier /= 0) stop 'User data setup failed' + +! SELECT SOLVER TYPE ACCORDING TO FILE INPUT +! SPGMR SOLVER +if ( solverType == 1 ) then +sunls => FSUNLinSol_SPGMR(sunvec_u, SUN_PREC_NONE, int(lookBack, kind=C_INT), ctx) +! SPGMR SOLVER WITH BANDED PRECONDITIONER +else if ( solverType == 2 ) then +sunmat_A => FSUNBandMatrix(neq, preconBandUpper, preconBandLower, ctx) +sunls => FSUNLinSol_Band(sunvec_u, sunmat_A, ctx) +! DENSE SOLVER +else if ( solverType == 3 ) then +! Create dense SUNMatrix for use in linear solves +sunmat_A => FSUNDenseMatrix(neq, neq, ctx) +sunls => FSUNLinSol_Dense(sunvec_u, sunmat_A, ctx) +! UNEXPECTED SOLVER TYPE +else +write (stderr,*) 'Error with solverType input, input = ', solverType +write (stderr,*) 'Available options are 1, 2, 3.' +stop +end if + + +! Attach the matrix and linear solver +ier = FCVodeSetLinearSolver(cvode_mem, sunls, sunmat_A); +! ERROR HANDLING +if ( ier /= 0 ) then +write (stderr,*) ' SUNDIALS_ERROR: FCVodeSetLinearSolver returned ier = ', ier +call FCVodeFree( cvode_mem ) +stop +end if + +! ***************************************************************** +! RUN MODEL +! ***************************************************************** + +write (*, '(A)') '-----------' +write (*, '(A)') ' Model run' +write (*, '(A)') '-----------' + +elapsed = int( t - modelStartTime ) + +do while ( currentNumTimestep < maxNumTimesteps ) + +call calcCurrentDateParameters( t ) + +call outputPhotoRateCalcParameters( t ) + +! Output Jacobian matrix (output frequency set in +! model.parameters) +if ( outputJacobian .eqv. .true. ) then + if ( mod( elapsed, jacobianOutputStepSize ) == 0 ) then + call jfy( numReac, speciesConcs, t ) end if - - - - write (*, '(A30, 1P e15.3) ') ' t0 = ', t - write (*,*) - - ! create and initialize CVode memory - cvode_mem = FCVodeCreate(meth, ctx) - if (.not. c_associated(cvode_mem)) print *, 'ERROR: cvode_mem = NULL' - - ier = FCVodeInit(cvode_mem, c_funloc(rhs_fn), t, sunvec_u) - if (ier /= 0) then - print *, 'Error in FCVodeInit, ierr = ', ier, '; halting' - stop 1 - end if - - ier = FCVodeSStolerances(cvode_mem, rtol, atol) - if (ier /= 0) then - print *, 'Error in FCVodeSStolerances, ierr = ', ier, '; halting' - stop 1 - end if - - ier = FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) - if (ier /= 0) then - print *, 'Error in FCVodeSetMaxNumSteps, ierr = ', ier, '; halting' - stop 1 - end if - - ier = FCVodeSetMaxStep(cvode_mem, real(maxStep, kind=C_DOUBLE)) - write (*, '(A, I0)') ' setting maxstep ier = ', ier - write (*,*) - - udata%ipar = ipar - udata%rpar = rpar - - ier = FCVodeSetUserData(cvode_mem, c_loc(udata)) - if (ier /= 0) stop 'User data setup failed' - - ! SELECT SOLVER TYPE ACCORDING TO FILE INPUT - ! SPGMR SOLVER - if ( solverType == 1 ) then - sunls => FSUNLinSol_SPGMR(sunvec_u, SUN_PREC_NONE, int(lookBack, kind=C_INT), ctx) - ! SPGMR SOLVER WITH BANDED PRECONDITIONER - else if ( solverType == 2 ) then - sunmat_A => FSUNBandMatrix(neq, preconBandUpper, preconBandLower, ctx) - sunls => FSUNLinSol_Band(sunvec_u, sunmat_A, ctx) - ! DENSE SOLVER - else if ( solverType == 3 ) then - ! Create dense SUNMatrix for use in linear solves - sunmat_A => FSUNDenseMatrix(neq, neq, ctx) - sunls => FSUNLinSol_Dense(sunvec_u, sunmat_A, ctx) - ! UNEXPECTED SOLVER TYPE - else - write (stderr,*) 'Error with solverType input, input = ', solverType - write (stderr,*) 'Available options are 1, 2, 3.' - stop - end if - - - ! Attach the matrix and linear solver - ier = FCVodeSetLinearSolver(cvode_mem, sunls, sunmat_A); - ! ERROR HANDLING - if ( ier /= 0 ) then - write (stderr,*) ' SUNDIALS_ERROR: FCVodeSetLinearSolver returned ier = ', ier - call FCVodeFree( cvode_mem ) - stop - end if - - ! ***************************************************************** - ! RUN MODEL - ! ***************************************************************** - - write (*, '(A)') '-----------' - write (*, '(A)') ' Model run' - write (*, '(A)') '-----------' - - elapsed = int( t - modelStartTime ) - - do while ( currentNumTimestep < maxNumTimesteps ) - - call calcCurrentDateParameters( t ) - - call outputPhotoRateCalcParameters( t ) - - ! Output Jacobian matrix (output frequency set in - ! model.parameters) - if ( outputJacobian .eqv. .true. ) then - if ( mod( elapsed, jacobianOutputStepSize ) == 0 ) then - call jfy( numReac, speciesConcs, t ) - end if - end if - - ! Get concentrations for unconstrained species - t_arr(1) =t - ier = FCVode(cvode_mem, tout, sunvec_u, t_arr, itask) - if ( ier /= 0 ) then - write (*, '(A, I0)') ' ier POST FCVODE()= ', ier - end if - t = t_arr(1) - flush(6) - - time = nint( t ) - elapsed = time - modelStartTime - - write (*, '(A, I0)') ' time = ', time - - ! Get concentrations for constrained species and add to array for - ! output - call addConstrainedSpeciesToProbSpec( z, getConstrainedConcs(), getConstrainedSpecies(), speciesConcs ) - - ! Output rates of production and loss (output frequency set in - ! model.parameters) - if ( mod( elapsed, ratesOutputStepSize ) == 0 ) then - call outputRates( detailedRatesSpecies, prodDetailedRatesSpecies, prodDetailedRatesSpeciesLengths, t, reactionRates, 1_SI ) - call outputRates( detailedRatesSpecies, reacDetailedRatesSpecies, reacDetailedRatesSpeciesLengths, t, reactionRates, 0_SI ) - end if - - call outputSpeciesOfInterest( t, speciesOfInterest, speciesConcs ) - call outputPhotolysisRates( t ) - - ! Output reaction rates - if ( mod( elapsed, irOutStepSize ) == 0 ) then - call outputreactionRates( time ) - end if - - ! Output envVar values - ro2 = ro2sum( speciesConcs ) - call outputEnvVar( t ) - - ! Error handling - if ( ier < 0 ) then - fmt = "(///' SUNDIALS_ERROR: FCVODE() returned ier = ', I5, /, 'Linear Solver returned ier = ', I5) " - ! free memory - call FCVodeFree( cvode_mem ) - stop - end if - - ! increment time - tout = tout + timestepSize - currentNumTimestep = currentNumTimestep + 1 - - end do - - ! Output final model concentrations, in a usable format for model - ! restart - call outputFinalModelState( getSpeciesList(), speciesConcs ) - - write (*, '(A)') '------------------' - write (*, '(A)') ' Final statistics' - write (*, '(A)') '------------------' - call PrintFinalStats( cvode_mem ) - write (*,*) - - - call SYSTEM_CLOCK( runEnd, clockRate ) - runTime = ( runEnd - runStart ) / clockRate - write (*, '(A, I0)') ' Runtime = ', runTime - write (*, '(A)') ' Deallocating memory.' - - ! ***************************************************************** - ! STOP MODEL - ! ***************************************************************** - - ! deallocate CVODE internal data +end if + +! Get concentrations for unconstrained species +t_arr(1) =t +ier = FCVode(cvode_mem, tout, sunvec_u, t_arr, itask) +if ( ier /= 0 ) then + write (*, '(A, I0)') ' ier POST FCVODE()= ', ier +end if +t = t_arr(1) +flush(6) + +time = nint( t ) +elapsed = time - modelStartTime + +write (*, '(A, I0)') ' time = ', time + +! Get concentrations for constrained species and add to array for +! output +call addConstrainedSpeciesToProbSpec( z, getConstrainedConcs(), getConstrainedSpecies(), speciesConcs ) + +! Output rates of production and loss (output frequency set in +! model.parameters) +if ( mod( elapsed, ratesOutputStepSize ) == 0 ) then + call outputRates( detailedRatesSpecies, prodDetailedRatesSpecies, prodDetailedRatesSpeciesLengths, t, reactionRates, 1_SI ) + call outputRates( detailedRatesSpecies, reacDetailedRatesSpecies, reacDetailedRatesSpeciesLengths, t, reactionRates, 0_SI ) +end if + +call outputSpeciesOfInterest( t, speciesOfInterest, speciesConcs ) +call outputPhotolysisRates( t ) + +! Output reaction rates +if ( mod( elapsed, irOutStepSize ) == 0 ) then + call outputreactionRates( time ) +end if + +! Output envVar values +ro2 = ro2sum( speciesConcs ) +call outputEnvVar( t ) + +! Error handling +if ( ier < 0 ) then + fmt = "(///' SUNDIALS_ERROR: FCVODE() returned ier = ', I5, /, 'Linear Solver returned ier = ', I5) " + ! free memory call FCVodeFree( cvode_mem ) - deallocate (speciesConcs, z) - deallocate (reacDetailedRatesSpecies, prodDetailedRatesSpecies) - deallocate (detailedRatesSpeciesName, speciesOfInterest) - deallocate (reactionRates) - deallocate (clhs, clcoeff, crhs, crcoeff) - deallocate (reacDetailedRatesSpeciesLengths, prodDetailedRatesSpeciesLengths) - - ! deallocate data allocated in inputFunctions.f90 - ! deallocate arrays from module constraints_mod - call deallocateConstrainedConcs() - call deallocateConstrainedSpecies() - deallocate (dataX, dataY, dataFixedY) - deallocate (speciesNumberOfPoints) - - ! deallocate arrays from module species_mod - call deallocateSpeciesList() - - ! deallocate arrays from module env_vars_mod - deallocate (envVarTypesNum, envVarNames, envVarTypes, envVarFixedValues) - deallocate (envVarX, envVarY, envVarNumberOfPoints) - - ! deallocate arrays from module photolysis_rates_mod - if ( PR_type == 2 .or. PR_type == 3) then - deallocate (photoX, photoY, photoNumberOfPoints) - end if - - ! Close output files and end program - close (50) - close (51) - close (52) - close (53) - close (55) - close (56) - close (57) - close (58) - close (59) - close (60) - closure=dlclose(handle) - stop +end if + +! increment time +tout = tout + timestepSize +currentNumTimestep = currentNumTimestep + 1 + +end do + +! Output final model concentrations, in a usable format for model +! restart +call outputFinalModelState( getSpeciesList(), speciesConcs ) + +write (*, '(A)') '------------------' +write (*, '(A)') ' Final statistics' +write (*, '(A)') '------------------' +call PrintFinalStats( cvode_mem ) +write (*,*) + + +call SYSTEM_CLOCK( runEnd, clockRate ) +runTime = ( runEnd - runStart ) / clockRate +write (*, '(A, I0)') ' Runtime = ', runTime +write (*, '(A)') ' Deallocating memory.' + +! ***************************************************************** +! STOP MODEL +! ***************************************************************** + +! deallocate CVODE internal data +call FCVodeFree( cvode_mem ) +deallocate (speciesConcs, z) +deallocate (reacDetailedRatesSpecies, prodDetailedRatesSpecies) +deallocate (detailedRatesSpeciesName, speciesOfInterest) +deallocate (reactionRates) +deallocate (clhs, clcoeff, crhs, crcoeff) +deallocate (reacDetailedRatesSpeciesLengths, prodDetailedRatesSpeciesLengths) + +! deallocate data allocated in inputFunctions.f90 +! deallocate arrays from module constraints_mod +call deallocateConstrainedConcs() +call deallocateConstrainedSpecies() +deallocate (dataX, dataY, dataFixedY) +deallocate (speciesNumberOfPoints) + +! deallocate arrays from module species_mod +call deallocateSpeciesList() + +! deallocate arrays from module env_vars_mod +deallocate (envVarTypesNum, envVarNames, envVarTypes, envVarFixedValues) +deallocate (envVarX, envVarY, envVarNumberOfPoints) + +! deallocate arrays from module photolysis_rates_mod +if ( PR_type == 2 .or. PR_type == 3) then +deallocate (photoX, photoY, photoNumberOfPoints) +end if + +! Close output files and end program +close (50) +close (51) +close (52) +close (53) +close (55) +close (56) +close (57) +close (58) +close (59) +close (60) +closure=dlclose(handle) + +stop END PROGRAM ATCHEM2 @@ -621,91 +621,91 @@ END PROGRAM ATCHEM2 ! -------------------------------------------------------- ! ! Fortran routine for right-hand side function. subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) - use types_mod - use constraints_mod, only : getNumberOfConstrainedSpecies, numberOfVariableConstrainedSpecies, dataFixedY, & +use types_mod +use constraints_mod, only : getNumberOfConstrainedSpecies, numberOfVariableConstrainedSpecies, dataFixedY, & getConstrainedSpecies, setConstrainedConcs - use reaction_structure_mod, only : clhs, clcoeff, crhs, crcoeff - use interpolation_method_mod, only : getSpeciesInterpMethod - use interpolation_functions_mod, only : getVariableConstrainedSpeciesConcentrationAtT, getConstrainedPhotoRatesAtT - use constraint_functions_mod, only : addConstrainedSpeciesToProbSpec, removeConstrainedSpeciesFromProbSpec - use solver_functions_mod, only : resid - implicit none - - real(kind=DP), intent(in) :: t, y(*) - real(kind=DP), intent(out) :: ydot(*) - integer(kind=NPI), intent(in) :: ipar(*) - real(kind=DP), intent(in) :: rpar(*) - integer(kind=NPI), intent(out) :: ier - integer(kind=NPI) :: numConSpec, np, numReac, i - real(kind=DP) :: dummy - real(kind=DP), allocatable :: dy(:), z(:), constrainedConcs(:) - - numConSpec = getNumberOfConstrainedSpecies() - np = ipar(1) + numConSpec - numReac = ipar(2) - dummy = rpar(1) - !TODO: these should be the same size on every call? If so, then better to allocate once and re-use - ! rather than re-allocating each time. - allocate (dy(np), z(np), constrainedConcs(numConSpec)) - - ! for each constrained species... - do i = 1, numConSpec - ! if it's a variable-concentration constrained species, - if ( i <= numberOfVariableConstrainedSpecies ) then - call getVariableConstrainedSpeciesConcentrationAtT( t, i, constrainedConcs(i) ) - else - constrainedConcs(i) = dataFixedY(i - numberOfVariableConstrainedSpecies) - end if - end do - - call setConstrainedConcs( constrainedConcs ) - - call addConstrainedSpeciesToProbSpec( y, constrainedConcs, getConstrainedSpecies(), z ) - - call resid( numReac, t, z, dy, clhs, clcoeff, crhs, crcoeff ) - - call removeConstrainedSpeciesFromProbSpec( dy, getConstrainedSpecies(), ydot ) - - deallocate (dy, z) - ier = 0 - - return +use reaction_structure_mod, only : clhs, clcoeff, crhs, crcoeff +use interpolation_method_mod, only : getSpeciesInterpMethod +use interpolation_functions_mod, only : getVariableConstrainedSpeciesConcentrationAtT, getConstrainedPhotoRatesAtT +use constraint_functions_mod, only : addConstrainedSpeciesToProbSpec, removeConstrainedSpeciesFromProbSpec +use solver_functions_mod, only : resid +implicit none + +real(kind=DP), intent(in) :: t, y(*) +real(kind=DP), intent(out) :: ydot(*) +integer(kind=NPI), intent(in) :: ipar(*) +real(kind=DP), intent(in) :: rpar(*) +integer(kind=NPI), intent(out) :: ier +integer(kind=NPI) :: numConSpec, np, numReac, i +real(kind=DP) :: dummy +real(kind=DP), allocatable :: dy(:), z(:), constrainedConcs(:) + +numConSpec = getNumberOfConstrainedSpecies() +np = ipar(1) + numConSpec +numReac = ipar(2) +dummy = rpar(1) +!TODO: these should be the same size on every call? If so, then better to allocate once and re-use +! rather than re-allocating each time. +allocate (dy(np), z(np), constrainedConcs(numConSpec)) + +! for each constrained species... +do i = 1, numConSpec +! if it's a variable-concentration constrained species, +if ( i <= numberOfVariableConstrainedSpecies ) then + call getVariableConstrainedSpeciesConcentrationAtT( t, i, constrainedConcs(i) ) +else + constrainedConcs(i) = dataFixedY(i - numberOfVariableConstrainedSpecies) +end if +end do + +call setConstrainedConcs( constrainedConcs ) + +call addConstrainedSpeciesToProbSpec( y, constrainedConcs, getConstrainedSpecies(), z ) + +call resid( numReac, t, z, dy, clhs, clcoeff, crhs, crcoeff ) + +call removeConstrainedSpeciesFromProbSpec( dy, getConstrainedSpecies(), ydot ) + +deallocate (dy, z) +ier = 0 + +return end subroutine FCVFUN ! ******************************************************************** ! integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) - use, intrinsic :: iso_c_binding - use types_mod - use fsundials_core_mod - use fnvector_serial_mod - use cvode_rhs_mod, only : UserData - use types_mod +use, intrinsic :: iso_c_binding +use types_mod +use fsundials_core_mod +use fnvector_serial_mod +use cvode_rhs_mod, only : UserData +use types_mod - implicit none +implicit none - real(c_double), value, intent(in) :: t - type(N_Vector), intent(inout) :: y - type(N_Vector), intent(inout) :: ydot - type(c_ptr), value, intent(in) :: user_data +real(c_double), value, intent(in) :: t +type(N_Vector), intent(inout) :: y +type(N_Vector), intent(inout) :: ydot +type(c_ptr), value, intent(in) :: user_data - ! Local pointers to vector data - real(c_double), pointer :: ydata(:) - real(c_double), pointer :: ydotdata(:) +! Local pointers to vector data +real(c_double), pointer :: ydata(:) +real(c_double), pointer :: ydotdata(:) - type(UserData), pointer :: ud - integer(kind=NPI) :: ier +type(UserData), pointer :: ud +integer(kind=NPI) :: ier - integer(kind=NPI) :: ipar_f(10) - integer :: i - call c_f_pointer( user_data, ud ) - ipar_f = [(int(ud%ipar(i), kind=NPI), i=1, size(ud%ipar))] +integer(kind=NPI) :: ipar_f(10) +integer :: i +call c_f_pointer( user_data, ud ) +ipar_f = [(int(ud%ipar(i), kind=NPI), i=1, size(ud%ipar))] - ! Get raw data arrays from N_Vector - ydata => FN_VGetArrayPointer(y) - ydotdata => FN_VGetArrayPointer(ydot) +! Get raw data arrays from N_Vector +ydata => FN_VGetArrayPointer(y) +ydotdata => FN_VGetArrayPointer(ydot) - call FCVFUN( t, ydata, ydotdata, ipar_f, ud%rpar, ier ) +call FCVFUN( t, ydata, ydotdata, ipar_f, ud%rpar, ier ) - rhs_fn = ier ! 0 = success +rhs_fn = ier ! 0 = success end function rhs_fn From ce3bb37aea0923f852041b15cbd93c9aa7cbb7ea Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Wed, 14 Jan 2026 15:21:32 +0000 Subject: [PATCH 09/11] restoring output function --- src/atchem2.f90 | 3 ++ src/outputFunctions.f90 | 80 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 33bd8aed..4a858981 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -529,6 +529,9 @@ end function rhs_fn call outputreactionRates( time ) end if +! Output CVODE solver parameters and timestep sizes +call outputSolverParameters( t, cvode_mem, solverType ) + ! Output envVar values ro2 = ro2sum( speciesConcs ) call outputEnvVar( t ) diff --git a/src/outputFunctions.f90 b/src/outputFunctions.f90 index ca5d5b29..fc48265f 100644 --- a/src/outputFunctions.f90 +++ b/src/outputFunctions.f90 @@ -63,6 +63,86 @@ subroutine outputEnvVar( t ) return end subroutine outputEnvVar + + ! ----------------------------------------------------------------- + ! Write parameters output by CVODE solver to file. + subroutine outputSolverParameters( t, cvode_mem, solver_type ) + use, intrinsic :: iso_fortran_env, only : stderr => error_unit + use types_mod + use iso_c_binding + use fcvode_mod + + real(kind=DP), intent(in) :: t + integer(c_int) :: retval ! error flag + type(c_ptr), intent(in) :: cvode_mem ! solver memory structure + integer(kind=SI), intent(in) :: solver_type + + integer(kind=SI) :: i + logical :: first_time = .true. + + real*8 :: hcur(1), hlast(1) + integer(c_long) :: lenrw(1), leniw(1), nst(1), nfe(1), netf(1), nsetups(1), nje(1), nni(1), ncfn(1), nor(1) + integer(c_long) :: lenrwls(1), leniwls(1), lsflag(1), nfels(1), njtv(1), npe(1), nps(1), nli(1), ncfl(1), nfeDQ(1) + integer(c_int) :: qu(1), qcur(1) + retval = FCVodeGetWorkSpace(cvode_mem, lenrw, leniw) + retval = FCVodeGetCurrentStep(cvode_mem, hcur) + retval = FCVodeGetLastStep(cvode_mem, hlast) + retval = FCVodeGetNumSteps(cvode_mem, nst) + retval = FCVodeGetNumRhsEvals(cvode_mem, nfe) + retval = FCVodeGetNumErrTestFails(cvode_mem, netf) + retval = FCVodeGetNumLinSolvSetups(cvode_mem, nsetups) + retval = FCVodeGetNumNonlinSolvIters(cvode_mem, nni) + retval = FCVodeGetNumNonlinSolvConvFails(cvode_mem, ncfn) + retval = FCVodeGetNumStabLimOrderReds(cvode_mem, nor) + retval = FCVodeGetLastOrder(cvode_mem, qu) + retval = FCVodeGetCurrentOrder(cvode_mem, qcur) + retval = FCVodeGetLinWorkSpace(cvode_mem, lenrwls, leniwls) + retval = FCVodeGetNumLinRhsEvals(cvode_mem, nfels) + + if ( ( solver_type == 1 ) .or. ( solver_type == 2 ) ) then + ! CVSPILS type solver + if ( first_time .eqv. .true. ) then + write (57, '(A9, 2A17, 20A9) ') 't', 'currentStepSize', 'previousStepSize', 'LENRW', 'LENIW', 'NST', 'NFE', & + 'NETF', 'NCFN', 'NNI', 'NSETUPS', 'QU', 'QCUR', 'NOR', 'LENRWLS', 'LENIWLS', & + 'LS_FLAG', 'NFELS', 'NJTV', 'NPE', 'NPS', 'NLI', 'NCFL' + first_time = .false. + end if + retval = FCVodeGetNumLinIters(cvode_mem, nli) + retval = FCVodeGetNumLinConvFails(cvode_mem, ncfl) + retval = FCVodeGetNumPrecEvals(cvode_mem, npe) + retval = FCVodeGetNumPrecSolves(cvode_mem, nps) + retval = FCVodeGetNumJtimesEvals(cvode_mem, njtv) + retval = FCVodeGetLastLinFlag(cvode_mem, lsflag) + + write (57, '(1P e9.2, 2 (ES17.8E3), 20I9) ') t, hcur, hlast, lenrw, leniw, nst, nfe, & + netf, ncfn, nni, nsetups, qu, qcur, nor, lenrwls, leniwls, & + lsflag, nfels, njtv, npe, nps, nli, ncfl + + else if ( solver_type == 3 ) then + ! CVDLS type solver + if ( first_time .eqv. .true. ) then + write (57, '(A9, 2A17, 16A9) ') 't', 'currentStepSize', 'previousStepSize', 'LENRW', 'LENIW', 'NST', 'NFE', & + 'NETF', 'NCFN', 'NNI', 'NSETUPS', 'QU', 'QCUR', 'NOR', 'LENRWLS', 'LENIWLS', & + 'LS_FLAG', 'NFELS', 'NJE' + first_time = .false. + end if + retval = FCVodeGetNumJacEvals(cvode_mem, nje) + retval = FCVDiagGetNumRhsEvals(cvode_mem, nfeDQ) + retval = FCVDiagGetLastFlag(cvode_mem, lsflag) + + write (57, '(1P e9.2, 2 (ES17.8E3), 16I9) ') t, hcur, hlast, lenrw, leniw, nst, nfe, & + netf, ncfn, nni, nsetups, qu, qcur, nor, lenrwls, leniwls, & + lsflag, nfels, nje + + else + write (stderr,*) 'outputSolverParameters(): Error with solver_type = ', solver_type + write (stderr,*) 'Available options are 1, 2, 3.' + stop + end if + + return + end subroutine outputSolverParameters + subroutine PrintFinalStats( cvode_mem ) !======= Inclusions =========== From b75ec7d1ad44c8032be9246defaecbbd758164df Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Wed, 14 Jan 2026 15:59:54 +0000 Subject: [PATCH 10/11] changes to the expected output of overview out files --- .../model_tests/env_model_1/env_model_1.out.cmp | 15 +++++++++------ .../model_tests/env_model_2/env_model_2.out.cmp | 15 +++++++++------ .../model_tests/env_model_3/env_model_3.out.cmp | 15 +++++++++------ .../model_tests/env_model_4/env_model_4.out.cmp | 15 +++++++++------ tests/model_tests/firstorder/firstorder.out.cmp | 15 +++++++++------ .../model_tests/secondorder/secondorder.out.cmp | 15 +++++++++------ .../spec_model_1/spec_model_1.out.cmp | 15 +++++++++------ .../spec_model_func/spec_model_func.out.cmp | 16 ++++++++++------ .../spec_model_kpp/spec_model_kpp.out.cmp | 16 ++++++++++------ .../spec_model_stoich/spec_model_stoich.out.cmp | 15 +++++++++------ tests/model_tests/static/static.out.cmp | 15 +++++++++------ 11 files changed, 101 insertions(+), 66 deletions(-) diff --git a/tests/model_tests/env_model_1/env_model_1.out.cmp b/tests/model_tests/env_model_1/env_model_1.out.cmp index 7e6928e0..4984ffa0 100644 --- a/tests/model_tests/env_model_1/env_model_1.out.cmp +++ b/tests/model_tests/env_model_1/env_model_1.out.cmp @@ -161,7 +161,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 5.040E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -191,14 +190,18 @@ AtChem2 v1.3-dev time = 63600 time = 64200 time = 64800 - ------------------ Final statistics ------------------ - No. steps = 380 No. f-s = 502 No. J-s = 514 No. LU-s = 95 - No. nonlinear iterations = 501 - No. nonlinear convergence failures = 0 - No. error test failures = 28 + Total internal steps taken = 358 + Total rhs function calls = 467 + Total Jacobian function calls = 7 + Total root function calls = 0 + Total LU function calls = 94 + Num error test failures = 29 + Num nonlinear solver iters = 466 + Num nonlinear solver fails = 0 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/env_model_2/env_model_2.out.cmp b/tests/model_tests/env_model_2/env_model_2.out.cmp index 5cbb76e1..4807258a 100644 --- a/tests/model_tests/env_model_2/env_model_2.out.cmp +++ b/tests/model_tests/env_model_2/env_model_2.out.cmp @@ -163,7 +163,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 5.040E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -193,14 +192,18 @@ AtChem2 v1.3-dev time = 63600 time = 64200 time = 64800 - ------------------ Final statistics ------------------ - No. steps = 359 No. f-s = 493 No. J-s = 471 No. LU-s = 102 - No. nonlinear iterations = 492 - No. nonlinear convergence failures = 0 - No. error test failures = 36 + Total internal steps taken = 363 + Total rhs function calls = 493 + Total Jacobian function calls = 7 + Total root function calls = 0 + Total LU function calls = 102 + Num error test failures = 37 + Num nonlinear solver iters = 492 + Num nonlinear solver fails = 0 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/env_model_3/env_model_3.out.cmp b/tests/model_tests/env_model_3/env_model_3.out.cmp index 9a2045f9..4b93db77 100644 --- a/tests/model_tests/env_model_3/env_model_3.out.cmp +++ b/tests/model_tests/env_model_3/env_model_3.out.cmp @@ -163,7 +163,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 5.040E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -193,14 +192,18 @@ AtChem2 v1.3-dev time = 63600 time = 64200 time = 64800 - ------------------ Final statistics ------------------ - No. steps = 229 No. f-s = 277 No. J-s = 157 No. LU-s = 32 - No. nonlinear iterations = 273 - No. nonlinear convergence failures = 0 - No. error test failures = 9 + Total internal steps taken = 233 + Total rhs function calls = 276 + Total Jacobian function calls = 4 + Total root function calls = 0 + Total LU function calls = 29 + Num error test failures = 8 + Num nonlinear solver iters = 273 + Num nonlinear solver fails = 0 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/env_model_4/env_model_4.out.cmp b/tests/model_tests/env_model_4/env_model_4.out.cmp index 09d7681b..ba84f5f4 100644 --- a/tests/model_tests/env_model_4/env_model_4.out.cmp +++ b/tests/model_tests/env_model_4/env_model_4.out.cmp @@ -163,7 +163,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 5.040E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -193,14 +192,18 @@ AtChem2 v1.3-dev time = 63600 time = 64200 time = 64800 - ------------------ Final statistics ------------------ - No. steps = 229 No. f-s = 277 No. J-s = 157 No. LU-s = 32 - No. nonlinear iterations = 273 - No. nonlinear convergence failures = 0 - No. error test failures = 9 + Total internal steps taken = 233 + Total rhs function calls = 276 + Total Jacobian function calls = 4 + Total root function calls = 0 + Total LU function calls = 29 + Num error test failures = 8 + Num nonlinear solver iters = 273 + Num nonlinear solver fails = 0 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/firstorder/firstorder.out.cmp b/tests/model_tests/firstorder/firstorder.out.cmp index dc5fd78e..2b566b92 100644 --- a/tests/model_tests/firstorder/firstorder.out.cmp +++ b/tests/model_tests/firstorder/firstorder.out.cmp @@ -160,7 +160,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 0.000E+00 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -1666,14 +1665,18 @@ AtChem2 v1.3-dev time = 1498 time = 1499 time = 1500 - ------------------ Final statistics ------------------ - No. steps = 54 No. f-s = 77 No. J-s = 56 No. LU-s = 17 - No. nonlinear iterations = 73 - No. nonlinear convergence failures = 0 - No. error test failures = 3 + Total internal steps taken = 54 + Total rhs function calls = 76 + Total Jacobian function calls = 2 + Total root function calls = 0 + Total LU function calls = 17 + Num error test failures = 3 + Num nonlinear solver iters = 73 + Num nonlinear solver fails = 0 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/secondorder/secondorder.out.cmp b/tests/model_tests/secondorder/secondorder.out.cmp index 2a027fae..0cdd7b87 100644 --- a/tests/model_tests/secondorder/secondorder.out.cmp +++ b/tests/model_tests/secondorder/secondorder.out.cmp @@ -160,7 +160,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 0.000E+00 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -1666,14 +1665,18 @@ AtChem2 v1.3-dev time = 1498 time = 1499 time = 1500 - ------------------ Final statistics ------------------ - No. steps = 67 No. f-s = 85 No. J-s = 72 No. LU-s = 18 - No. nonlinear iterations = 81 - No. nonlinear convergence failures = 0 - No. error test failures = 1 + Total internal steps taken = 73 + Total rhs function calls = 89 + Total Jacobian function calls = 2 + Total root function calls = 0 + Total LU function calls = 19 + Num error test failures = 1 + Num nonlinear solver iters = 86 + Num nonlinear solver fails = 0 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/spec_model_1/spec_model_1.out.cmp b/tests/model_tests/spec_model_1/spec_model_1.out.cmp index 7b30406d..756cd081 100644 --- a/tests/model_tests/spec_model_1/spec_model_1.out.cmp +++ b/tests/model_tests/spec_model_1/spec_model_1.out.cmp @@ -162,7 +162,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 2.340E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -198,14 +197,18 @@ AtChem2 v1.3-dev time = 48600 time = 49500 time = 50400 - ------------------ Final statistics ------------------ - No. steps = 503 No. f-s = 619 No. J-s = 970 No. LU-s = 108 - No. nonlinear iterations = 618 - No. nonlinear convergence failures = 0 - No. error test failures = 30 + Total internal steps taken = 557 + Total rhs function calls = 723 + Total Jacobian function calls = 12 + Total root function calls = 0 + Total LU function calls = 135 + Num error test failures = 40 + Num nonlinear solver iters = 722 + Num nonlinear solver fails = 3 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/spec_model_func/spec_model_func.out.cmp b/tests/model_tests/spec_model_func/spec_model_func.out.cmp index 59f9215a..ba8871a2 100644 --- a/tests/model_tests/spec_model_func/spec_model_func.out.cmp +++ b/tests/model_tests/spec_model_func/spec_model_func.out.cmp @@ -162,7 +162,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 2.340E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -197,15 +196,20 @@ AtChem2 v1.3-dev time = 47700 time = 48600 time = 49500 +Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL time = 50400 - ------------------ Final statistics ------------------ - No. steps = 536 No. f-s = 663 No. J-s = 1004 No. LU-s = 117 - No. nonlinear iterations = 662 - No. nonlinear convergence failures = 0 - No. error test failures = 32 + Total internal steps taken = 537 + Total rhs function calls = 687 + Total Jacobian function calls = 12 + Total root function calls = 0 + Total LU function calls = 131 + Num error test failures = 37 + Num nonlinear solver iters = 686 + Num nonlinear solver fails = 3 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp b/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp index 65d50ee7..6a5493a7 100644 --- a/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp +++ b/tests/model_tests/spec_model_kpp/spec_model_kpp.out.cmp @@ -162,7 +162,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 2.340E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -197,15 +196,20 @@ AtChem2 v1.3-dev time = 47700 time = 48600 time = 49500 +Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL time = 50400 - ------------------ Final statistics ------------------ - No. steps = 536 No. f-s = 663 No. J-s = 1001 No. LU-s = 117 - No. nonlinear iterations = 662 - No. nonlinear convergence failures = 0 - No. error test failures = 32 + Total internal steps taken = 532 + Total rhs function calls = 666 + Total Jacobian function calls = 12 + Total root function calls = 0 + Total LU function calls = 114 + Num error test failures = 29 + Num nonlinear solver iters = 665 + Num nonlinear solver fails = 3 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp b/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp index d62ee738..ae75e9d0 100644 --- a/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp +++ b/tests/model_tests/spec_model_stoich/spec_model_stoich.out.cmp @@ -162,7 +162,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 2.340E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -198,14 +197,18 @@ AtChem2 v1.3-dev time = 48600 time = 49500 time = 50400 - ------------------ Final statistics ------------------ - No. steps = 503 No. f-s = 619 No. J-s = 947 No. LU-s = 108 - No. nonlinear iterations = 618 - No. nonlinear convergence failures = 0 - No. error test failures = 30 + Total internal steps taken = 524 + Total rhs function calls = 651 + Total Jacobian function calls = 11 + Total root function calls = 0 + Total LU function calls = 110 + Num error test failures = 25 + Num nonlinear solver iters = 650 + Num nonlinear solver fails = 3 + Runtime = 0 Deallocating memory. diff --git a/tests/model_tests/static/static.out.cmp b/tests/model_tests/static/static.out.cmp index fc4af145..6a805e98 100644 --- a/tests/model_tests/static/static.out.cmp +++ b/tests/model_tests/static/static.out.cmp @@ -160,7 +160,6 @@ AtChem2 v1.3-dev numberOfConstrainedSpecies = 0 t0 = 4.320E+04 - setting maxnumsteps ier = 0 setting maxstep ier = 0 ----------- @@ -266,14 +265,18 @@ AtChem2 v1.3-dev time = 2983200 time = 3013200 time = 3043200 - ------------------ Final statistics ------------------ - No. steps = 30001 No. f-s = 30005 No. J-s = 0 No. LU-s = 1501 - No. nonlinear iterations = 30001 - No. nonlinear convergence failures = 0 - No. error test failures = 0 + Total internal steps taken = 30001 + Total rhs function calls = 30004 + Total Jacobian function calls = 500 + Total root function calls = 0 + Total LU function calls = 1501 + Num error test failures = 0 + Num nonlinear solver iters = 30001 + Num nonlinear solver fails = 0 + Runtime = 0 Deallocating memory. From 222fdbd149b15e2ad4cf6a21ec30d92c40158c63 Mon Sep 17 00:00:00 2001 From: Neil Butcher Date: Thu, 15 Jan 2026 08:51:31 +0000 Subject: [PATCH 11/11] hopeful unit tests fix --- tests/makefile.tests | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/makefile.tests b/tests/makefile.tests index 7db9e066..709b4f24 100644 --- a/tests/makefile.tests +++ b/tests/makefile.tests @@ -41,7 +41,7 @@ $(UNITTESTDIR)/fruit_basket_gen.f90 : $(unittest_code) # build `fruit_driver.exe` from the individual unit tests $(fruit_driver) : $(all_unittest_code) - $(FORT_COMP) -o $(fruit_driver) -J$(OBJ) -I$(OBJ) $(all_unittest_code) $(FFLAGS) $(LDFLAGS) + $(FORT_COMP) -o $(fruit_driver) -J$(OBJ) -I$(OBJ) -I$(CVODEOBJDIR) $(all_unittest_code) $(FFLAGS) $(LDFLAGS) # ==================== Model tests ==================== #