72 REAL,
DIMENSION(:),
POINTER :: b_low,b_high,c
73 REAL,
DIMENSION(:,:),
POINTER :: a
103 TYPE(
dict_typ),
POINTER :: config,IO
104 INTEGER,
INTENT(IN) :: ttype
105 CHARACTER(LEN=32),
INTENT(IN) :: tname
107 INTEGER :: err,k,ShowBT
109 CALL this%timedisc_modeuler%InitTimedisc(mesh,physics,config,io,ttype,tname)
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), &
116 CALL this%Error(
"timedisc_rkfehlberg::InitTimedisc",
"Unable to allocate memory.")
120 CALL this%SetButcherTableau()
121 CALL getattr(config,
"ShowButcherTableau", showbt, 0)
122 IF (showbt.EQ.1)
CALL this%ShowButcherTableau()
125 this%coeff(1)%p => this%rhs
127 CALL physics%new_statevector(this%coeff(k)%p,
conservative)
128 this%coeff(k)%p%data1d(:) = 0.0
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")
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)
151 TYPE(
dict_typ),
POINTER :: config,IO
154 CALL getattr(config,
"order", this%order, 5)
157 SELECT CASE(this%GetOrder())
163 CALL this%Error(
"InitTimedisc_rkfehlberg",
"time order must be 3 or 5")
170 SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
175 CLASS(
sources_list),
ALLOCATABLE,
INTENT(INOUT) :: Sources
177 REAL,
INTENT(IN) :: time
178 REAL,
INTENT(INOUT) :: dt, err
192 CALL this%ComputeCVar_rkfehlberg(mesh,physics,fluxes,dt,m,this%coeff,this%cvar,this%ctmp)
195 CALL this%ComputeRHS(mesh,physics,sources,fluxes,t+this%c(m)*dt,dt,&
200 this%ctmp = this%cvar
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)
219 IF(mesh%INUM.GT.1)
THEN
222 DO k=mesh%KMIN,mesh%KMAX
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)
233 IF(mesh%JNUM.GT.1)
THEN
235 DO k=mesh%KMIN,mesh%KMAX
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)
246 IF(mesh%KNUM.GT.1)
THEN
248 DO j=mesh%JMIN,mesh%JMAX
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)
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)
284 REAL ,
INTENT(IN) :: dt
285 INTEGER ,
INTENT(IN) :: m
286 TYPE(
coeff_type),
DIMENSION(m),
INTENT(IN) :: coeff
296 cnew%data1d(:) = cvar%data1d(:) - this%a(m,1)*dt*coeff(1)%p%data1d(:)
299 cnew%data1d(:) = cnew%data1d(:) - this%a(m,mm)*dt*coeff(mm)%p%data1d(:)
310 SELECT CASE(this%GetOrder())
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((/ &
318 -1.0, 2.0, 0.0 /),(/this%m,this%m/)))
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/),&
333 CALL this%Error(
"timedisc_rkfehlberg::SetButcherTableau",
"only order 3 or 5 supported")
344 CHARACTER(LEN=64) :: sformat
345 CHARACTER(LEN=128) :: buffer
347 INTENT(INOUT) :: this
349 CALL this%Info(repeat(
"+",(this%m+1)*12+4))
350 CALL this%Info(
"Butcher tableau")
353 WRITE (sformat,
'(A)')
'(ES12.3,A)'
354 WRITE (buffer,fmt=sformat) this%c(i),
" | "
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)
359 CALL this%Info(buffer)
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)
380 DEALLOCATE(this%b_high,this%b_low,this%c,this%a)
382 DEALLOCATE(this%coeff(k)%p)
384 DEALLOCATE(this%coeff)
385 CALL this%timedisc_modeuler%Finalize()
Dictionary for generic data types.
base module for numerical flux functions
base class for mesh arrays
derived class for compound of mesh arrays
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)
subroutine finalize(this)
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)
container class to manage the list of source terms