From 13e4fc9c3e90594470ffa2d8149acd27f9c2df8d Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Tue, 16 Dec 2025 16:16:58 +0000 Subject: [PATCH 1/9] updating the version of sundials --- docker/install.sh | 2 +- tools/install/Makefile.skel | 4 ++-- tools/install/install_cvode.sh | 3 ++- 3 files changed, 5 insertions(+), 4 deletions(-) mode change 100755 => 100644 docker/install.sh diff --git a/docker/install.sh b/docker/install.sh old mode 100755 new mode 100644 index 1c8d0ef7a..391806f91 --- a/docker/install.sh +++ b/docker/install.sh @@ -48,7 +48,7 @@ cd atchem # Change atchem dependancy paths and create Makefile from skeleton sed 's,cvode/lib,/atchem-lib/cvode/lib,g' tools/install/Makefile.skel > ./Makefile -sed -i 's,openlibm-0.8.1,/atchem-lib/openlibm-0.8.1,g' ./Makefile +sed -i 's,openlibm/lib,/atchem-lib/openlibm/,g' ./Makefile # Fix python command to match installed version sed -i "s/python/python3/g" ./build/build_atchem2.sh diff --git a/tools/install/Makefile.skel b/tools/install/Makefile.skel index 40a2e1ff6..f72b75350 100644 --- a/tools/install/Makefile.skel +++ b/tools/install/Makefile.skel @@ -22,7 +22,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 -OPENLIBMDIR = openlibm +OPENLIBMDIR = openlibm/lib FRUITDIR = fruit_3.4.3 # =========================================================================== # @@ -82,7 +82,7 @@ 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_cvode -lsundials_fnvecserial_mod -lsundials_nvecserial -ldl # ---------------------------------------------------- # EXECUTABLE setup diff --git a/tools/install/install_cvode.sh b/tools/install/install_cvode.sh index 3a9cec862..090ca546f 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" From 8651211d562a44bae5db1c026ee23d01ea621506 Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Tue, 16 Dec 2025 18:15:05 +0000 Subject: [PATCH 2/9] changes to use new sundials --- src/atchem2.f90 | 169 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 144 insertions(+), 25 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 496921665..87859d053 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -41,6 +41,15 @@ PROGRAM ATCHEM2 use output_functions_mod use constraint_functions_mod, only : addConstrainedSpeciesToProbSpec, removeConstrainedSpeciesFromProbSpec use solver_functions_mod, only : jfy, proc + + ! 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 implicit none ! interface to linux API @@ -110,6 +119,23 @@ 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) :: sunctx + type(N_Vector) :: n_vector_z + type(c_ptr) :: cvode_mem + type(SUNLinearSolver) :: LS + type(SUNMatrix) :: A + real(c_double) :: t_arr(1) + + type :: UserData + integer :: neq + integer :: numReac + real(DP) :: rpar + end type UserData + + type(UserData), target :: udata + ! ***************************************************************** ! Explicit declaration of FCVFUN() interface, which is a ! user-supplied function to CVODE. @@ -269,10 +295,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,49 +376,73 @@ end subroutine FCVFUN ! CONFIGURE SOLVER ! ***************************************************************** + + ier = FSUNContext_Create(0, sunctx) + if (ier /= 0) then + write(*,*) 'SUNContext creation failed, ier = ', ier + stop + end if + 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) + n_vector_z = FN_VNew_Serial(neq, sunctx) + if (.not. associated(FN_VGetArrayPointer(n_vector_z))) then + write(stderr,*) 'SUNDIALS_ERROR: FN_VNew_Serial failed' stop end if + ! Initialise the cvode + 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) + cvode_mem = FCVodeCreate(meth, sunctx) + !call FCVMALLOC( , itmeth, iatol, & + ! iout, rout, ipar, rpar, ier ) + if ( .not. c_associated(cvode_mem) ) then + write (stderr,*) 'SUNDIALS_ERROR: FCVodeCreate failed' stop end if - call FCVSETIIN( 'MAX_NSTEPS', maxNumInternalSteps, ier ) + call CVodeInit(cvode_mem, f_rhs, t, n_vector_z, ier) + if (ier /= 0) stop 'CVodeInit failed' + + call CVodeSStolerances(cvode_mem, rtol, atol, ier) + if (ier /= 0) stop 'Tolerance setup failed' + + ier= FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) 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%neq = neq + udata%numReac = numReac + udata%rpar = rpar(1) + call CVodeSetUserData(cvode_mem, c_loc(udata), ier) + 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 ) + LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_NONE, int(lookBack, kind=C_INT), sunctx) + call CVodeSetLinearSolver(cvode_mem, LS, ier) ! SPGMR SOLVER WITH BANDED PRECONDITIONER else if ( solverType == 2 ) then - call FCVSPGMR( 1, 1, lookBack, deltaMain, ier ) - call FCVBPINIT( neq, preconBandUpper, preconBandLower, ier ) + A = FSUNBandMatrix(neq, preconBandUpper, preconBandLower, sunctx) + LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_LEFT, int(lookBack, kind=C_INT), sunctx) + call CVodeSetLinearSolver(cvode_mem, LS, ier) if ( ier /= 0 ) then write (stderr,*) 'SUNDIALS_ERROR: preconditioner returned ier = ', ier ; - call FCVFREE() + call CVodeFree(cvode_mem) stop end if ! DENSE SOLVER else if ( solverType == 3 ) then - call FCVDENSE( neq, ier ) + A = FSUNDenseMatrix(neq, neq, sunctx) + LS = FSUNLinSol_Dense(n_vector_z, A, sunctx) + call CVodeSetLinearSolver(cvode_mem, LS, ier) ! UNEXPECTED SOLVER TYPE else write (stderr,*) 'Error with solverType input, input = ', solverType @@ -404,14 +452,14 @@ end subroutine FCVFUN ! ERROR HANDLING if ( ier /= 0 ) then write (stderr,*) ' SUNDIALS_ERROR: SOLVER returned ier = ', ier - call FCVFREE() + call CVodeFree(cvode_mem) stop end if if ( ier /= 0 ) then write (stderr, 40) ier 40 format (///' SUNDIALS_ERROR: FCVDENSE() returned ier = ', I5) - call FCVFREE() + call CVodeFree(cvode_mem) stop end if @@ -440,10 +488,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, n_vector_z, 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 ) @@ -482,7 +532,7 @@ end subroutine FCVFUN fmt = "(///' SUNDIALS_ERROR: FCVODE() returned ier = ', I5, /, 'Linear Solver returned ier = ', I5) " write (stderr, fmt) ier, iout (15) ! free memory - call FCVFREE() + call CVodeFree(cvode_mem) stop end if @@ -521,7 +571,7 @@ end subroutine FCVFUN ! ***************************************************************** ! deallocate CVODE internal data - call FCVFREE() + call CVodeFree(cvode_mem) deallocate (speciesConcs, z) deallocate (reacDetailedRatesSpecies, prodDetailedRatesSpecies) deallocate (detailedRatesSpeciesName, speciesOfInterest) @@ -626,3 +676,72 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) end subroutine FCVFUN ! ******************************************************************** ! +integer(c_int) function f_rhs(t, y, ydot, user_data) bind(C) + use iso_c_binding + use fnvector_serial_mod + 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 + + ! --- CVODE arguments --- + real(c_double), value, intent(in) :: t + type(N_Vector), intent(in) :: y + type(N_Vector), intent(in) :: ydot + type(c_ptr), value , intent(in) :: user_data + + ! --- local pointers --- + real(DP), pointer :: yv(:), ydv(:) + type(UserData), pointer :: udata + + ! --- local variables --- + integer(NPI) :: numConSpec, np, numReac, i + real(DP), allocatable :: dy(:), z(:), constrainedConcs(:) + + ! access user data + call c_f_pointer(user_data, udata) + numReac = udata%numReac + + ! access N_Vector data + call N_VGetArrayPointer(y, yv) + call N_VGetArrayPointer(ydot, ydv) + + numConSpec = getNumberOfConstrainedSpecies() + np = size(yv) + numConSpec + + allocate (dy(np), z(np), constrainedConcs(numConSpec)) + + ! for each constrained species + do i = 1, numConSpec + 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(yv, constrainedConcs, & + getConstrainedSpecies(), z) + + call resid(numReac, t, z, dy, clhs, clcoeff, crhs, crcoeff) + + call removeConstrainedSpeciesFromProbSpec(dy, & + getConstrainedSpecies(), ydv) + + deallocate(dy, z, constrainedConcs) + + f_rhs = 0 ! success +end function f_rhs \ No newline at end of file From c420cafd71c1f785d633b6363a83350b45044db0 Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Wed, 17 Dec 2025 10:57:39 +0000 Subject: [PATCH 3/9] changes to make sundials work --- src/atchem2.f90 | 242 +++++++++++++++++------------------- tools/install/Makefile.skel | 14 ++- 2 files changed, 127 insertions(+), 129 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 87859d053..504f5ad54 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -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,7 @@ 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 @@ -125,15 +137,9 @@ function dlclose( handle ) bind ( c, name="dlclose" ) type(N_Vector) :: n_vector_z type(c_ptr) :: cvode_mem type(SUNLinearSolver) :: LS - type(SUNMatrix) :: A + type(SUNMatrix), pointer :: A real(c_double) :: t_arr(1) - type :: UserData - integer :: neq - integer :: numReac - real(DP) :: rpar - end type UserData - type(UserData), target :: udata ! ***************************************************************** @@ -162,6 +168,20 @@ 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 :: t + type(N_Vector) :: y + type(N_Vector) :: ydot + type(c_ptr), value :: user_data + end function rhs_fn + + end interface ! ***************************************************************** @@ -386,11 +406,11 @@ end subroutine FCVFUN ipar(1) = neq ipar(2) = numReac - n_vector_z = FN_VNew_Serial(neq, sunctx) - if (.not. associated(FN_VGetArrayPointer(n_vector_z))) then - write(stderr,*) 'SUNDIALS_ERROR: FN_VNew_Serial failed' - stop - end if + !n_vector_z = FN_VNew_Serial(neq, sunctx) + !if (.not. associated(FN_VGetArrayPointer(n_vector_z))) then + ! write(stderr,*) 'SUNDIALS_ERROR: FN_VNew_Serial failed' + ! stop + !end if ! Initialise the cvode @@ -404,64 +424,65 @@ end subroutine FCVFUN stop end if - call CVodeInit(cvode_mem, f_rhs, t, n_vector_z, ier) - if (ier /= 0) stop 'CVodeInit failed' - - call CVodeSStolerances(cvode_mem, rtol, atol, ier) - if (ier /= 0) stop 'Tolerance setup failed' - - ier= FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) - write (*, '(A, I0)') ' setting maxnumsteps ier = ', ier - - ier= FCVodeSetMaxStep(cvode_mem, real(maxStep, kind=C_DOUBLE)) - write (*, '(A, I0)') ' setting maxstep ier = ', ier - write (*,*) - - udata%neq = neq - udata%numReac = numReac - udata%rpar = rpar(1) - call CVodeSetUserData(cvode_mem, c_loc(udata), ier) - if (ier /= 0) stop 'User data setup failed' + !ier = FCVodeInit(cvode_mem, c_funloc(rhs_fn), t, n_vector_z) + !if (ier /= 0) stop 'CVodeInit failed' +! + !ier = FCVodeSStolerances(cvode_mem, rtol, atol) + !if (ier /= 0) stop 'Tolerance setup failed' +! + !ier = FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) + !write (*, '(A, I0)') ' setting maxnumsteps ier = ', 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 - LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_NONE, int(lookBack, kind=C_INT), sunctx) - call CVodeSetLinearSolver(cvode_mem, LS, ier) - ! SPGMR SOLVER WITH BANDED PRECONDITIONER - else if ( solverType == 2 ) then - A = FSUNBandMatrix(neq, preconBandUpper, preconBandLower, sunctx) - LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_LEFT, int(lookBack, kind=C_INT), sunctx) - call CVodeSetLinearSolver(cvode_mem, LS, ier) - if ( ier /= 0 ) then - write (stderr,*) 'SUNDIALS_ERROR: preconditioner returned ier = ', ier ; - call CVodeFree(cvode_mem) - stop - end if - ! DENSE SOLVER - else if ( solverType == 3 ) then - A = FSUNDenseMatrix(neq, neq, sunctx) - LS = FSUNLinSol_Dense(n_vector_z, A, sunctx) - call CVodeSetLinearSolver(cvode_mem, LS, ier) - ! 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 CVodeFree(cvode_mem) - stop - end if - - if ( ier /= 0 ) then - write (stderr, 40) ier - 40 format (///' SUNDIALS_ERROR: FCVDENSE() returned ier = ', I5) - call CVodeFree(cvode_mem) - stop - end if + !if ( solverType == 1 ) then + ! A => null() + ! LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_NONE, int(lookBack, kind=C_INT), sunctx) + ! ier = FCVodeSetLinearSolver(cvode_mem, LS, A) + ! ! SPGMR SOLVER WITH BANDED PRECONDITIONER + !else if ( solverType == 2 ) then + ! A = FSUNBandMatrix(neq, preconBandUpper, preconBandLower, sunctx) + ! LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_LEFT, int(lookBack, kind=C_INT), sunctx) + ! ier = FCVodeSetLinearSolver(cvode_mem, LS, A) + ! if ( ier /= 0 ) then + ! write (stderr,*) 'SUNDIALS_ERROR: preconditioner returned ier = ', ier ; + ! call FCVodeFree(cvode_mem) + ! stop + ! end if + ! ! DENSE SOLVER + !else if ( solverType == 3 ) then + ! A = FSUNDenseMatrix(neq, neq, sunctx) + ! LS = FSUNLinSol_Dense(n_vector_z, A, sunctx) + ! ier = FCVodeSetLinearSolver(cvode_mem, LS, A) + ! ! 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 FCVodeFree(cvode_mem) + ! stop + !end if +! + !if ( ier /= 0 ) then + ! write (stderr, 40) ier + ! 40 format (///' SUNDIALS_ERROR: FCVDENSE() returned ier = ', I5) + ! call FCVodeFree(cvode_mem) + ! stop + !end if ! ***************************************************************** ! RUN MODEL @@ -532,7 +553,7 @@ end subroutine FCVFUN fmt = "(///' SUNDIALS_ERROR: FCVODE() returned ier = ', I5, /, 'Linear Solver returned ier = ', I5) " write (stderr, fmt) ier, iout (15) ! free memory - call CVodeFree(cvode_mem) + call FCVodeFree(cvode_mem) stop end if @@ -571,7 +592,7 @@ end subroutine FCVFUN ! ***************************************************************** ! deallocate CVODE internal data - call CVodeFree(cvode_mem) + call FCVodeFree(cvode_mem) deallocate (speciesConcs, z) deallocate (reacDetailedRatesSpecies, prodDetailedRatesSpecies) deallocate (detailedRatesSpeciesName, speciesOfInterest) @@ -676,72 +697,39 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) end subroutine FCVFUN ! ******************************************************************** ! -integer(c_int) function f_rhs(t, y, ydot, user_data) bind(C) - use iso_c_binding - use fnvector_serial_mod +integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) + use, intrinsic :: iso_c_binding 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 - - ! --- CVODE arguments --- - real(c_double), value, intent(in) :: t - type(N_Vector), intent(in) :: y - type(N_Vector), intent(in) :: ydot - type(c_ptr), value , intent(in) :: user_data - - ! --- local pointers --- - real(DP), pointer :: yv(:), ydv(:) - type(UserData), pointer :: udata - - ! --- local variables --- - integer(NPI) :: numConSpec, np, numReac, i - real(DP), allocatable :: dy(:), z(:), constrainedConcs(:) + use fsundials_core_mod + use fnvector_serial_mod + use cvode_rhs_mod, only : UserData - ! access user data - call c_f_pointer(user_data, udata) - numReac = udata%numReac + implicit none - ! access N_Vector data - call N_VGetArrayPointer(y, yv) - call N_VGetArrayPointer(ydot, ydv) + real(c_double), value :: t + type(N_Vector) :: y + type(N_Vector) :: ydot + type(c_ptr), value :: user_data - numConSpec = getNumberOfConstrainedSpecies() - np = size(yv) + numConSpec + ! Local pointers to vector data + real(c_double), pointer :: ydata(:) + real(c_double), pointer :: ydotdata(:) - allocate (dy(np), z(np), constrainedConcs(numConSpec)) - - ! for each constrained species - do i = 1, numConSpec - if (i <= numberOfVariableConstrainedSpecies) then - call getVariableConstrainedSpeciesConcentrationAtT(t, i, constrainedConcs(i)) - else - constrainedConcs(i) = dataFixedY(i - numberOfVariableConstrainedSpecies) - end if - end do + type(UserData), pointer :: ud + integer(kind=NPI) :: ier - call setConstrainedConcs(constrainedConcs) + 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))] - call addConstrainedSpeciesToProbSpec(yv, constrainedConcs, & - getConstrainedSpecies(), z) - call resid(numReac, t, z, dy, clhs, clcoeff, crhs, crcoeff) + ! Get raw data arrays from N_Vector + ydata => FN_VGetArrayPointer(y) + ydotdata => FN_VGetArrayPointer(ydot) - call removeConstrainedSpeciesFromProbSpec(dy, & - getConstrainedSpecies(), ydv) - deallocate(dy, z, constrainedConcs) + call FCVFUN( t, ydata, ydotdata, ipar_f, ud%rpar, ier ) - f_rhs = 0 ! success -end function f_rhs \ No newline at end of file + rhs_fn = ier ! 0 = success +end function rhs_fn \ No newline at end of file diff --git a/tools/install/Makefile.skel b/tools/install/Makefile.skel index f72b75350..13c3ff649 100644 --- a/tools/install/Makefile.skel +++ b/tools/install/Makefile.skel @@ -22,7 +22,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 -OPENLIBMDIR = openlibm/lib +OPENLIBMDIR = openlibm FRUITDIR = fruit_3.4.3 # =========================================================================== # @@ -82,7 +82,17 @@ endif # set the CVODE and openlibm compilation flags LDFLAGS = -L$(CVODELIBDIR) -L$(OPENLIBMDIR) -Wl,$(RPATH_OPTION),/usr/lib/:$(CVODELIBDIR):$(OPENLIBMDIR) \ - -lopenlibm -lsundials_fcvode_mod -lsundials_cvode -lsundials_fnvecserial_mod -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 From cb6386c695f033d28168c228aabefd3c3a977131 Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Wed, 17 Dec 2025 16:22:22 +0000 Subject: [PATCH 4/9] changed output method and parameters --- src/atchem2.f90 | 219 +++++++++++++++++++++------------------- src/outputFunctions.f90 | 113 ++++++++++++++------- 2 files changed, 195 insertions(+), 137 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 504f5ad54..6e6136993 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -62,6 +62,7 @@ PROGRAM ATCHEM2 use fsunlinsol_dense_mod use fsunmatrix_dense_mod use fsunmatrix_band_mod + use fsunlinsol_band_mod implicit none ! interface to linux API @@ -98,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 @@ -133,11 +134,11 @@ function dlclose( handle ) bind ( c, name="dlclose" ) !sundials declarations - type(c_ptr) :: sunctx - type(N_Vector) :: n_vector_z - type(c_ptr) :: cvode_mem - type(SUNLinearSolver) :: LS - type(SUNMatrix), pointer :: A + 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 @@ -147,7 +148,7 @@ function dlclose( handle ) bind ( c, name="dlclose" ) ! user-supplied function to CVODE. interface - subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) + subroutine FCVFUN(t, y, ydot, ipar, rpar, ier ) use types_mod use species_mod use constraints_mod @@ -155,17 +156,12 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) 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(*) + real(kind=DP), intent(in) :: t + real(kind=DP), intent(inout) :: y(*) + real(kind=DP), intent(inout) :: 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) @@ -175,13 +171,12 @@ integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) implicit none - real(c_double), value :: t - type(N_Vector) :: y - type(N_Vector) :: ydot - type(c_ptr), value :: 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 end function rhs_fn - end interface ! ***************************************************************** @@ -191,9 +186,7 @@ end function rhs_fn 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' @@ -397,12 +390,14 @@ end function rhs_fn ! ***************************************************************** - ier = FSUNContext_Create(0, sunctx) - if (ier /= 0) then - write(*,*) 'SUNContext creation failed, ier = ', ier - stop - end if + !ier = FSUNContext_Create(0, sunctx) + !if (ier /= 0) then + ! write(*,*) 'SUNContext creation failed, ier = ', ier + ! stop + !end if + ! create the SUNDIALS context + ier = FSUNContext_Create(SUN_COMM_NULL, ctx) ipar(1) = neq ipar(2) = numReac @@ -412,71 +407,91 @@ end function rhs_fn ! stop !end if + ! 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 + ! Initialise the cvode write (*, '(A30, 1P e15.3) ') ' t0 = ', t write (*,*) - cvode_mem = FCVodeCreate(meth, sunctx) - !call FCVMALLOC( , itmeth, iatol, & - ! iout, rout, ipar, rpar, ier ) - if ( .not. c_associated(cvode_mem) ) then - write (stderr,*) 'SUNDIALS_ERROR: FCVodeCreate failed' - stop - end if + !cvode_mem = FCVodeCreate(meth, sunctx) + !!call FCVMALLOC( , itmeth, iatol, & + !! rout, ipar, rpar, ier ) + !if ( .not. c_associated(cvode_mem) ) then + ! write (stderr,*) 'SUNDIALS_ERROR: FCVodeCreate failed' + ! stop + !end if + + ! 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, n_vector_z) !if (ier /= 0) stop 'CVodeInit failed' + + 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) stop 'Tolerance setup failed' -! + + 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)) !write (*, '(A, I0)') ' setting maxnumsteps ier = ', 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' + + 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 - ! A => null() - ! LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_NONE, int(lookBack, kind=C_INT), sunctx) - ! ier = FCVodeSetLinearSolver(cvode_mem, LS, A) - ! ! SPGMR SOLVER WITH BANDED PRECONDITIONER - !else if ( solverType == 2 ) then - ! A = FSUNBandMatrix(neq, preconBandUpper, preconBandLower, sunctx) - ! LS = FSUNLinSol_SPGMR(n_vector_z, SUN_PREC_LEFT, int(lookBack, kind=C_INT), sunctx) - ! ier = FCVodeSetLinearSolver(cvode_mem, LS, A) - ! if ( ier /= 0 ) then - ! write (stderr,*) 'SUNDIALS_ERROR: preconditioner returned ier = ', ier ; - ! call FCVodeFree(cvode_mem) - ! stop - ! end if - ! ! DENSE SOLVER - !else if ( solverType == 3 ) then - ! A = FSUNDenseMatrix(neq, neq, sunctx) - ! LS = FSUNLinSol_Dense(n_vector_z, A, sunctx) - ! ier = FCVodeSetLinearSolver(cvode_mem, LS, A) - ! ! 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 FCVodeFree(cvode_mem) - ! stop - !end if -! + if ( solverType == 1 ) then + sunmat_A => null() + 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 + ! ERROR HANDLING + + !if ( ier /= 0 ) then ! write (stderr, 40) ier ! 40 format (///' SUNDIALS_ERROR: FCVDENSE() returned ier = ', I5) @@ -484,6 +499,15 @@ end function rhs_fn ! stop !end if + + ! Attach the matrix and linear solver + ier = FCVodeSetLinearSolver(cvode_mem, sunls, sunmat_A); + if ( ier /= 0 ) then + write (stderr,*) ' SUNDIALS_ERROR: FCVodeSetLinearSolver returned ier = ', ier + call FCVodeFree(cvode_mem) + stop + end if + ! ***************************************************************** ! RUN MODEL ! ***************************************************************** @@ -510,7 +534,7 @@ end function rhs_fn ! Get concentrations for unconstrained species t_arr(1) =t - ier = FCVode(cvode_mem, tout, n_vector_z, t_arr, itask) + ier = FCVode(cvode_mem, tout, sunvec_u, t_arr, itask) if ( ier /= 0 ) then write (*, '(A, I0)') ' ier POST FCVODE()= ', ier end if @@ -541,9 +565,6 @@ end function rhs_fn 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 ) @@ -551,7 +572,6 @@ end function rhs_fn ! 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 FCVodeFree(cvode_mem) stop @@ -566,21 +586,13 @@ end function rhs_fn ! 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 @@ -644,7 +656,7 @@ END PROGRAM ATCHEM2 ! -------------------------------------------------------- ! ! Fortran routine for right-hand side function. -subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) +subroutine FCVFUN(t, y, ydot, ipar, rpar, ier ) use types_mod use constraints_mod, only : getNumberOfConstrainedSpecies, numberOfVariableConstrainedSpecies, dataFixedY, & getConstrainedSpecies, setConstrainedConcs @@ -655,8 +667,9 @@ subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) use solver_functions_mod, only : resid implicit none - real(kind=DP), intent(in) :: t, y(*) - real(kind=DP), intent(out) :: ydot(*) + real(kind=DP), intent(in) :: t + real(kind=DP), intent(inout) :: y(*) + real(kind=DP), intent(inout) :: ydot(*) integer(kind=NPI), intent(in) :: ipar(*) real(kind=DP), intent(in) :: rpar(*) integer(kind=NPI), intent(out) :: ier @@ -703,13 +716,14 @@ integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) use fsundials_core_mod use fnvector_serial_mod use cvode_rhs_mod, only : UserData + use types_mod implicit none - real(c_double), value :: t - type(N_Vector) :: y - type(N_Vector) :: ydot - type(c_ptr), value :: 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(:) @@ -723,7 +737,6 @@ integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) 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) @@ -732,4 +745,4 @@ integer(c_int) function rhs_fn(t, y, ydot, user_data) bind(C) call FCVFUN( t, ydata, ydotdata, ipar_f, ud%rpar, ier ) rhs_fn = ier ! 0 = success -end function rhs_fn \ No newline at end of file +end function rhs_fn diff --git a/src/outputFunctions.f90 b/src/outputFunctions.f90 index b96930869..e61dc288d 100644 --- a/src/outputFunctions.f90 +++ b/src/outputFunctions.f90 @@ -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 c490810d3ea18123b33d3dcce8afd3d202ff930a Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Wed, 17 Dec 2025 16:42:50 +0000 Subject: [PATCH 5/9] fixed the 3 types of linear solver --- src/atchem2.f90 | 46 +++------------------------------------------- 1 file changed, 3 insertions(+), 43 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 6e6136993..fd9adf563 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -62,7 +62,7 @@ PROGRAM ATCHEM2 use fsunlinsol_dense_mod use fsunmatrix_dense_mod use fsunmatrix_band_mod - use fsunlinsol_band_mod + use fsunlinsol_band_mod implicit none ! interface to linux API @@ -389,24 +389,11 @@ end function rhs_fn ! CONFIGURE SOLVER ! ***************************************************************** - - !ier = FSUNContext_Create(0, sunctx) - !if (ier /= 0) then - ! write(*,*) 'SUNContext creation failed, ier = ', ier - ! stop - !end if - ! create the SUNDIALS context ier = FSUNContext_Create(SUN_COMM_NULL, ctx) ipar(1) = neq ipar(2) = numReac - !n_vector_z = FN_VNew_Serial(neq, sunctx) - !if (.not. associated(FN_VGetArrayPointer(n_vector_z))) then - ! write(stderr,*) 'SUNDIALS_ERROR: FN_VNew_Serial failed' - ! stop - !end if - ! create SUNDIALS N_Vector sunvec_u => FN_VMake_Serial(neq, z, ctx) if (.not. associated(sunvec_u)) then @@ -414,42 +401,26 @@ end function rhs_fn stop 1 end if - ! Initialise the cvode + write (*, '(A30, 1P e15.3) ') ' t0 = ', t write (*,*) - !cvode_mem = FCVodeCreate(meth, sunctx) - !!call FCVMALLOC( , itmeth, iatol, & - !! rout, ipar, rpar, ier ) - !if ( .not. c_associated(cvode_mem) ) then - ! write (stderr,*) 'SUNDIALS_ERROR: FCVodeCreate failed' - ! stop - !end if ! 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, n_vector_z) - !if (ier /= 0) stop 'CVodeInit failed' 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) stop 'Tolerance setup failed' 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)) - !write (*, '(A, I0)') ' setting maxnumsteps ier = ', ier ier = FCVodeSetMaxNumSteps(cvode_mem, int(maxNumInternalSteps, kind=C_LONG)) if (ier /= 0) then @@ -467,12 +438,9 @@ end function rhs_fn 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 - sunmat_A => null() sunls => FSUNLinSol_SPGMR(sunvec_u, SUN_PREC_NONE, int(lookBack, kind=C_INT), ctx) ! SPGMR SOLVER WITH BANDED PRECONDITIONER else if ( solverType == 2 ) then @@ -489,19 +457,11 @@ end function rhs_fn write (stderr,*) 'Available options are 1, 2, 3.' stop end if - ! ERROR HANDLING - - - !if ( ier /= 0 ) then - ! write (stderr, 40) ier - ! 40 format (///' SUNDIALS_ERROR: FCVDENSE() returned ier = ', I5) - ! call FCVodeFree(cvode_mem) - ! 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) From cada090edd7967a259cc847415d0e72e5ac7f554 Mon Sep 17 00:00:00 2001 From: Neil Butcher Date: Wed, 17 Dec 2025 17:01:49 +0000 Subject: [PATCH 6/9] revert change to docker file --- docker/install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) mode change 100644 => 100755 docker/install.sh diff --git a/docker/install.sh b/docker/install.sh old mode 100644 new mode 100755 index 391806f91..1c8d0ef7a --- a/docker/install.sh +++ b/docker/install.sh @@ -48,7 +48,7 @@ cd atchem # Change atchem dependancy paths and create Makefile from skeleton sed 's,cvode/lib,/atchem-lib/cvode/lib,g' tools/install/Makefile.skel > ./Makefile -sed -i 's,openlibm/lib,/atchem-lib/openlibm/,g' ./Makefile +sed -i 's,openlibm-0.8.1,/atchem-lib/openlibm-0.8.1,g' ./Makefile # Fix python command to match installed version sed -i "s/python/python3/g" ./build/build_atchem2.sh From 041c9941d980f2bfe1227a55f7a349164a06ef79 Mon Sep 17 00:00:00 2001 From: Neil Butcher Date: Thu, 18 Dec 2025 08:46:36 +0000 Subject: [PATCH 7/9] We needed to add cvode object files --- tools/install/Makefile.skel | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tools/install/Makefile.skel b/tools/install/Makefile.skel index 13c3ff649..70da69688 100644 --- a/tools/install/Makefile.skel +++ b/tools/install/Makefile.skel @@ -22,6 +22,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 @@ -116,7 +117,7 @@ SRCS = $(CORE_SRCS) $(SRC)/atchem2.f90 # the executable is rebuilt every time a fortran source file 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 dcb71626b0395b192101f8c4306f132ad3842a42 Mon Sep 17 00:00:00 2001 From: Neil Butcher Date: Thu, 18 Dec 2025 08:54:55 +0000 Subject: [PATCH 8/9] changing the readme and docker file script to try to inform people about the new setting for fortran object files --- README.md | 2 +- docker/install.sh | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index bd4d311cb..89787ea31 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,7 @@ AtChem2 requires a **Fortran** compiler (GNU `gfortran` or Intel `ifort`), the * The latest stable version of AtChem2 can be downloaded from the [Releases page](https://github.com/AtChem/AtChem2/releases), and is associated with a **doi number** for referencing in publications. -After installing the required dependencies using the scripts in the `tools/install/` directory, copy the file `tools/install/Makefile.skel` to the _Main Directory_ and rename it `Makefile`. Set the variables `CVODELIBDIR`, `OPENLIBMDIR` and `FRUITDIR` in the `Makefile` to the full paths of CVODE, openlibm and (if installed) FRUIT. +After installing the required dependencies using the scripts in the `tools/install/` directory, copy the file `tools/install/Makefile.skel` to the _Main Directory_ and rename it `Makefile`. Set the variables `CVODELIBDIR`, `CVODEOBJDIR`, `OPENLIBMDIR` and `FRUITDIR` in the `Makefile` to the full paths of CVODE (libraries and fortran object files), openlibm and (if installed) FRUIT. Alternatively, and _optionally_, AtChem2 can be run as a [Docker](https://www.docker.com/) container. Currently the containerized version is available only for AtChem2 v1.2.2, thanks to the [Uni York group](https://github.com/wacl-york/AtChem2/pkgs/container/atchem2). diff --git a/docker/install.sh b/docker/install.sh index 1c8d0ef7a..78fb9c397 100755 --- a/docker/install.sh +++ b/docker/install.sh @@ -47,7 +47,8 @@ cd atchem ./tools/install/install_openlibm.sh /atchem-lib/ # Change atchem dependancy paths and create Makefile from skeleton -sed 's,cvode/lib,/atchem-lib/cvode/lib,g' tools/install/Makefile.skel > ./Makefile +sed 's,cvode/lib,/atchem-lib/cvode/lib64,g' tools/install/Makefile.skel > ./Makefile +sed -i 's,cvode/fortran,/atchem-lib/cvode/fortran,g' ./Makefile sed -i 's,openlibm-0.8.1,/atchem-lib/openlibm-0.8.1,g' ./Makefile # Fix python command to match installed version From 977dc15f45a8b4763bdb2250c49b879fc225e2fc Mon Sep 17 00:00:00 2001 From: "Neil Butcher (Advanced Research Computing)" Date: Thu, 18 Dec 2025 09:15:46 +0000 Subject: [PATCH 9/9] reverted some of the changes to the old fortran method --- src/atchem2.f90 | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/atchem2.f90 b/src/atchem2.f90 index fd9adf563..f26d6f502 100644 --- a/src/atchem2.f90 +++ b/src/atchem2.f90 @@ -148,7 +148,7 @@ function dlclose( handle ) bind ( c, name="dlclose" ) ! user-supplied function to CVODE. interface - subroutine FCVFUN(t, y, ydot, ipar, rpar, ier ) + subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) use types_mod use species_mod use constraints_mod @@ -156,12 +156,17 @@ subroutine FCVFUN(t, y, ydot, ipar, rpar, ier ) use interpolation_functions_mod, only : getVariableConstrainedSpeciesConcentrationAtT use constraint_functions_mod - real(kind=DP), intent(in) :: t - real(kind=DP), intent(inout) :: y(*) - real(kind=DP), intent(inout) :: ydot(*) + ! 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) @@ -616,7 +621,7 @@ END PROGRAM ATCHEM2 ! -------------------------------------------------------- ! ! Fortran routine for right-hand side function. -subroutine FCVFUN(t, y, ydot, ipar, rpar, ier ) +subroutine FCVFUN( t, y, ydot, ipar, rpar, ier ) use types_mod use constraints_mod, only : getNumberOfConstrainedSpecies, numberOfVariableConstrainedSpecies, dataFixedY, & getConstrainedSpecies, setConstrainedConcs @@ -627,9 +632,8 @@ subroutine FCVFUN(t, y, ydot, ipar, rpar, ier ) use solver_functions_mod, only : resid implicit none - real(kind=DP), intent(in) :: t - real(kind=DP), intent(inout) :: y(*) - real(kind=DP), intent(inout) :: ydot(*) + 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