1! ------------------------------------------------------------------ 2! Programmer(s): Daniel R. Reynolds @ SMU 3! ------------------------------------------------------------------ 4! SUNDIALS Copyright Start 5! Copyright (c) 2002-2020, Lawrence Livermore National Security 6! and Southern Methodist University. 7! All rights reserved. 8! 9! See the top-level LICENSE and NOTICE files for details. 10! 11! SPDX-License-Identifier: BSD-3-Clause 12! SUNDIALS Copyright End 13! ------------------------------------------------------------------ 14! Example problem: 15! 16! The following is a simple example problem with analytical 17! solution, 18! dy/dt = lamda*y + 1/(1+t^2) - lamda*atan(t) 19! 20! The stiffness of the problem is directly proportional to the 21! value of "lamda", which is specified through an input file. The 22! real part of lamda should be negative to result in a well-posed 23! ODE; for lambdas with magnitude larger than 100 the problem 24! becomes quite stiff. 25! 26! Here we choose lamda = -0.1 + 10 i. 27! 28! This program solves the problem with the ERK method and 29! a user-supplied complex vector module. 30! Output is printed every 1.0 units of time (10 total). 31! Run statistics (optional outputs) are printed at the end. 32! ------------------------------------------------------------------ 33 34module ode_mod 35 36 !======= Inclusions =========== 37 use, intrinsic :: iso_c_binding 38 39 !======= Declarations ========= 40 implicit none 41 integer(c_long), parameter :: neq = 1 42 integer(c_int), parameter :: Nt = 10 43 complex(c_double_complex), parameter :: lambda = (-1d-2, 10.d0) 44 real(c_double), parameter :: T0 = 0.d0 45 real(c_double), parameter :: Tf = 10.d0 46 real(c_double), parameter :: dtmax = 0.01d0 47 real(c_double), parameter :: reltol = 1.d-6 48 real(c_double), parameter :: abstol = 1.d-10 49 50contains 51 52 ! ---------------------------------------------------------------- 53 ! Rhs provides the right hand side function for the ODE: 54 ! dy/dt = f(t,y). 55 ! 56 ! Return values: 57 ! 0 = success, 58 ! 1 = recoverable error, 59 ! -1 = non-recoverable error 60 ! ---------------------------------------------------------------- 61 integer(c_int) function Rhs(tn, sunvec_y, sunvec_f, user_data) & 62 result(ierr) bind(C,name='Rhs') 63 64 !======= Inclusions =========== 65 use, intrinsic :: iso_c_binding 66 use fnvector_complex_mod 67 68 !======= Declarations ========= 69 implicit none 70 71 ! calling variables 72 real(c_double), value :: tn ! current time 73 type(N_Vector) :: sunvec_y ! solution N_Vector 74 type(N_Vector) :: sunvec_f ! rhs N_Vector 75 type(c_ptr) :: user_data ! user-defined data 76 77 ! local variables 78 type(FVec), pointer :: y, f ! ptrs to Fortran vector data 79 80 !======= Internals ============ 81 82 ! extract Fortran vector structures to work with 83 y => FN_VGetFVec(sunvec_y) 84 f => FN_VGetFVec(sunvec_f) 85 86 ! evaluate ODE RHS and return success 87 f%data(1) = lambda*y%data(1) 88 ierr = 0 89 90 end function Rhs 91 92 ! ---------------------------------------------------------------- 93 ! Sol provides the analytical solution to the ODE. 94 ! ---------------------------------------------------------------- 95 complex(c_double_complex) function Sol(tn) result(y) 96 97 !======= Inclusions =========== 98 use, intrinsic :: iso_c_binding 99 use fnvector_complex_mod 100 101 !======= Declarations ========= 102 implicit none 103 104 ! calling variables 105 real(c_double), value :: tn ! current time 106 107 !======= Internals ============ 108 109 ! fill analytical solution for this value of tn, and return 110 y = 2.d0*exp(lambda*tn) 111 return 112 113 end function Sol 114 115end module ode_mod 116 117program main 118 119 !======= Inclusions =========== 120 use, intrinsic :: iso_c_binding 121 122 use farkode_mod ! Fortran interface to the ARKode module 123 use farkode_arkstep_mod ! Fortran interface to the ARKStep module 124 use fnvector_complex_mod ! Custom complex N_Vector 125 use ode_mod ! ODE functions 126 127 !======= Declarations ========= 128 implicit none 129 130 ! local variables 131 integer(c_int) :: ierr, iout 132 real(c_double) :: tcur(1), tout, dTout, yerr, yerrI, yerr2 133 134 type(N_Vector), pointer :: sunvec_y ! sundials vector 135 type(c_ptr) :: arkode_mem ! ARKODE memory 136 137 ! solution vector 138 type(FVec), pointer :: y 139 140 !======= Internals ============ 141 142 ! initial problem output 143 print *, " " 144 print *, "Analytical ODE test problem:" 145 print '(2(a,f5.2),a)', " lambda = (", real(lambda), " , ", imag(lambda), " ) " 146 print '(2(a,es8.1))', " reltol = ",reltol,", abstol = ",abstol 147 148 ! initialize SUNDIALS solution vector 149 sunvec_y => FN_VNew_Complex(neq) 150 if (.not. associated(sunvec_y)) then 151 print *, 'ERROR: sunvec = NULL' 152 stop 1 153 end if 154 y => FN_VGetFVec(sunvec_y) 155 156 ! Set initial conditions into y 157 y%data(1) = Sol(T0) 158 159 ! create ARKStep memory 160 arkode_mem = FARKStepCreate(c_funloc(Rhs), c_null_funptr, T0, sunvec_y) 161 if (.not. c_associated(arkode_mem)) then 162 print *,'ERROR: arkode_mem = NULL' 163 stop 1 164 end if 165 166 ! main time-stepping loop: calls FARKStepEvolve to perform the integration, then 167 ! prints results. Stops when the final time has been reached 168 tcur(1) = T0 169 dTout = (Tf-T0)/Nt 170 tout = T0+dTout 171 yerrI = 0.d0 172 yerr2 = 0.d0 173 print *, " " 174 print *, " t real(u) imag(u) error" 175 print *, " -------------------------------------------" 176 print '(5x,f4.1,2(2x,es9.2),2x,es8.1)', tcur(1), real(y%data(1)), imag(y%data(1)), 0.d0 177 do iout = 1,Nt 178 179 ! call integrator 180 ierr = FARKStepEvolve(arkode_mem, tout, sunvec_y, tcur, ARK_NORMAL) 181 if (ierr /= 0) then 182 write(*,*) 'Error in FARKStepEvolve, ierr = ', ierr, '; halting' 183 stop 1 184 endif 185 186 ! compute/accumulate solution error 187 yerr = abs( y%data(1) - Sol(tcur(1)) ) 188 yerrI = max(yerrI, yerr) 189 yerr2 = yerr2 + yerr**2 190 191 ! print solution statistics 192 print '(5x,f4.1,2(2x,es9.2),2x,es8.1)', tcur(1), real(y%data(1)), imag(y%data(1)), yerr 193 194 ! update output time 195 tout = min(tout + dTout, Tf) 196 197 end do 198 yerr2 = dsqrt( yerr2 / Nt ) 199 print *, " -------------------------------------------" 200 201 ! diagnostics output 202 call ARKStepStats(arkode_mem) 203 print '(2(a,es9.2))', " Error: max = ", yerrI, ", rms = ", yerr2 204 print *, ' ' 205 206 ! clean up 207 call FARKStepFree(arkode_mem) 208 call FN_VDestroy(sunvec_y) 209 210end program main 211 212 213! ---------------------------------------------------------------- 214! ARKStepStats 215! 216! Print ARKODE statstics to standard out 217! ---------------------------------------------------------------- 218subroutine ARKStepStats(arkode_mem) 219 220 !======= Inclusions =========== 221 use iso_c_binding 222 use farkode_mod 223 use farkode_arkstep_mod 224 225 !======= Declarations ========= 226 implicit none 227 228 type(c_ptr), intent(in) :: arkode_mem ! solver memory structure 229 230 integer(c_int) :: ierr ! error flag 231 integer(c_long) :: nsteps(1) ! num steps 232 integer(c_long) :: nst_a(1) ! num steps attempted 233 integer(c_long) :: nfe(1) ! num explicit function evals 234 integer(c_long) :: nfi(1) ! num implicit function evals 235 integer(c_long) :: netfails(1) ! num error test fails 236 237 !======= Internals ============ 238 239 ierr = FARKStepGetNumSteps(arkode_mem, nsteps) 240 ierr = FARKStepGetNumStepAttempts(arkode_mem, nst_a) 241 ierr = FARKStepGetNumRhsEvals(arkode_mem, nfe, nfi) 242 ierr = FARKStepGetNumErrTestFails(arkode_mem, netfails) 243 244 print *, ' ' 245 print *, 'Final Solver Statistics:' 246 print '(4x,2(A,i4),A)' ,'Internal solver steps = ',nsteps(1),', (attempted = ',nst_a(1),')' 247 print '(4x,A,i5)' ,'Total RHS evals = ',nfe(1) 248 print '(4x,A,i5)' ,'Total number of error test failures =',netfails(1) 249 250 return 251 252end subroutine ARKStepStats 253