Skip to content

Commit

Permalink
Allow time integrator initialise with a method name.
Browse files Browse the repository at this point in the history
  • Loading branch information
semi-h committed Jul 5, 2024
1 parent a37a022 commit a04c75a
Show file tree
Hide file tree
Showing 3 changed files with 26 additions and 20 deletions.
36 changes: 19 additions & 17 deletions src/time_integrator.f90
Original file line number Diff line number Diff line change
Expand Up @@ -61,10 +61,10 @@ function init(backend, allocator, method, nvars)
type(time_intg_t) :: init
class(base_backend_t), pointer :: backend
class(allocator_t), pointer :: allocator
integer, intent(in), optional :: method
character(3), intent(in), optional :: method
integer, intent(in), optional :: nvars

integer :: i, j
integer :: i, j, stat

! initialize Runge-Kutta coefficients
! rk1
Expand Down Expand Up @@ -110,24 +110,28 @@ function init(backend, allocator, method, nvars)
init%allocator => allocator

if (present(method)) then
init%method = method
init%sname = method
else
init%method = 3
init%sname = 'AB3'
end if
if (init%method < 5) then
! method 1 to 4 -> AB1 to AB4
init%order = init%method
init%nstep = init%method

if (init%sname(1:2) == 'AB') then
read (init%sname(3:3), *, iostat=stat) init%order
if (stat /= 0) error stop 'Error reading AB integration order'
if (init%order >= 5) error stop 'Integration order >4 is not supported'
init%nstep = init%order
init%nstage = 1
init%nolds = init%nstep - 1
write (init%sname, "(A2,I1)") "AB", init%order
else
! method 5 to 8 -> RK1 to RK4
init%order = init%method - 4
else if (init%sname(1:2) == 'RK') then
read (init%sname(3:3), *, iostat=stat) init%order
if (stat /= 0) error stop 'Error reading RK integration order'
if (init%order >= 5) error stop 'Integration order >4 is not supported'
init%nstep = 1
init%nstage = init%method - 4
init%nstage = init%order
init%nolds = init%nstage
write (init%sname, "(A2,I1)") "RK", init%order
else
print *, 'Integration method '//init%sname//' is not defined'
error stop
end if

if (present(nvars)) then
Expand All @@ -151,8 +155,6 @@ function init(backend, allocator, method, nvars)
end do
end do

print *, init%sname, ' time integrator instantiated'

end function init

subroutine step(self, u, v, w, du, dv, dw, dt)
Expand All @@ -174,7 +176,7 @@ subroutine step(self, u, v, w, du, dv, dw, dt)
self%deriv(2)%ptr => dv
self%deriv(3)%ptr => dw

if (self%method < 5) then
if (self%sname(1:2) == 'AB' ) then
call self%adams_bashforth(dt)
else
call self%runge_kutta(dt)
Expand Down
5 changes: 3 additions & 2 deletions src/xcompact.f90
Original file line number Diff line number Diff line change
Expand Up @@ -103,8 +103,9 @@ program xcompact
if (nrank == 0) print *, 'OpenMP backend instantiated'
#endif

time_integrator = time_intg_t(allocator=allocator, backend=backend)
if (nrank == 0) print *, 'time integrator instantiated'
time_integrator = time_intg_t(allocator=allocator, backend=backend, &
method='AB3')
if (nrank == 0) print *, time_integrator%sname//' time intg instantiated'
solver = solver_t(backend, mesh, time_integrator, host_allocator)
if (nrank == 0) print *, 'solver instantiated'

Expand Down
5 changes: 4 additions & 1 deletion tests/test_time_integrator.f90
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ program test_omp_adamsbashforth
logical :: allpass = .true.
integer :: i, j, k, stage, istartup, nrank, nproc, ierr
integer :: nstep0 = 64, nstep, nrun = 4, nmethod = 8
character(len=3) :: method(8)
real(dp), allocatable, dimension(:) :: err
real(dp), allocatable, dimension(:) :: norm
real(dp) :: dt0 = 0.01_dp, dt, order
Expand Down Expand Up @@ -95,12 +96,14 @@ program test_omp_adamsbashforth

allocate (norm(nrun))

method = ['AB1', 'AB2', 'AB3', 'AB4', 'RK1', 'RK2', 'RK3', 'RK4']

! compute l2 norm for various step sizes
do k = 1, nmethod

! initialize time-integrator
time_integrator = time_intg_t(allocator=allocator, &
backend=backend, method=k)
backend=backend, method=method(k))

dt = dt0
nstep = nstep0
Expand Down

0 comments on commit a04c75a

Please sign in to comment.