timedisc_rkfehlberg.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: timedisc_rkfehlberg .f90 #
5!# #
6!# Copyright (C) 2011-2024 #
7!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
8!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
9!# #
10!# This program is free software; you can redistribute it and/or modify #
11!# it under the terms of the GNU General Public License as published by #
12!# the Free Software Foundation; either version 2 of the License, or (at #
13!# your option) any later version. #
14!# #
15!# This program is distributed in the hope that it will be useful, but #
16!# WITHOUT ANY WARRANTY; without even the implied warranty of #
17!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18!# NON INFRINGEMENT. See the GNU General Public License for more #
19!# details. #
20!# #
21!# You should have received a copy of the GNU General Public License #
22!# along with this program; if not, write to the Free Software #
23!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24!# #
25!#############################################################################
26
27!----------------------------------------------------------------------------!
41!----------------------------------------------------------------------------!
53 USE common_dict
54#ifdef PARALLEL
55#ifdef HAVE_MPI_MOD
56 USE mpi
57#endif
58#endif
59 IMPLICIT NONE
60#ifdef PARALLEL
61#ifdef HAVE_MPIF_H
62 include 'mpif.h'
63#endif
64#endif
65 !--------------------------------------------------------------------------!
66 PRIVATE
68 CLASS(marray_compound), POINTER :: p => null()
69 END TYPE coeff_type
71 TYPE(coeff_type), DIMENSION(:), ALLOCATABLE :: coeff
72 REAL, DIMENSION(:), POINTER :: b_low,b_high,c
73 REAL, DIMENSION(:,:), POINTER :: a
74 CONTAINS
75 PROCEDURE :: inittimedisc
77 PROCEDURE :: solveode
78 PROCEDURE :: finalize
80 generic :: computecvar => computecvar_rkfehlberg
81 PROCEDURE :: setbutchertableau
82 PROCEDURE :: showbutchertableau
83 END TYPE timedisc_rkfehlberg
84 CHARACTER(LEN=32), PARAMETER :: odesolver_name = "Runge-Kutta Fehlberg"
85
86 !--------------------------------------------------------------------------!
87 PUBLIC :: &
88 ! types
89 coeff_type, &
91 !--------------------------------------------------------------------------!
92
93CONTAINS
94
95 SUBROUTINE inittimedisc(this,Mesh,Physics,config,IO,ttype,tname)
98 IMPLICIT NONE
99 !------------------------------------------------------------------------!
100 CLASS(timedisc_rkfehlberg), INTENT(INOUT) :: this
101 CLASS(mesh_base), INTENT(INOUT) :: Mesh
102 CLASS(physics_base), INTENT(IN) :: Physics
103 TYPE(dict_typ), POINTER :: config,IO
104 INTEGER, INTENT(IN) :: ttype
105 CHARACTER(LEN=32), INTENT(IN) :: tname
106 !------------------------------------------------------------------------!
107 INTEGER :: err,k,ShowBT
108 !------------------------------------------------------------------------!
109 CALL this%timedisc_modeuler%InitTimedisc(mesh,physics,config,io,ttype,tname)
110
111 ALLOCATE(this%b_high(this%m),this%b_low(this%m),&
112 this%c(this%m),this%a(this%m,this%m), &
113 this%coeff(this%m), &
114 stat = err)
115 IF (err.NE.0) THEN
116 CALL this%Error("timedisc_rkfehlberg::InitTimedisc", "Unable to allocate memory.")
117 END IF
118
119 ! init Butcher tableau
120 CALL this%SetButcherTableau()
121 CALL getattr(config, "ShowButcherTableau", showbt, 0)
122 IF (showbt.EQ.1) CALL this%ShowButcherTableau()
123
124 ! init coefficient compounds
125 this%coeff(1)%p => this%rhs
126 DO k=2,this%m
127 CALL physics%new_statevector(this%coeff(k)%p,conservative)
128 this%coeff(k)%p%data1d(:) = 0.0
129 END DO
130
131 IF ((this%tol_rel.LT.0.0).OR.minval(this%tol_abs(:)).LT.0.0) &
132 CALL this%Error("timedisc_rkfehlberg::InitTimedisc", &
133 "error tolerance levels must be greater than 0")
134
135 IF (this%tol_rel.GT.1.0) THEN
136 CALL this%Warning("timedisc_rkfehlberg::InitTimedisc", &
137 "adaptive step size control disabled (tol_rel>1)",0)
138 ELSE IF(this%tol_rel.GE.0.01 .AND. this%order .GE. 5) THEN
139 CALL this%Warning("timedisc_rkfehlberg::InitTimedisc", &
140 "You chose a relatively high tol_rel (in comparison to order)",0)
141 END IF
142
143 END SUBROUTINE
144
145 SUBROUTINE inittimedisc_rkfehlberg(this,Mesh,Physics,config,IO)
146 IMPLICIT NONE
147 !------------------------------------------------------------------------!
148 CLASS(timedisc_rkfehlberg), INTENT(INOUT) :: this
149 CLASS(mesh_base), INTENT(INOUT) :: Mesh
150 CLASS(physics_base), INTENT(IN) :: Physics
151 TYPE(dict_typ), POINTER :: config,IO
152 !------------------------------------------------------------------------!
153 ! set default order
154 CALL getattr(config, "order", this%order, 5)
155
156 ! set number of coefficients
157 SELECT CASE(this%GetOrder())
158 CASE(3)
159 this%m = 3
160 CASE(5)
161 this%m = 6
162 CASE DEFAULT
163 CALL this%Error("InitTimedisc_rkfehlberg","time order must be 3 or 5")
164 END SELECT
165
166 CALL this%InitTimedisc(mesh,physics,config,io,rk_fehlberg,odesolver_name)
167 END SUBROUTINE inittimedisc_rkfehlberg
168
169
170 SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
171 IMPLICIT NONE
172 CLASS(timedisc_rkfehlberg),INTENT(INOUT):: this
173 CLASS(mesh_base), INTENT(IN) :: Mesh
174 CLASS(physics_base), INTENT(INOUT) :: Physics
175 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
176 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
177 REAL, INTENT(IN) :: time
178 REAL, INTENT(INOUT) :: dt, err
179 !------------------------------------------------------------------------!
180 INTEGER :: i,j,k,l,m
181 REAL :: t
182 !------------------------------------------------------------------------!
183 t = time
184 ! ATTENTION: rhs and therefore coeff(1)%p should already be up to date.
185 ! In the very first step this is done in FirstStepFosite, any later step
186 ! must update coeff(1)%p at the end of the whole time step update. This is
187 ! done in timedisc_base::AcceptTimestep or timedisc_base::RejectTimestep.
188 ! Don't forget: this%rhs => this%coeff(1)%p.
189!NEC$ SHORTLOOP
190 DO m=2,this%m
191 ! time step update of cell mean values
192 CALL this%ComputeCVar_rkfehlberg(mesh,physics,fluxes,dt,m,this%coeff,this%cvar,this%ctmp)
193 ! compute right-hand-side
194 ! coeff_m is k_m/dt from Butcher tableau
195 CALL this%ComputeRHS(mesh,physics,sources,fluxes,t+this%c(m)*dt,dt,&
196 this%ptmp,this%ctmp,check_nothing,this%coeff(m)%p)
197 END DO
198
199 !reset ctmp
200 this%ctmp = this%cvar
201!NEC$ SHORTLOOP
202 DO m=1,this%m
203 ! compute two solutions with different numerical orders
204 ! y_n+1 = y_n + SUM(b_i*k_i) = y_n + SUM(b_i*dt*coeff_i)
205 DO concurrent(i=1:SIZE(this%ctmp%data1d))
206 this%ctmp%data1d(i) = this%ctmp%data1d(i) &
207 - this%b_low(m)*dt*this%coeff(m)%p%data1d(i)
208 this%cvar%data1d(i) = this%cvar%data1d(i) &
209 - this%b_high(m)*dt*this%coeff(m)%p%data1d(i)
210 END DO
211 END DO
212
213 ! at the boundary the this%rhs contains the boundary fluxes
214 ! (see subroutine ComputeRHS_modeuler )
215!NEC$ SHORTLOOP
216 DO m=1,this%m
217!NEC$ SHORTLOOP
218 DO l=1,physics%VNUM
219 IF(mesh%INUM.GT.1) THEN
220 ! western and eastern
221!NEC$ IVDEP
222 DO k=mesh%KMIN,mesh%KMAX
223!NEC$ IVDEP
224 DO j=mesh%JMIN,mesh%JMAX
225 fluxes%bxflux(j,k,1,l) = fluxes%bxflux(j,k,1,l) &
226 - dt*this%b_high(m)*this%coeff(m)%p%data4d(mesh%IMIN-mesh%IP1,j,k,l)
227 fluxes%bxflux(j,k,2,l) = fluxes%bxflux(j,k,2,l) &
228 - dt*this%b_high(m)*this%coeff(m)%p%data4d(mesh%IMAX+mesh%IP1,j,k,l)
229 END DO
230 END DO
231 END IF
232 ! southern and northern
233 IF(mesh%JNUM.GT.1) THEN
234!NEC$ IVDEP
235 DO k=mesh%KMIN,mesh%KMAX
236!NEC$ IVDEP
237 DO i=mesh%IMIN,mesh%IMAX
238 fluxes%byflux(k,i,1,l) = fluxes%byflux(k,i,1,l) &
239 - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,mesh%JMIN-mesh%JP1,k,l)
240 fluxes%byflux(k,i,2,l) = fluxes%byflux(k,i,2,l) &
241 - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,mesh%JMAX+mesh%JP1,k,l)
242 END DO
243 END DO
244 END IF
245 ! bottom and top
246 IF(mesh%KNUM.GT.1) THEN
247!NEC$ IVDEP
248 DO j=mesh%JMIN,mesh%JMAX
249!NEC$ IVDEP
250 DO i=mesh%IMIN,mesh%IMAX
251 fluxes%bzflux(i,j,1,l) = fluxes%bzflux(i,j,1,l) &
252 - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,j,mesh%KMIN-mesh%KP1,l)
253 fluxes%bzflux(i,j,2,l) = fluxes%bzflux(i,j,2,l) &
254 - dt*this%b_high(m)*this%coeff(m)%p%data4d(i,j,mesh%KMAX+mesh%KP1,l)
255 END DO
256 END DO
257 END IF
258 END DO
259 END DO
260
261 ! compute an error estimate based on the two independent numerical
262 ! solutions err and dt
263 IF (this%tol_rel.LE.1.0) THEN
264 err = this%ComputeError(mesh,physics,this%cvar,this%ctmp)
265 dt = this%AdjustTimestep(err,dt)
266 END IF
267
268 END SUBROUTINE solveode
269
277 SUBROUTINE computecvar_rkfehlberg(this,Mesh,Physics,Fluxes,dt,m,coeff,cvar,cnew)
278 IMPLICIT NONE
279 !------------------------------------------------------------------------!
280 CLASS(timedisc_rkfehlberg),INTENT(INOUT) :: this
281 CLASS(mesh_base) ,INTENT(IN) :: Mesh
282 CLASS(physics_base) ,INTENT(INOUT) :: Physics
283 CLASS(fluxes_base) ,INTENT(INOUT) :: Fluxes
284 REAL ,INTENT(IN) :: dt
285 INTEGER ,INTENT(IN) :: m
286 TYPE(coeff_type), DIMENSION(m), INTENT(IN) :: coeff
287 CLASS(marray_compound) ,INTENT(INOUT) :: cvar,cnew
288 !------------------------------------------------------------------------!
289 INTEGER :: mm
290 !------------------------------------------------------------------------!
291 ! compute y + SUM(a_mi*k_i) = y + SUM(a_mi*dt*rhs_i) : i in [1,m-1]
292 ! see Runge-Kutta method (Butcher tableau)
293 IF (m.EQ.1) THEN
294 cnew = cvar
295 ELSE
296 cnew%data1d(:) = cvar%data1d(:) - this%a(m,1)*dt*coeff(1)%p%data1d(:)
297!NEC$ SHORTLOOP
298 DO mm=2,m-1
299 cnew%data1d(:) = cnew%data1d(:) - this%a(m,mm)*dt*coeff(mm)%p%data1d(:)
300 END DO
301 END IF
302 END SUBROUTINE computecvar_rkfehlberg
303
305 SUBROUTINE setbutchertableau(this)
306 IMPLICIT NONE
307 !------------------------------------------------------------------------!
308 CLASS(timedisc_rkfehlberg) :: this
309 !------------------------------------------------------------------------!
310 SELECT CASE(this%GetOrder())
311 CASE(3) ! RK-Fehlberg rkf23
312 this%b_high = (/ 1.0/6.0, 2.0/3.0, 1.0/6.0 /)
313 this%b_low = (/ 0.0, 1.0, 0.0 /)
314 this%c = (/ 0.0, 0.5, 1.0 /)
315 this%a = transpose(reshape((/ &
316 0.0, 0.0, 0.0, &
317 0.5, 0.0, 0.0, &
318 -1.0, 2.0, 0.0 /),(/this%m,this%m/)))
319 CASE(5) ! RK-Fehlberg rkf45
320 this%b_high = (/ 16.0/135.0, 0.0, 6656.0/12825.0, &
321 28561.0/56430.0, -9.0/50.0, 2.0/55.0 /)
322 this%b_low = (/ 25.0/216.0, 0.0, 1408.0/2565.0, 2197.0/4104.0, -0.2, 0.0 /)
323 this%c = (/ 0.0, 0.25, 3.0/8.0, 12.0/13.0, 1.0, 0.5 /)
324 this%a = transpose(reshape((/ &
325 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
326 0.25, 0.0, 0.0, 0.0, 0.0, 0.0, &
327 3.0/32.0, 9.0/32.0, 0.0, 0.0, 0.0, 0.0, &
328 1932.0/2197.0, -7200.0/2197.0, 7296.0/2197.0, 0.0, 0.0, 0.0, &
329 439.0/216.0, -8.0, 3680.0/513.0, -845.0/4104.0, 0.0, 0.0, &
330 -8.0/27.0, 2.0, -3544.0/2565.0, 1859.0/4104.0, -11.0/40.0, 0.0/),&
331 (/this%m,this%m/)))
332 CASE DEFAULT
333 CALL this%Error("timedisc_rkfehlberg::SetButcherTableau","only order 3 or 5 supported")
334 END SELECT
335 END SUBROUTINE setbutchertableau
336
337
338 SUBROUTINE showbutchertableau(this)
339 IMPLICIT NONE
340 !------------------------------------------------------------------------!
341 CLASS(timedisc_rkfehlberg) :: this
342 !------------------------------------------------------------------------!
343 INTEGER :: i
344 CHARACTER(LEN=64) :: sformat
345 CHARACTER(LEN=128) :: buffer
346 !------------------------------------------------------------------------!
347 INTENT(INOUT) :: this
348 !------------------------------------------------------------------------!
349 CALL this%Info(repeat("+",(this%m+1)*12+4))
350 CALL this%Info("Butcher tableau")
351 DO i=1,this%m
352 IF (i == 1) THEN
353 WRITE (sformat,'(A)') '(ES12.3,A)'
354 WRITE (buffer,fmt=sformat) this%c(i), " | "
355 ELSE
356 WRITE (sformat,'(A,I1,A)') '(ES12.3,A,',i-1,'(ES12.3))'
357 WRITE (buffer,fmt=sformat) this%c(i), " | " , this%a(i,1:i-1)
358 END IF
359 CALL this%Info(buffer)
360 END DO
361 CALL this%Info(repeat("-",(this%m+1)*12+4))
362 WRITE (sformat,'(A,I1,A)') '(A,',this%m,'(ES12.3))'
363 WRITE (buffer,fmt=sformat) " high order | ", this%b_high(:)
364 CALL this%Info(buffer)
365 WRITE (buffer,fmt=sformat) " low order | ", this%b_low(:)
366 CALL this%Info(buffer)
367 CALL this%Info(repeat("+",(this%m+1)*12+4))
368 WRITE (buffer,fmt=sformat) "c_i-SUM(a_ij)| ", this%c(:)-sum(this%a(:,:),dim=2)
369 CALL this%Info(buffer)
370 END SUBROUTINE showbutchertableau
371
372
373 SUBROUTINE finalize(this)
374 IMPLICIT NONE
375 !------------------------------------------------------------------------!
376 CLASS(timedisc_rkfehlberg) :: this
377 !------------------------------------------------------------------------!
378 INTEGER :: k
379 !------------------------------------------------------------------------!
380 DEALLOCATE(this%b_high,this%b_low,this%c,this%a)
381 DO k=2,this%m
382 DEALLOCATE(this%coeff(k)%p)
383 END DO
384 DEALLOCATE(this%coeff)
385 CALL this%timedisc_modeuler%Finalize()
386 END SUBROUTINE finalize
387
Dictionary for generic data types.
Definition: common_dict.f90:61
base module for numerical flux functions
Definition: fluxes_base.f90:39
base class for mesh arrays
Definition: marray_base.f90:36
derived class for compound of mesh arrays
basic mesh module
Definition: mesh_base.f90:72
Basic physics module.
@, public conservative
physics module for 1D,2D and 3D non-isothermal Euler equations
physics module for 1D,2D and 3D isothermal Euler equations
generic source terms module providing functionaly common to all source terms
module to manage list of source terms
integer, parameter, public check_nothing
subroutine inittimedisc(this, Mesh, Physics, config, IO, ttype, tname)
integer, parameter, public rk_fehlberg
subroutines for modified Euler i.e. Runge-Kutta methods
subroutine solveode(this, Mesh, Physics, Sources, Fluxes, time, dt, err)
character(len=32), parameter odesolver_name
subroutines for Runge-Kutta Fehlberg method
subroutine computecvar_rkfehlberg(this, Mesh, Physics, Fluxes, dt, m, coeff, cvar, cnew)
perfroms the time step update using the RHS
subroutine inittimedisc_rkfehlberg(this, Mesh, Physics, config, IO)
subroutine setbutchertableau(this)
set coefficients for RK-Fehlberg schemes
subroutine showbutchertableau(this)
mesh data structure
Definition: mesh_base.f90:122
container class to manage the list of source terms