1! ------------------------------------------------------------------ 2! Programmer(s): David J. Gardner, Cody J. Balos @ LLNL 3! Daniel R. Reynolds @ SMU 4! ------------------------------------------------------------------ 5! SUNDIALS Copyright Start 6! Copyright (c) 2002-2021, Lawrence Livermore National Security 7! and Southern Methodist University. 8! All rights reserved. 9! 10! See the top-level LICENSE and NOTICE files for details. 11! 12! SPDX-License-Identifier: BSD-3-Clause 13! SUNDIALS Copyright End 14! ------------------------------------------------------------------ 15! This test simulates a brusselator problem from chemical kinetics. 16! This is an ODE system with 3 components, Y = [u,v,w], satisfying 17! the equations, 18! 19! du/dt = a - (w+1)*u + v*u^2 20! dv/dt = w*u - v*u^2 21! dw/dt = (b-w)/ep - w*u 22! 23! for t in the interval [0.0, 10.0], with initial conditions 24! Y0 = [3.9, 1.1, 2.8], and parameter values a=1.2, b=2.5, and 25! ep=1.0e-5 26! 27! Here, all three components exhibit a rapid transient change during 28! the first 0.2 time units, followed by a slow and smooth evolution. 29! ------------------------------------------------------------------ 30 31module ode_mod 32 33 !======= Inclusions =========== 34 use, intrinsic :: iso_c_binding 35 36 !======= Declarations ========= 37 implicit none 38 39 ! number of equations 40 integer(c_long), parameter :: neq = 3 41 42 ! ODE parameters 43 double precision, parameter :: a = 1.2d0 44 double precision, parameter :: b = 2.5d0 45 double precision, parameter :: ep = 1.0d-5 46 47contains 48 49 ! ---------------------------------------------------------------- 50 ! RhsFn provides the right hand side function for the 51 ! ODE: dy/dt = f(t,y) 52 ! 53 ! Return values: 54 ! 0 = success, 55 ! 1 = recoverable error, 56 ! -1 = non-recoverable error 57 ! ---------------------------------------------------------------- 58 integer(c_int) function RhsFn(tn, sunvec_y, sunvec_f, user_data) & 59 result(ierr) bind(C,name='RhsFn') 60 61 !======= Inclusions =========== 62 use, intrinsic :: iso_c_binding 63 use fsundials_nvector_mod 64 65 !======= Declarations ========= 66 implicit none 67 68 ! calling variables 69 real(c_double), value :: tn ! current time 70 type(N_Vector) :: sunvec_y ! solution N_Vector 71 type(N_Vector) :: sunvec_f ! rhs N_Vector 72 type(c_ptr), value :: user_data ! user-defined data 73 74 ! pointers to data in SUNDIALS vectors 75 real(c_double), pointer :: yvec(:) 76 real(c_double), pointer :: fvec(:) 77 78 !======= Internals ============ 79 80 ! get data arrays from SUNDIALS vectors 81 yvec => FN_VGetArrayPointer(sunvec_y) 82 fvec => FN_VGetArrayPointer(sunvec_f) 83 84 ! fill RHS vector 85 fvec(1) = a - (yvec(3) + 1.0d0) * yvec(1) + yvec(2) * yvec(1) * yvec(1) 86 fvec(2) = yvec(3) * yvec(1) - yvec(2) * yvec(1) * yvec(1) 87 fvec(3) = (b-yvec(3))/ep - yvec(3) * yvec(1) 88 89 ! return success 90 ierr = 0 91 return 92 93 end function RhsFn 94 95 ! ---------------------------------------------------------------- 96 ! JacFn: The Jacobian of the ODE hand side function J = df/dy 97 ! 98 ! Return values: 99 ! 0 = success, 100 ! 1 = recoverable error, 101 ! -1 = non-recoverable error 102 ! ---------------------------------------------------------------- 103 integer(c_int) function JacFn(tn, sunvec_y, sunvec_f, sunmat_J, & 104 user_data, tmp1, tmp2, tmp3) & 105 result(ierr) bind(C,name='JacFn') 106 107 !======= Inclusions =========== 108 use, intrinsic :: iso_c_binding 109 use fsundials_nvector_mod 110 use fsunmatrix_dense_mod 111 use fsundials_matrix_mod 112 113 !======= Declarations ========= 114 implicit none 115 116 ! calling variables 117 real(c_double), value :: tn ! current time 118 type(N_Vector) :: sunvec_y ! current solution N_Vector 119 type(N_Vector) :: sunvec_f ! current rhs N_Vector 120 type(SUNMatrix) :: sunmat_J ! Jacobian SUNMatrix 121 type(c_ptr), value :: user_data ! user-defined data 122 type(N_Vector) :: tmp1, tmp2, tmp3 ! workspace N_Vectors 123 124 ! pointers to data in SUNDIALS vectors 125 real(c_double), pointer :: yvec(:) 126 127 ! pointer to data in SUNDIALS matrix 128 real(c_double), pointer :: Jmat(:) 129 130 !======= Internals ============ 131 132 ! get data array from SUNDIALS vector 133 yvec => FN_VGetArrayPointer(sunvec_y) 134 135 ! get data arrays from SUNDIALS vectors 136 Jmat => FSUNDenseMatrix_Data(sunmat_J) 137 138 ! fill Jacobian matrix 139 Jmat = [-(yvec(3)+1.0d0) + 2.0d0*yvec(1)*yvec(2),& 140 yvec(3) - 2.0d0*yvec(1)*yvec(2), -yvec(3),& 141 yvec(1)*yvec(1), -yvec(1)*yvec(1), 0.0d0,& 142 -yvec(1), yvec(1), -1.0d0/ep - yvec(1)] 143 144 ! return success 145 ierr = 0 146 return 147 148 end function JacFn 149 150end module ode_mod 151 152 153program main 154 155 !======= Inclusions =========== 156 use, intrinsic :: iso_c_binding 157 158 use fcvode_mod ! Fortran interface to CVODE 159 use fnvector_serial_mod ! Fortran interface to serial N_Vector 160 use fsunmatrix_dense_mod ! Fortran interface to dense SUNMatrix 161 use fsunlinsol_dense_mod ! Fortran interface to dense SUNLinearSolver 162 use fsundials_linearsolver_mod ! Fortran interface to generic SUNLinearSolver 163 use fsundials_matrix_mod ! Fortran interface to generic SUNMatrix 164 use fsundials_nvector_mod ! Fortran interface to generic N_Vector 165 use ode_mod ! ODE functions 166 167 !======= Declarations ========= 168 implicit none 169 170 ! local variables 171 real(c_double) :: tstart ! initial time 172 real(c_double) :: tend ! final time 173 real(c_double) :: rtol, atol ! relative and absolute tolerance 174 real(c_double) :: dtout ! output time interval 175 real(c_double) :: tout ! output time 176 real(c_double) :: tcur(1) ! current time 177 178 integer(c_int) :: ierr ! error flag from C functions 179 integer(c_int) :: nout ! number of outputs 180 181 integer :: outstep ! output loop counter 182 183 type(c_ptr) :: cvode_mem ! CVODE memory 184 type(N_Vector), pointer :: sunvec_y ! sundials vector 185 type(SUNMatrix), pointer :: sunmat_A ! sundials matrix 186 type(SUNLinearSolver), pointer :: sunlinsol_LS ! sundials linear solver 187 188 ! solution vector, neq is set in the ode_mod module 189 real(c_double) :: yvec(neq) 190 191 !======= Internals ============ 192 193 ! initialize ODE 194 tstart = 0.0d0 195 tend = 10.0d0 196 tcur = tstart 197 tout = tstart 198 dtout = (tend-tstart)/10.d0 199 nout = ceiling(tend/dtout) 200 201 ! initialize solution vector 202 yvec(1) = 3.9d0 203 yvec(2) = 1.1d0 204 yvec(3) = 2.8d0 205 206 ! create SUNDIALS N_Vector 207 sunvec_y => FN_VMake_Serial(neq, yvec) 208 if (.not. associated(sunvec_y)) then 209 print *, 'ERROR: sunvec = NULL' 210 stop 1 211 end if 212 213 ! create a dense matrix 214 sunmat_A => FSUNDenseMatrix(neq, neq) 215 if (.not. associated(sunmat_A)) then 216 print *, 'ERROR: sunmat = NULL' 217 stop 1 218 end if 219 220 ! create a dense linear solver 221 sunlinsol_LS => FSUNDenseLinearSolver(sunvec_y, sunmat_A) 222 if (.not. associated(sunlinsol_LS)) then 223 print *, 'ERROR: sunlinsol = NULL' 224 stop 1 225 end if 226 227 ! create CVode memory 228 cvode_mem = FCVodeCreate(CV_BDF) 229 if (.not. c_associated(cvode_mem)) then 230 print *, 'ERROR: cvode_mem = NULL' 231 stop 1 232 end if 233 234 ! initialize CVode 235 ierr = FCVodeInit(cvode_mem, c_funloc(RhsFn), tstart, sunvec_y) 236 if (ierr /= 0) then 237 print *, 'Error in FCVodeInit, ierr = ', ierr, '; halting' 238 stop 1 239 end if 240 241 ! set relative and absolute tolerances 242 rtol = 1.0d-5 243 atol = 1.0d-10 244 245 ierr = FCVodeSStolerances(cvode_mem, rtol, atol) 246 if (ierr /= 0) then 247 print *, 'Error in FCVodeSStolerances, ierr = ', ierr, '; halting' 248 stop 1 249 end if 250 251 ! attach linear solver 252 ierr = FCVodeSetLinearSolver(cvode_mem, sunlinsol_LS, sunmat_A); 253 if (ierr /= 0) then 254 print *, 'Error in FCVodeSetLinearSolver, ierr = ', ierr, '; halting' 255 stop 1 256 end if 257 258 ! set Jacobian routine 259 ierr = FCVodeSetJacFn(cvode_mem, c_funloc(JacFn)) 260 if (ierr /= 0) then 261 print *, 'Error in FCVodeSetJacFn, ierr = ', ierr, '; halting' 262 stop 1 263 end if 264 265 ! start time stepping 266 print *, ' ' 267 print *, 'Finished initialization, starting time steps' 268 print *, ' ' 269 print *, ' t u v w' 270 print *, '----------------------------------------------------' 271 print '(1x,4(es12.5,1x))', tcur, yvec(1), yvec(2), yvec(3) 272 do outstep = 1,nout 273 274 ! call CVode 275 tout = min(tout + dtout, tend) 276 ierr = FCVode(cvode_mem, tout, sunvec_y, tcur, CV_NORMAL) 277 if (ierr /= 0) then 278 print *, 'Error in FCVODE, ierr = ', ierr, '; halting' 279 stop 1 280 endif 281 282 ! output current solution 283 print '(1x,4(es12.5,1x))', tcur, yvec(1), yvec(2), yvec(3) 284 285 enddo 286 287 ! diagnostics output 288 call CVodeStats(cvode_mem) 289 290 ! clean up 291 call FCVodeFree(cvode_mem) 292 ierr = FSUNLinSolFree(sunlinsol_LS) 293 call FSUNMatDestroy(sunmat_A) 294 call FN_VDestroy(sunvec_y) 295 296end program Main 297 298 299! ---------------------------------------------------------------- 300! CVodeStats 301! 302! Print CVODE statstics to standard out 303! ---------------------------------------------------------------- 304subroutine CVodeStats(cvode_mem) 305 306 !======= Inclusions =========== 307 use iso_c_binding 308 use fcvode_mod 309 310 !======= Declarations ========= 311 implicit none 312 313 type(c_ptr), intent(in) :: cvode_mem ! solver memory structure 314 315 integer(c_int) :: ierr ! error flag 316 317 integer(c_long) :: nsteps(1) ! num steps 318 integer(c_long) :: nfevals(1) ! num function evals 319 integer(c_long) :: nlinsetups(1) ! num linear solver setups 320 integer(c_long) :: netfails(1) ! num error test fails 321 322 integer(c_int) :: qlast(1) ! method order in last step 323 integer(c_int) :: qcur(1) ! method order for next step 324 325 real(c_double) :: hinused(1) ! initial step size 326 real(c_double) :: hlast(1) ! last step size 327 real(c_double) :: hcur(1) ! step size for next step 328 real(c_double) :: tcur(1) ! internal time reached 329 330 integer(c_long) :: nniters(1) ! nonlinear solver iterations 331 integer(c_long) :: nncfails(1) ! nonlinear solver fails 332 333 integer(c_long) :: njevals(1) ! num Jacobian evaluations 334 335 !======= Internals ============ 336 337 ! general solver statistics 338 ierr = FCVodeGetIntegratorStats(cvode_mem, nsteps, nfevals, nlinsetups, & 339 netfails, qlast, qcur, hinused, hlast, hcur, tcur) 340 if (ierr /= 0) then 341 print *, 'Error in FCVodeGetIntegratorStats, ierr = ', ierr, '; halting' 342 stop 1 343 end if 344 345 ! nonlinear solver statistics 346 ierr = FCVodeGetNonlinSolvStats(cvode_mem, nniters, nncfails) 347 if (ierr /= 0) then 348 print *, 'Error in FCVodeGetNonlinSolvStats, ierr = ', ierr, '; halting' 349 stop 1 350 end if 351 352 ! nonlinear solver statistics 353 ierr = FCVodeGetNumJacEvals(cvode_mem, njevals) 354 if (ierr /= 0) then 355 print *, 'Error in FCVodeGetNumJacEvals, ierr = ', ierr, '; halting' 356 stop 1 357 end if 358 359 print *, ' ' 360 print *, ' General Solver Stats:' 361 print '(4x,A,i9)' ,'Total internal steps taken =',nsteps 362 print '(4x,A,i9)' ,'Total rhs function calls =',nfevals 363 print '(4x,A,i9)' ,'Num lin solver setup calls =',nlinsetups 364 print '(4x,A,i9)' ,'Num error test failures =',netfails 365 print '(4x,A,i9)' ,'Last method order =',qlast 366 print '(4x,A,i9)' ,'Next method order =',qcur 367 print '(4x,A,es12.5)','First internal step size =',hinused 368 print '(4x,A,es12.5)','Last internal step size =',hlast 369 print '(4x,A,es12.5)','Next internal step size =',hcur 370 print '(4x,A,es12.5)','Current internal time =',tcur 371 print '(4x,A,i9)' ,'Num nonlinear solver iters =',nniters 372 print '(4x,A,i9)' ,'Num nonlinear solver fails =',nncfails 373 print '(4x,A,i9)' ,'Num Jacobian evaluations =',njevals 374 print *, ' ' 375 376 return 377 378end subroutine CVodeStats 379