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