74 REAL,
PARAMETER ::
eta(3,3) = &
75 reshape((/ 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.75, 1.0/3.0 /),(/3,3/))
76 REAL,
PARAMETER ::
zeta(3,3) = &
77 reshape((/ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5 /),(/3,3/))
96 TYPE(
dict_typ),
POINTER :: config, IO
99 CALL getattr(config,
"order", this%order, 3)
103 SELECT CASE(this%GetOrder())
109 this%tol_abs(:) = 1.0
110 CALL this%Warning(
"InitTimedisc_modeuler", &
111 "adaptive step size control not supported in 1st order scheme",0)
113 IF (this%tol_rel.GT.1.0) &
114 CALL this%Warning(
"InitTimedisc_modeuler", &
115 "adaptive step size control disabled (tol_rel>1)",0)
117 CALL this%Error(
"InitTimedisc_modeuler",
"time order must be one of 1,2,3")
119 IF ((this%tol_rel.LT.0.0).OR.minval(this%tol_abs(:)).LT.0.0) &
120 CALL this%Error(
"InitTimedisc_modeuler", &
121 "error tolerance levels must be greater than 0")
125 SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
130 CLASS(
sources_list),
ALLOCATABLE,
INTENT(INOUT) :: Sources
132 REAL,
INTENT(IN) :: time
133 REAL,
INTENT(INOUT) :: dt, err
141 TYPE(var_typ) :: p(4),c(4)
144 order = this%GetOrder()
146 IF (this%tol_rel.GE.1.0)
THEN
151 t = time+
zeta(n,order)*dt
153 CALL this%ComputeCVar(mesh,physics,fluxes,
eta(n,order), &
154 t,dt,this%cold,this%pvar,this%cvar,this%rhs,this%cvar)
157 CALL this%ComputeRHS(mesh,physics,sources,fluxes,t,dt,&
164 p(1)%var => this%pvar
165 p(2)%var => this%ptmp
166 p(3)%var => this%pvar
167 p(4)%var => this%pvar
168 c(1)%var => this%cvar
169 c(2)%var => this%ctmp
170 c(3)%var => this%cvar
171 c(4)%var => this%cvar
175 t = time+
zeta(n,order)*dt
177 CALL this%ComputeCVar(mesh,physics,fluxes,
eta(n,order), &
178 t,dt,this%cold,p(n)%var,c(n)%var,this%rhs,c(n+1)%var)
181 IF (n.EQ.2.AND.order.EQ.3) &
182 this%ctmp%data4d(:,:,:,:) =
updatetimestep_modeuler(
eta(2,2),dt,this%cold%data4d(:,:,:,:), &
183 this%ctmp%data4d(:,:,:,:),this%rhs%data4d(:,:,:,:))
186 CALL this%ComputeRHS(mesh,physics,sources,fluxes,t,dt,p(n+1)%var,c(n+1)%var,&
190 err = this%ComputeError(mesh,physics,this%cvar,this%ctmp)
191 dt = this%AdjustTimestep(err,dt)
201 cold,pvar,cvar,rhs,cnew)
213 INTENT(IN) :: eta,time,dt
221 cold%data2d(:,l),cvar%data2d(:,l),rhs%data2d(:,l))
223 IF(mesh%INUM.GT.1)
THEN
226 DO k=mesh%KMIN,mesh%KMAX
228 DO j=mesh%JMIN,mesh%JMAX
231 fluxes%bxflux(j,k,1,l),rhs%data4d(mesh%IMIN-mesh%IP1,j,k,l))
233 fluxes%bxflux(j,k,2,l),rhs%data4d(mesh%IMAX+mesh%IP1,j,k,l))
237 IF(mesh%JNUM.GT.1)
THEN
240 DO i=mesh%IMIN,mesh%IMAX
242 DO k=mesh%KMIN,mesh%KMAX
245 fluxes%byflux(k,i,1,l),rhs%data4d(i,mesh%JMIN-mesh%Jp1,k,l))
247 fluxes%byflux(k,i,2,l),rhs%data4d(i,mesh%JMAX+mesh%Jp1,k,l))
252 IF(mesh%KNUM.GT.1)
THEN
255 DO j=mesh%JMIN,mesh%JMAX
258 DO i=mesh%IMIN,mesh%IMAX
260 fluxes%bzflux(i,j,1,l),rhs%data4d(i,j,mesh%KMIN-mesh%Kp1,l))
262 fluxes%bzflux(i,j,2,l),rhs%data4d(i,j,mesh%KMAX+mesh%Kp1,l))
275 CALL this%Finalize_base()
282 REAL,
INTENT(IN) :: eta_n,dt,y0,yn,rhs
290 y = yn-dt*rhs+eta_n*(y0-yn+dt*rhs)
Dictionary for generic data types.
base module for numerical flux functions
base class for mesh arrays
derived class for compound of mesh arrays
generic source terms module providing functionaly common to all source terms
module to manage list of source terms
integer, parameter, public dtcause_smallerr
integer, parameter, public check_tmin
integer, parameter, public check_csound
integer, parameter, public check_invalid
integer, parameter, public check_nothing
integer, parameter, public check_rhomin
integer, parameter, public dtcause_erradj
integer, parameter, public check_all
integer, parameter, public modified_euler
integer, parameter, public dtcause_cfl
integer, parameter, public check_pmin
subroutines for modified Euler i.e. Runge-Kutta methods
real, dimension(3, 3), parameter zeta
subroutine inittimedisc_modeuler(this, Mesh, Physics, config, IO)
elemental real function updatetimestep_modeuler(eta_n, dt, y0, yn, rhs)
subroutine computecvar_modeuler(this, Mesh, Physics, Fluxes, eta, time, dt, cold, pvar, cvar, rhs, cnew)
performs the time step update using the RHS
real, dimension(3, 3), parameter eta
subroutine solveode(this, Mesh, Physics, Sources, Fluxes, time, dt, err)
subroutine finalize(this)
character(len=32), parameter odesolver_name
container class to manage the list of source terms