ode.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: ode.f90 #
5!# #
6!# Copyright (C) 2006-2021 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Björn Sperling <sperling@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!----------------------------------------------------------------------------!
31
32!! References:
33!! \cite sedov1945 Sedov, L. I.: Unsteady motions of compressible fluids,
34!! J. Appl. Math. Mech. 9 (1945)
35!! \cite sedov1959 Sedov, L. I.: Similarity and Dimensional Methods in Mechanics
36!! Academic Press Ltd., New York (1959)
37!----------------------------------------------------------------------------!
38PROGRAM ode_test
39 USE fosite_mod
40#include "tap.h"
41 IMPLICIT NONE
42 !--------------------------------------------------------------------------!
43 ! simulation parameters
44 REAL, PARAMETER :: tsim = 0.05 ! simulation stop time
45 REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats
46 ! initial condition (dimensionless units)
47 REAL, PARAMETER :: rho0 = 1.0 ! ambient density
48 REAL, PARAMETER :: p0 = 1.0e-05 ! ambient pressure
49 REAL, PARAMETER :: e1 = 1.0 ! initial energy input
50 ! Spatial with of the initial pulse should be at least 5 cells;
51 ! if you wish to compare the results on different grids
52 ! R0 should be of the same order
53 REAL, PARAMETER :: r0 = 3.0e-2
54 ! mesh settings
55 INTEGER, PARAMETER :: mgeo = cartesian ! geometry
56 INTEGER, PARAMETER :: res = 64 ! x- & y-resolution
57 ! output parameters
58 INTEGER, PARAMETER :: onum = 5 ! number of output data sets
59 CHARACTER(LEN=256), PARAMETER & ! output data dir
60 :: odir = './'
61 CHARACTER(LEN=256), PARAMETER & ! output data file name
62 :: ofname = 'ode'
63 ! ode solvers
64 INTEGER, PARAMETER :: num_tests = 8
65 !--------------------------------------------------------------------------!
66 CLASS(fosite),ALLOCATABLE :: sim
67 INTEGER :: k
68 CHARACTER(LEN=80) :: testinfo(num_tests)
69 !--------------------------------------------------------------------------!
70
71 tap_plan(num_tests)
72
73 DO k=1,num_tests
74 ALLOCATE(sim)
75 CALL sim%InitFosite()
76
77 CALL makeconfig(sim, sim%config)
78
79 ! specific test settings
80 SELECT CASE(k)
81 CASE(1) ! Modified Euler, Order 2
82 CALL setattr(sim%config, "/timedisc/method", modified_euler)
83 CALL setattr(sim%config, "/timedisc/order", 2)
84 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-mdeu2") )
85 CASE(2) ! Modified Euler, Order 3
86 CALL setattr(sim%config, "/timedisc/method", modified_euler)
87 CALL setattr(sim%config, "/timedisc/order", 3)
88 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-mdeu3") )
89 CASE(3) ! RK-Fehlberg, Order 3
90 CALL setattr(sim%config, "/timedisc/method", rk_fehlberg)
91 CALL setattr(sim%config, "/timedisc/order", 3)
92 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-rkfeh3") )
93 CASE(4) ! RK-Fehlberg, Order 5
94 CALL setattr(sim%config, "/timedisc/method", rk_fehlberg)
95 CALL setattr(sim%config, "/timedisc/order", 5)
96 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-rkfeh5") )
97 CASE(5) ! Cash-Karp, Order 5
98 CALL setattr(sim%config, "/timedisc/method", cash_karp)
99 CALL setattr(sim%config, "/timedisc/order", 5)
100 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-cshka5") )
101 CASE(6) ! Dormand-Prince, Order 5
102 CALL setattr(sim%config, "/timedisc/method", dormand_prince)
103 CALL setattr(sim%config, "/timedisc/order", 5)
104 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-drmpr5") )
105 CASE(7) ! SSPRK, Order 5
106 CALL setattr(sim%config, "/timedisc/method", ssprk)
107 CALL setattr(sim%config, "/timedisc/order", 3)
108 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-ssprk3") )
109 CASE(8) ! SSPRK, Order 5
110 CALL setattr(sim%config, "/timedisc/method", ssprk)
111 CALL setattr(sim%config, "/timedisc/order", 5)
112 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-ssprk5") )
113 CASE DEFAULT
114 CALL sim%Error("ode","Unknown test number!")
115 END SELECT
116
117 ! CALL PrintDict(Sim%config)
118 CALL sim%Setup()
119 ! set initial condition
120 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
121
122 CALL sim%Run()
123
124 IF (.NOT. sim%aborted) THEN
125 tap_check(.true.,"Simulation finished")
126 CALL sim%ComputeRunTime()
127 WRITE(testinfo(k),'(A,I3,A38,A,I2,A11,F5.2,A2)') "#", k, &
128 ". ODE solver: " // adjustl(sim%Timedisc%GetName()) // ", ", &
129 "order", sim%Timedisc%GetOrder(), ", runtime: ", sim%run_time, " s"
130 ELSE
131 WRITE(testinfo(k),'(A)') "#", k, ". ODE solver: " // adjustl(sim%Timedisc%GetName()) // "FAILED"
132 END IF
133
134 IF(k.LT.num_tests) THEN
135 CALL sim%Finalize(.false.)
136 DEALLOCATE(sim)
137 END IF
138 END DO
139
140 ! print test summary
141 CALL sim%Finalize()
142 CALL sim%Info(repeat("*",67))
143 DO k=1,num_tests
144 CALL sim%Info(testinfo(k))
145 END DO
146 CALL sim%Info(repeat("*",67))
147
148 DEALLOCATE(sim)
149
150 tap_done
151
152CONTAINS
153
154 SUBROUTINE makeconfig(Sim, config)
155 IMPLICIT NONE
156 !------------------------------------------------------------------------!
157 CLASS(fosite) :: Sim
158 TYPE(dict_typ),POINTER :: config
159 !------------------------------------------------------------------------!
160 ! Local variable declaration
161 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
162 timedisc, fluxes
163 !------------------------------------------------------------------------!
164 INTENT(INOUT) :: sim
165 !------------------------------------------------------------------------!
166 ! mesh settings
167 mesh => dict("meshtype" / midpoint, &
168 "geometry" / cartesian, &
169 "inum" / res, &
170 "jnum" / res, &
171 "knum" / 1, &
172 "xmin" / (-0.3), &
173 "xmax" / 0.3, &
174 "ymin" / (-0.3), &
175 "ymax" / 0.3, &
176 "zmin" / 0.0, &
177 "zmax" / 0.0, &
178 "gparam" / 1.0)
179
180 ! boundary conditions
181 boundary => dict("western" / no_gradients, &
182 "eastern" / no_gradients, &
183 "southern" / no_gradients, &
184 "northern" / no_gradients, &
185 "bottomer" / no_gradients, &
186 "topper" / no_gradients)
187
188 ! physics settings
189 physics => dict("problem" / euler, &
190 "gamma" / gamma) ! ratio of specific heats !
191
192 ! flux calculation and reconstruction method
193 fluxes => dict("order" / linear, &
194 "fluxtype" / kt, &
195 "variables" / conservative, & ! vars. to use for reconstruction!
196 "limiter" / monocent, & ! one of: minmod, monocent,... !
197 "theta" / 1.2) ! optional parameter for limiter !
198
199 ! time discretization settings
200 timedisc => dict( &
201 "method" / modified_euler, & ! overwritten in main program
202 "cfl" / 0.3, &
203 "ShowButcherTableau" / 1, &
204 "stoptime" / tsim, &
205 "dtlimit" / 1.0e-15, &
206 "tol_rel" / 0.01, &
207 "tol_abs" / (/0.0,1e-3,1e-3,0.0/), &
208 "output/error" / 1, &
209 "maxiter" / 1000000)
210
211 ! initialize data input/output
212 datafile => dict("fileformat" / vtk, &
213! datafile => Dict("fileformat" / GNUPLOT, "filecycles" / 0, &
214 "filename" / (trim(odir) // trim(ofname)), & ! overwritte in main program
215 "count" / onum)
216
217 config => dict("mesh" / mesh, &
218 "physics" / physics, &
219 "boundary" / boundary, &
220 "fluxes" / fluxes, &
221 "timedisc" / timedisc, &
222 "datafile" / datafile)
223 END SUBROUTINE makeconfig
224
225
226 SUBROUTINE initdata(Mesh,Physics,Timedisc)
227 USE physics_euler_mod, ONLY : physics_euler
228 IMPLICIT NONE
229 !------------------------------------------------------------------------!
230 CLASS(physics_base), INTENT(IN) :: Physics
231 CLASS(mesh_base), INTENT(IN) :: Mesh
232 CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
233 !------------------------------------------------------------------------!
234 ! Local variable declaration
235 INTEGER :: n
236 REAL :: P1
237 !------------------------------------------------------------------------!
238 ! isothermal modules are excluded
239 SELECT TYPE (phys => physics)
240 CLASS IS(physics_euler)
241 SELECT TYPE (pvar => timedisc%pvar)
242 CLASS IS(statevector_euler)
243 ! uniform density everywhere
244 pvar%density%data1d(:) = rho0
245 ! vanishing initial velocities
246 pvar%velocity%data1d(:) = 0.0
247 ! set initial peak pressure P1 inside sphere with radius R0 centered on the origin
248 n = 2 ! 2 for 2D
249 p1 = 3.*(phys%gamma - 1.0)*e1 / ((n + 1)*pi*r0**n)
250 WHERE (mesh%radius%bcenter(:,:,:).LE.r0)
251 ! behind the shock front
252 pvar%pressure%data3d(:,:,:) = p1
253 ELSEWHERE
254 ! in front of the shock front (ambient medium)
255 pvar%pressure%data3d(:,:,:) = p0
256 END WHERE
257 CLASS DEFAULT
258 ! abort
259 CALL phys%Error("ode:InitData","statevector must be of class euler")
260 END SELECT
261 CLASS DEFAULT
262 ! abort
263 CALL phys%Error("ode:InitData","physics not supported")
264 END SELECT
265
266 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
267 CALL mesh%Info(" DATA-----> initial condition: 2D Sedov explosion")
268
269 END SUBROUTINE initdata
270
271
272END PROGRAM ode_test
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
makeconfig
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
fosite_mod
Definition: fosite.f90:41
physics_euler_mod
physics module for 1D,2D and 3D non-isothermal Euler equations
Definition: physics_euler.f90:41
ode_test
program ode_test
Definition: ode.f90:38
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
physics_euler_mod::physics_euler
Definition: physics_euler.f90:54
physics_euler_mod::statevector_euler
Definition: physics_euler.f90:105
Generated by   doxygen 1.9.4