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 diff --git a/src/atchem2.f90 b/src/atchem2.f90 index 496921665..f26d6f502 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,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 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. diff --git a/tools/install/Makefile.skel b/tools/install/Makefile.skel index 40a2e1ff6..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 @@ -82,7 +83,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 @@ -106,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 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"