KHI2D_sphere.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: KHI2D_sphere.f90 #
5!# #
6!# Copyright (C) 2006-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Jubin Lirawi <jlirawi@astrophysik.uni-kiel.de> #
9!# Lars Boesch <lboesch@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!----------------------------------------------------------------------------!
37!----------------------------------------------------------------------------!
38PROGRAM khi
39 USE fosite_mod
40#include "tap.h"
41 IMPLICIT NONE
42 !--------------------------------------------------------------------------!
43 ! simulation parameters ! !
44 REAL, PARAMETER :: tsim = 10.0 ! simulation time !
45 REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats !
46 REAL, PARAMETER :: re = 1.0e+4 ! Reynolds number (HUGE(RE) disables viscosity)
47 REAL, PARAMETER :: omega = 0.1 ! angular speed of rotational frame
48 ! initial condition ! !
49 REAL, PARAMETER :: rho0 = 1.0 ! outer region: density !
50 REAL, PARAMETER :: v0 = 0.0 ! velocity !
51 REAL, PARAMETER :: p0 = 2.5 ! pressure !
52 REAL, PARAMETER :: rho1 = 2.0 ! inner region: density !
53 REAL, PARAMETER :: v1 = 1.0 ! velocity !
54 REAL, PARAMETER :: p1 = p0 ! pressure !
55 ! mesh settings ! !
56 INTEGER, PARAMETER :: mgeo = spherical_planet ! geometry of the mesh !
57 INTEGER, PARAMETER :: xres = 20 ! x-resolution !
58 INTEGER, PARAMETER :: yres = 40 ! y-resolution !
59 REAL, PARAMETER :: gpar = 1.0 ! radius of sphere !
60 ! output file parameter ! !
61 INTEGER, PARAMETER :: onum = 10 ! number of output data sets !
62 CHARACTER(LEN=256), PARAMETER & ! output data dir !
63 :: odir = './' ! !
64 CHARACTER(LEN=256), PARAMETER & ! output data file name !
65 :: ofname = 'KHI2D_sphere' ! !
66 !--------------------------------------------------------------------------!
67 CLASS(fosite), ALLOCATABLE :: sim
68 LOGICAL :: ok
69 !--------------------------------------------------------------------------!
70 tap_plan(1)
71
72 ALLOCATE(sim)
73 CALL sim%InitFosite()
74 CALL makeconfig(sim,sim%config)
75 CALL sim%Setup()
76 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
77
78 CALL sim%Run()
79 ok = .NOT.sim%aborted
80 CALL sim%Finalize()
81 DEALLOCATE(sim)
82
83 tap_check(ok, "stoptime reached")
84 tap_done
85
86CONTAINS
87
88 SUBROUTINE makeconfig(Sim, config)
89 IMPLICIT NONE
90 !------------------------------------------------------------------------!
91 CLASS(fosite),INTENT(INOUT) :: Sim
92 TYPE(dict_typ),POINTER :: config
93 !------------------------------------------------------------------------!
94 ! Local variable declaration
95 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
96 timedisc, fluxes, vis, rotframe
97 REAL :: dynvis
98 !------------------------------------------------------------------------!
99 ! mesh settings
100 mesh => dict( &
101 "meshtype" / midpoint, &
102 "geometry" / mgeo, &
103 "fargo/method" / 0, &
104 "decomposition" / (/ -1, 1, 1/), &
105 "inum" / xres, & ! resolution in x, !
106 "jnum" / yres, & ! y and !
107 "knum" / 1, & ! z direction !
108 "xmin" / (0.1), &
109 "xmax" / (pi-0.1), &
110 "ymin" / (0.0), &
111 "ymax" / (2.0*pi), &
112 "zmin" / gpar, &
113 "zmax" / gpar, &
114! "zmin" / GPAR, &
115! "zmax" / (GPAR + .1), &
116 "gparam" / gpar &
117 )
118
119 ! physics settings
120 physics => dict( &
121 "problem" / euler, &
122 "gamma" / gamma & ! ratio of specific heats !
123 )
124
125 ! flux calculation and reconstruction method
126 fluxes => dict( &
127 "fluxtype" / kt, &
128 "order" / linear, &
129 "variables" / conservative, & ! vars. to use for reconstruction !
130! "variables" / PRIMITIVE, & ! vars. to use for reconstruction !
131 "limiter" / monocent, & ! one of: minmod, monocent,... !
132! "output/pfluxes" / 1, &
133! "output/pstates" / 1, &
134! "output/cstates" / 1, &
135 "theta" / 1.2 & ! optional parameter for limiter !
136 )
137
138 ! boundary conditions
139 boundary => dict( &
140 "western" / reflecting, &
141 "eastern" / reflecting, &
142 "southern" / periodic, &
143 "northern" / periodic, &
144 "bottomer" / reflecting, &
145 "topper" / reflecting &
146 )
147
148 NULLIFY(vis,rotframe)
149 ! viscosity source term
150 ! compute dynamic viscosity constant using typical scales and Reynolds number
151 dynvis = abs(rho0 * 2*pi * (v0-v1) / re)
152 IF (re.LT.huge(re)) THEN
153 vis => dict( &
154 "stype" / viscosity, &
155 "vismodel" / molecular, &
156 "dynconst" / dynvis, &
157 "bulkconst" / (-2./3.*dynvis), &
158 "output/dynvis" / 0, &
159 "output/stress" / 0, &
160 "output/kinvis" / 0, &
161 "output/bulkvis" / 0 &
162 )
163 END IF
164 ! activate inertial forces due to rotating frame if OMEGA > 0
165 IF (omega.GT.tiny(omega)) THEN
166 rotframe => dict("stype" / rotating_frame, &
167 "disable_centaccel" / 1)
168 CALL setattr(mesh,"omega",omega)
169 END IF
170
171 ! time discretization settings
172 timedisc => dict( &
173 "method" / modified_euler, &
174! "method" / SSPRK, &
175 "order" / 3, &
176 "cfl" / 0.4, &
177 "stoptime" / tsim, &
178 "tol_rel" / 0.0095, &
179 "dtlimit" / 1.0e-15, &
180 "maxiter" / 100000000, &
181 "output/pressure" / 1, &
182 "output/density" / 1, &
183 "output/xvelocity" / 1, &
184 "output/yvelocity" / 1 &
185 )
186
187 ! initialize data input/output
188 datafile => dict(&
189! "fileformat" / VTK, &
190 "fileformat" / xdmf, &
191 "filename" / (trim(odir) // trim(ofname)), &
192 "count" / onum &
193 )
194
195 config => dict( &
196 "mesh" / mesh, &
197 "physics" / physics, &
198 "boundary" / boundary, &
199 "fluxes" / fluxes, &
200 "timedisc" / timedisc, &
201 "datafile" / datafile &
202 )
203 IF (ASSOCIATED(vis)) &
204 CALL setattr(config, "sources/vis", vis)
205 IF (ASSOCIATED(rotframe)) &
206 CALL setattr(config, "sources/rotframe", rotframe)
207 END SUBROUTINE makeconfig
208
209 SUBROUTINE initdata(Mesh,Physics,Timedisc)
210 USE sources_base_mod, ONLY : sources_base
211 USE sources_rotframe_mod, ONLY : sources_rotframe
212 IMPLICIT NONE
213 !------------------------------------------------------------------------!
214 CLASS(physics_base) :: Physics
215 CLASS(mesh_base) :: Mesh
216 CLASS(timedisc_base) :: Timedisc
217 !------------------------------------------------------------------------!
218 ! Local variable declaration
219 CLASS(sources_base), POINTER :: sp => null()
220 INTEGER :: k,n,clock
221 INTEGER, DIMENSION(:), ALLOCATABLE :: seed
222 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,2) :: dv
223 !------------------------------------------------------------------------!
224 INTENT(IN) :: mesh,physics
225 INTENT(INOUT) :: timedisc
226 !------------------------------------------------------------------------!
227 ! Seed the random number generator
228 ! initialize with a mix from current time and mpi rank
229 CALL system_clock(count=clock)
230 CALL random_seed(SIZE = n) ! determine the minimum size of the seed
231 seed = clock + (timedisc%getrank()+1) * (/(k-1, k=1,n)/)
232 CALL random_seed(put=seed)
233 CALL random_number(dv)
234
235 ! initial condition
236 SELECT TYPE(pvar => timedisc%pvar)
237 TYPE IS(statevector_euler) ! non-isothermal HD
238 ! flow along x-direction:
239 WHERE ((mesh%bcenter(:,:,:,1).LT.(mesh%ymin+0.33*pi)).OR. &
240 (mesh%bcenter(:,:,:,1).GT.(sim%Mesh%ymin+0.66*pi)))
241 pvar%density%data3d(:,:,:) = rho0
242 pvar%velocity%data4d(:,:,:,2) = v0
243 pvar%velocity%data4d(:,:,:,1) = 0.0
244 pvar%pressure%data3d(:,:,:) = p0
245 ELSEWHERE
246 pvar%density%data3d(:,:,:) = rho1
247 pvar%velocity%data4d(:,:,:,2) = v1
248 pvar%velocity%data4d(:,:,:,1) = 0.0
249 pvar%pressure%data3d(:,:,:) = p1
250 END WHERE
251 sp => sim%Sources%GetSourcesPointer(rotating_frame)
252 IF (ASSOCIATED(sp)) THEN
253 SELECT TYPE(sp)
254 CLASS IS(sources_rotframe)
255 CALL sp%Convert2RotatingFrame(mesh,physics,pvar)
256 END SELECT
257 END IF
258
259 SELECT CASE(mesh%fargo%GetType())
260 CASE(1,2)
261 ! 2D simulation -> 2nd velocity is the azimuthal velocity
262 timedisc%w(:,:) = pvar%velocity%data4d(:,mesh%JMIN,:,2)
263 END SELECT
264
265 ! add velocity perturbations
266 pvar%velocity%data4d(:,:,:,1) = pvar%velocity%data4d(:,:,:,1) &
267 + (dv(:,:,:,1)-0.5)*0.2
268 pvar%velocity%data4d(:,:,:,2) = pvar%velocity%data4d(:,:,:,2) &
269 + (dv(:,:,:,2)-0.5)*0.2
270 CLASS DEFAULT
271 CALL physics%Error("KHI2D_sphere::InitData","only non-isothermal HD supported")
272 END SELECT
273
274 ! initial condition
275 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
276 CALL mesh%Info(" DATA-----> initial condition: " // &
277 "spherical shear flow")
278 END SUBROUTINE initdata
279
280END PROGRAM khi
khi
program khi
Kelvin-Helmholtz instability.
Definition: KHI2D.f90:37
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
sources_base_mod
generic source terms module providing functionaly common to all source terms
Definition: sources_base.f90:42
sources_rotframe_mod
source terms module for inertial forces caused by a rotating grid
Definition: sources_rotframe.f90:39
omega
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
Definition: sources_viscosity.f90:432
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
sources_base_mod::sources_base
Definition: sources_base.f90:57
sources_rotframe_mod::sources_rotframe
Definition: sources_rotframe.f90:52
Generated by   doxygen 1.9.4