timedisc_modeuler.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: timedisc_modeuler.f90 #
5!# #
6!# Copyright (C) 2007-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
27
28!----------------------------------------------------------------------------!
39!----------------------------------------------------------------------------!
50 USE common_dict
51#ifdef PARALLEL
52#ifdef HAVE_MPI_MOD
53 USE mpi
54#endif
55#endif
56 IMPLICIT NONE
57#ifdef PARALLEL
58#ifdef HAVE_MPIF_H
59 include 'mpif.h'
60#endif
61#endif
62 !--------------------------------------------------------------------------!
63 PRIVATE
65 CONTAINS
68 generic :: computecvar => computecvar_modeuler
69 PROCEDURE :: solveode
70 PROCEDURE :: finalize
71 END TYPE timedisc_modeuler
72 !--------------------------------------------------------------------------!
73 CHARACTER(LEN=32), PARAMETER :: odesolver_name = "modified Euler"
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/))
78 !--------------------------------------------------------------------------!
79 PUBLIC :: &
80 ! types
82 ! constants
86 !--------------------------------------------------------------------------!
87
88CONTAINS
89
90 SUBROUTINE inittimedisc_modeuler(this,Mesh,Physics,config,IO)
91 IMPLICIT NONE
92 !------------------------------------------------------------------------!
93 CLASS(timedisc_modeuler), INTENT(INOUT) :: this
94 CLASS(mesh_base), INTENT(INOUT) :: Mesh
95 CLASS(physics_base), INTENT(IN) :: Physics
96 TYPE(dict_typ), POINTER :: config, IO
97 !------------------------------------------------------------------------!
98 ! set default order
99 CALL getattr(config, "order", this%order, 3)
100
101 CALL this%InitTimedisc(mesh,physics,config,io,modified_euler,odesolver_name)
102
103 SELECT CASE(this%GetOrder())
104 CASE(1)
105 ! set relative error tolarance to value > 1.0
106 ! to disable adaptive step size control
107 ! (not available for 1st order scheme)
108 this%tol_rel = 10.0
109 this%tol_abs(:) = 1.0
110 CALL this%Warning("InitTimedisc_modeuler", &
111 "adaptive step size control not supported in 1st order scheme",0)
112 CASE(2,3)
113 IF (this%tol_rel.GT.1.0) &
114 CALL this%Warning("InitTimedisc_modeuler", &
115 "adaptive step size control disabled (tol_rel>1)",0)
116 CASE DEFAULT
117 CALL this%Error("InitTimedisc_modeuler","time order must be one of 1,2,3")
118 END SELECT
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")
122
123 END SUBROUTINE inittimedisc_modeuler
124
125 SUBROUTINE solveode(this,Mesh,Physics,Sources,Fluxes,time,dt,err)
126 IMPLICIT NONE
127 CLASS(timedisc_modeuler), INTENT(INOUT) :: this
128 CLASS(mesh_base), INTENT(IN) :: Mesh
129 CLASS(physics_base), INTENT(INOUT) :: Physics
130 CLASS(sources_list), ALLOCATABLE, INTENT(INOUT) :: Sources
131 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
132 REAL, INTENT(IN) :: time
133 REAL, INTENT(INOUT) :: dt, err
134 !------------------------------------------------------------------------!
135 INTEGER :: n
136 INTEGER :: order
137 REAL :: t
138 TYPE var_typ
139 CLASS(marray_compound), POINTER :: var
140 END TYPE var_typ
141 TYPE(var_typ) :: p(4),c(4)
142 !------------------------------------------------------------------------!
143 t = time
144 order = this%GetOrder()
145 ! check if adaptive step size control is enabled
146 IF (this%tol_rel.GE.1.0) THEN
147 ! no adaptive step size control
148!NEC$ UNROLL(3)
149 DO n=1,order
150 ! update time variable
151 t = time+zeta(n,order)*dt
152 ! time step update of cvar and bfluxes
153 CALL this%ComputeCVar(mesh,physics,fluxes,eta(n,order), &
154 t,dt,this%cold,this%pvar,this%cvar,this%rhs,this%cvar)
155 ! compute right hand side for next time step update except for last step
156 IF (n.LT.order) &
157 CALL this%ComputeRHS(mesh,physics,sources,fluxes,t,dt,&
158 this%pvar,this%cvar,check_nothing,this%rhs)
159 END DO
160 err = 0.0
161 dt = huge(dt)
162 ELSE
163 ! with adaptive step size control (2nd / 3rd order scheme)
164 p(1)%var => this%pvar
165 p(2)%var => this%ptmp ! store intermediate result for error control
166 p(3)%var => this%pvar
167 p(4)%var => this%pvar
168 c(1)%var => this%cvar
169 c(2)%var => this%ctmp ! store intermediate result for error control
170 c(3)%var => this%cvar
171 c(4)%var => this%cvar
172!NEC$ UNROLL(3)
173 DO n=1,order
174 ! update time variable
175 t = time+zeta(n,order)*dt
176 ! time step update of cvar and bfluxes
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)
179 ! for 3rd order scheme compute the 2nd order result with the same RHS
180 ! and store it in this%ctmp, bfluxes are not required
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(:,:,:,:))
184 ! compute right hand side for next time step update except for last step
185 IF (n.LT.order) &
186 CALL this%ComputeRHS(mesh,physics,sources,fluxes,t,dt,p(n+1)%var,c(n+1)%var,&
187 check_nothing,this%rhs)
188 END DO
189
190 err = this%ComputeError(mesh,physics,this%cvar,this%ctmp)
191 dt = this%AdjustTimestep(err,dt)
192
193 END IF
194
195 END SUBROUTINE solveode
196
197
200 SUBROUTINE computecvar_modeuler(this,Mesh,Physics,Fluxes,eta,time,dt, &
201 cold,pvar,cvar,rhs,cnew)
202 IMPLICIT NONE
203 !------------------------------------------------------------------------!
204 CLASS(timedisc_modeuler), INTENT(INOUT) :: this
205 CLASS(mesh_base), INTENT(IN) :: Mesh
206 CLASS(physics_base), INTENT(INOUT) :: Physics
207 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
208 CLASS(marray_compound), INTENT(INOUT) :: cold,pvar,cvar,cnew,rhs
209 REAL :: eta,time,dt
210 !------------------------------------------------------------------------!
211 INTEGER :: i,j,k,l
212 !------------------------------------------------------------------------!
213 INTENT(IN) :: eta,time,dt
214 !------------------------------------------------------------------------!
215!NEC$ NOVECTOR
216 DO l=1,physics%VNUM
217 ! use collapsed arrays for better vector performance
218 ! ATTENTION: yields garbage in ghost zones
219 ! -> overwritten in CenterBoundary called in timedisc_base::ComputeRHS
220 cnew%data2d(:,l) = updatetimestep_modeuler(eta,dt, &
221 cold%data2d(:,l),cvar%data2d(:,l),rhs%data2d(:,l))
222
223 IF(mesh%INUM.GT.1) THEN
224 ! western and eastern boundary fluxes
225!NEC$ IVDEP
226 DO k=mesh%KMIN,mesh%KMAX
227!NEC$ IVDEP
228 DO j=mesh%JMIN,mesh%JMAX
229 ! time step update of boundary fluxes
230 fluxes%bxflux(j,k,1,l) = updatetimestep_modeuler(eta,dt,fluxes%bxfold(j,k,1,l), &
231 fluxes%bxflux(j,k,1,l),rhs%data4d(mesh%IMIN-mesh%IP1,j,k,l))
232 fluxes%bxflux(j,k,2,l) = updatetimestep_modeuler(eta,dt,fluxes%bxfold(j,k,2,l), &
233 fluxes%bxflux(j,k,2,l),rhs%data4d(mesh%IMAX+mesh%IP1,j,k,l))
234 END DO
235 END DO
236 END IF
237 IF(mesh%JNUM.GT.1) THEN
238 ! southern and northern boundary fluxes
239!NEC$ IVDEP
240 DO i=mesh%IMIN,mesh%IMAX
241!NEC$ IVDEP
242 DO k=mesh%KMIN,mesh%KMAX
243 ! time step update of boundary fluxes
244 fluxes%byflux(k,i,1,l) = updatetimestep_modeuler(eta,dt,fluxes%byfold(k,i,1,l), &
245 fluxes%byflux(k,i,1,l),rhs%data4d(i,mesh%JMIN-mesh%Jp1,k,l))
246 fluxes%byflux(k,i,2,l) = updatetimestep_modeuler(eta,dt,fluxes%byfold(k,i,2,l), &
247 fluxes%byflux(k,i,2,l),rhs%data4d(i,mesh%JMAX+mesh%Jp1,k,l))
248 END DO
249 END DO
250 END IF
251
252 IF(mesh%KNUM.GT.1) THEN
253 ! bottom and top boundary fluxes
254!NEC$ IVDEP
255 DO j=mesh%JMIN,mesh%JMAX
256 ! time step update of boundary fluxes
257!NEC$ IVDEP
258 DO i=mesh%IMIN,mesh%IMAX
259 fluxes%bzflux(i,j,1,l) = updatetimestep_modeuler(eta,dt,fluxes%bzfold(i,j,1,l), &
260 fluxes%bzflux(i,j,1,l),rhs%data4d(i,j,mesh%KMIN-mesh%Kp1,l))
261 fluxes%bzflux(i,j,2,l) = updatetimestep_modeuler(eta,dt,fluxes%bzfold(i,j,2,l), &
262 fluxes%bzflux(i,j,2,l),rhs%data4d(i,j,mesh%KMAX+mesh%Kp1,l))
263 END DO
264 END DO
265 END IF
266 END DO
267 END SUBROUTINE computecvar_modeuler
268
269
270 SUBROUTINE finalize(this)
271 IMPLICIT NONE
272 !------------------------------------------------------------------------!
273 CLASS(timedisc_modeuler) :: this
274 !------------------------------------------------------------------------!
275 CALL this%Finalize_base()
276 END SUBROUTINE finalize
277
278
279 ELEMENTAL FUNCTION updatetimestep_modeuler(eta_n,dt,y0,yn,rhs) RESULT(y)
280 IMPLICIT NONE
281 !------------------------------------------------------------------------!
282 REAL, INTENT(IN) :: eta_n,dt,y0,yn,rhs
283 REAL :: y
284 !------------------------------------------------------------------------!
285 ! ATTENTION:
286 ! The time step update is computed according to:
287 ! y = eta_n * y0 + (1.0 - eta_n) * (yn - dt * rhs)
288 ! but to minimize the truncation error it is essential to sort the terms
289 ! in this way:
290 y = yn-dt*rhs+eta_n*(y0-yn+dt*rhs)
291 END FUNCTION updatetimestep_modeuler
292
293END MODULE timedisc_modeuler_mod
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.
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)
character(len=32), parameter odesolver_name
mesh data structure
Definition: mesh_base.f90:122
container class to manage the list of source terms