KHI2D.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: KHI2D.f90 #
5!# #
6!# Copyright (C) 2006-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Jubin Lirawi <jlirawi@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!----------------------------------------------------------------------------!
36!----------------------------------------------------------------------------!
37PROGRAM khi
38 USE fosite_mod
39#include "tap.h"
40 IMPLICIT NONE
41 !--------------------------------------------------------------------------!
42 ! simulation parameters
43 REAL, PARAMETER :: tsim = 10.0 ! simulation time
44 REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats
45 REAL, PARAMETER :: re = 1.0e+4 ! Reynolds number (HUGE(RE) disables viscosity)
46! REAL, PARAMETER :: RE = HUGE(RE)
47 ! initial condition
48 REAL, PARAMETER :: rho0 = 1.0 ! outer region: density
49 REAL, PARAMETER :: v0 =-0.5 ! velocity
50 REAL, PARAMETER :: p0 = 2.5 ! pressure
51 REAL, PARAMETER :: rho1 = 2.0 ! inner region: density
52 REAL, PARAMETER :: v1 = 0.5 ! x-velocity
53 REAL, PARAMETER :: p1 = p0 ! pressure
54 ! mesh settings
55 INTEGER, PARAMETER :: mgeo = cartesian ! geometry of the mesh
56 INTEGER, PARAMETER :: res = 32 ! resolution
57 REAL, PARAMETER :: xyzlen= 1.0 ! spatial extend
58 ! output file parameter
59 INTEGER, PARAMETER :: onum = 10 ! number of output data sets
60 CHARACTER(LEN=256), PARAMETER & ! output data dir
61 :: odir = './'
62 CHARACTER(LEN=256), PARAMETER & ! output data file name
63 :: ofname = 'KHI2D'
64 !--------------------------------------------------------------------------!
65 CLASS(fosite), ALLOCATABLE :: sim
66 CLASS(marray_compound), POINTER :: pvar,pvar_init
67 INTEGER :: n,fargo,d,dt(8)
68 REAL :: sigma(2,3)
69 INTEGER, DIMENSION(:), ALLOCATABLE :: seed
70 CHARACTER(LEN=2),PARAMETER :: dir_name(2) = (/"xy","xz"/)
71 REAL, DIMENSION(:,:), ALLOCATABLE :: w
72 !--------------------------------------------------------------------------!
73 tap_plan(6)
74
75 ! allocate memory for dynamic types and arrays
76 CALL random_seed(size=n)
77 ALLOCATE(seed(n))
78 ! Seed the random number generator using current system date and time
79 call date_and_time(values=dt)
80 seed = sum(dt(1:8)) * (/(d-1, d=1,n)/)
81 DO d=1,2
82 DO fargo=0,2
83 ALLOCATE(sim,pvar,pvar_init)
84 ! simulate KHI with initial x-velocity along x-direction
85 CALL sim%InitFosite()
86 CALL makeconfig(sim,sim%config,dir_name(d))
87 CALL sim%Setup()
88 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, dir_name(d))
89 IF (ASSOCIATED(sim%Timedisc%w).AND..NOT.ALLOCATED(w)) ALLOCATE(w,source=transpose(sim%Timedisc%w))
90 ! store transposed initial data
91 CALL sim%Physics%new_statevector(pvar_init,primitive)
92 CALL transposedata(sim%Mesh,sim%Timedisc%pvar,pvar_init,dir_name(d))
93 CALL sim%Run()
94 ! store transposed result of the first run
95 CALL sim%Physics%new_statevector(pvar,primitive)
96 CALL transposedata(sim%Mesh,sim%Timedisc%pvar,pvar,dir_name(d))
97 ! finish the simulation
98 CALL sim%Finalize(mpifinalize_=.false.)
99 DEALLOCATE(sim)
100
101 ! simulate KHI with transposed initial data from the previous run
102 ALLOCATE(sim)
103 CALL sim%InitFosite()
104 CALL makeconfig(sim, sim%config,dir_name(d))
105 CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "-" // dir_name(d) &
106 // "-fargo" // achar(48+fargo) // "_transposed"))
107 CALL setattr(sim%config, "/mesh/fargo/direction", (1+d))
108#ifdef PARALLEL
109 SELECT CASE(dir_name(d))
110 CASE("xy")
111 ! decompose computational domain along y-direction only
112 CALL setattr(sim%config, "/mesh/decomposition", (/1,-1,1/))
113 CASE("xz")
114 ! decompose computational domain along z-direction only
115 CALL setattr(sim%config, "/mesh/decomposition", (/1,1,-1/))
116 END SELECT
117#endif
118 CALL sim%Setup()
119 sim%Timedisc%pvar%data1d(:) = pvar_init%data1d(:)
120 IF (ASSOCIATED(sim%Timedisc%w)) THEN
121 sim%Timedisc%w(:,:) = reshape(transpose(w(:,:)),shape(sim%Timedisc%w(:,:)))
122 END IF
123 SELECT TYPE(phys => sim%Physics)
124 TYPE IS(physics_euler)
125 CALL phys%Convert2Conservative(sim%Timedisc%pvar,sim%Timedisc%cvar)
126 CLASS DEFAULT
127 CALL phys%Error("KHI2D::InitData","only non-isothermal HD supported")
128 END SELECT
129 CALL sim%Run()
130 ! compare results
131 ! IMPORTANT: Use arrays with collapsed dimensions here! Since pvar has the shape
132 ! from the previous run (with already transposed data) the comparison would fail if the
133 ! resolution in these dimensions differ. This is of particular importance for parallel tests.
134 sigma(d,fargo+1) = sqrt(sum((sim%Timedisc%pvar%data1d(:)-pvar%data1d(:))**2)/SIZE(pvar%data1d))
135
136 IF (d.EQ.2.AND.fargo.EQ.2) THEN
137 CALL sim%Finalize(mpifinalize_=.true.)
138 ELSE
139 CALL sim%Finalize(mpifinalize_=.false.)
140 END IF
141 DEALLOCATE(pvar,pvar_init,sim)
142 IF (ALLOCATED(w)) DEALLOCATE(w)
143 END DO
144 END DO
145 DEALLOCATE(seed)
146 DO d=1,2
147 tap_check_small(sigma(d,1),tiny(sigma),dir_name(d) // "-symmetry test (w/o fargo)")
148 tap_check_small(sigma(d,2),tiny(sigma),dir_name(d) // "-symmetry test (fargo type 1)")
149 tap_check_small(sigma(d,3),tiny(sigma),dir_name(d) // "-symmetry test (fargo type 2)")
150 END DO
151 tap_done
152
153CONTAINS
154
155 SUBROUTINE makeconfig(Sim, config, dir)
156 IMPLICIT NONE
157 !------------------------------------------------------------------------!
158 CLASS(fosite),INTENT(INOUT) :: Sim
159 TYPE(dict_typ),POINTER :: config
160 CHARACTER(LEN=2), INTENT(IN) :: dir
161 !------------------------------------------------------------------------!
162 ! Local variable declaration
163 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
164 sources, timedisc, fluxes
165 REAL :: dynvis
166 !------------------------------------------------------------------------!
167 ! mesh settings
168 SELECT CASE(dir)
169 CASE("xy")
170 mesh => dict( &
171 "meshtype" / midpoint, &
172 "geometry" / mgeo, &
173 "fargo/method" / fargo, &
174 "fargo/direction" / 1, &
175 "decomposition" / (/-1,1,1/), &
176 "inum" / res, & ! resolution in x, !
177 "jnum" / res, & ! y and !
178 "knum" / 1, & ! z direction !
179 "xmin" / (-0.5*xyzlen), &
180 "xmax" / (0.5*xyzlen), &
181 "ymin" / (-0.5*xyzlen), &
182 "ymax" / (0.5*xyzlen), &
183 "zmin" / (-0.0), &
184 "zmax" / (0.0) )
185
186 CASE("xz")
187 mesh => dict( &
188 "meshtype" / midpoint, &
189 "geometry" / mgeo, &
190 "fargo/method" / fargo, &
191 "fargo/direction" / 1, &
192 "decomposition" / (/-1,1,1/), &
193 "inum" / res, & ! resolution in x, !
194 "knum" / res, & ! y and !
195 "jnum" / 1, & ! z direction !
196 "xmin" / (-0.5*xyzlen), &
197 "xmax" / (0.5*xyzlen), &
198 "zmin" / (-0.5*xyzlen), &
199 "zmax" / (0.5*xyzlen), &
200 "ymin" / (-0.0), &
201 "ymax" / (0.0) )
202 CASE DEFAULT
203 CALL sim%Error("KHI2D::InitData","only xy- and xz-symmetry test supported")
204 END SELECT
205
206 ! physics settings
207 physics => dict( &
208 "problem" / euler, &
209 "gamma" / gamma & ! ratio of specific heats !
210 )
211
212 ! flux calculation and reconstruction method
213 fluxes => dict( &
214 "fluxtype" / kt, &
215 "order" / linear, &
216 "variables" / conservative, & ! vars. to use for reconstruction !
217! "variables" / PRIMITIVE, & ! vars. to use for reconstruction !
218 "limiter" / monocent, & ! one of: minmod, monocent,... !
219! "output/pfluxes" / 1, &
220! "output/pstates" / 1, &
221! "output/cstates" / 1, &
222 "theta" / 1.2 & ! optional parameter for limiter !
223 )
224
225 ! boundary conditions
226 boundary => dict( &
227 "western" / periodic, &
228 "eastern" / periodic, &
229 "southern" / periodic, &
230 "northern" / periodic, &
231 "bottomer" / periodic, &
232 "topper" / periodic &
233 )
234
235 NULLIFY(sources)
236 ! viscosity source term
237 ! compute dynamic viscosity constant using typical scales and Reynolds number
238 dynvis = abs(rho0 * xyzlen * (v0-v1) / re)
239 IF (dynvis.GT.tiny(1.0)) THEN
240 sources => dict( &
241 "vis/stype" / viscosity, &
242 "vis/vismodel" / molecular, &
243 "vis/dynconst" / dynvis, &
244 "vis/bulkconst" / (-2./3.*dynvis), &
245 "vis/output/dynvis" / 0, &
246 "vis/output/stress" / 0, &
247 "vis/output/kinvis" / 0, &
248 "vis/output/bulkvis" / 0 &
249 )
250 END IF
251
252 ! time discretization settings
253 timedisc => dict( &
254 "method" / modified_euler, &
255 "order" / 3, &
256 "cfl" / 0.4, &
257 "stoptime" / tsim, &
258 "dtlimit" / 1.0e-10, &
259 "maxiter" / 100000, &
260! "output/geometrical_sources" / 1, &
261 "output/energy" / 1 &
262 )
263
264 ! initialize data input/output
265 datafile => dict(&
266 "fileformat" / vtk, &
267! "fileformat" / XDMF, &
268 "filename" / (trim(odir) // trim(ofname) // "-" // dir // "-fargo" // achar(48+fargo)), &
269 "count" / onum &
270 )
271
272 config => dict( &
273 "mesh" / mesh, &
274 "physics" / physics, &
275 "boundary" / boundary, &
276 "fluxes" / fluxes, &
277 "timedisc" / timedisc, &
278 "datafile" / datafile &
279 )
280 IF (ASSOCIATED(sources)) &
281 CALL setattr(config, "sources", sources)
282 END SUBROUTINE makeconfig
283
284
285 SUBROUTINE initdata(Mesh,Physics,Timedisc,dir)
286 USE physics_euler_mod, ONLY : statevector_euler
287 USE physics_eulerisotherm_mod, ONLY : statevector_eulerisotherm
288 IMPLICIT NONE
289 !------------------------------------------------------------------------!
290 CLASS(physics_base) :: Physics
291 CLASS(mesh_base) :: Mesh
292 CLASS(timedisc_base) :: Timedisc
293 CHARACTER(LEN=2) :: dir
294 !------------------------------------------------------------------------!
295 ! Local variable declaration
296 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,2) :: dv
297 !------------------------------------------------------------------------!
298 INTENT(IN) :: mesh,physics,dir
299 INTENT(INOUT) :: timedisc
300 !------------------------------------------------------------------------!
301 CALL random_seed(put=seed)
302 CALL random_number(dv)
303
304 ! initial condition
305 SELECT TYPE(pvar => timedisc%pvar)
306 TYPE IS(statevector_euler) ! non-isothermal HD
307 SELECT CASE(dir)
308 CASE("xy")
309 ! flow along x-direction extends in y-direction
310 WHERE ((mesh%bcenter(:,:,:,2).LT.(mesh%ymin+0.25*xyzlen)).OR. &
311 (mesh%bcenter(:,:,:,2).GT.(sim%Mesh%ymin+0.75*xyzlen)))
312 pvar%density%data3d(:,:,:) = rho0
313 pvar%velocity%data4d(:,:,:,1) = v0
314 pvar%velocity%data4d(:,:,:,2) = 0.0
315 pvar%pressure%data3d(:,:,:) = p0
316 ELSEWHERE
317 pvar%density%data3d(:,:,:) = rho1
318 pvar%velocity%data4d(:,:,:,1) = v1
319 pvar%velocity%data4d(:,:,:,2) = 0.0
320 pvar%pressure%data3d(:,:,:) = p1
321 END WHERE
322 ! set fargo background velocity to initial velocity along
323 ! x-direction
324 SELECT CASE(mesh%fargo%GetType())
325 CASE(1,2)
326 WHERE ((mesh%bcenter(mesh%IMIN,:,:,2).LT.(mesh%ymin+0.25*xyzlen)).OR. &
327 (mesh%bcenter(mesh%IMIN,:,:,2).GT.(sim%Mesh%ymin+0.75*xyzlen)))
328 timedisc%w(:,:) = v0
329 ELSEWHERE
330 timedisc%w(:,:) = v1
331 END WHERE
332 END SELECT
333 CASE("xz")
334 ! flow along x-direction extends in z-direction
335 WHERE ((mesh%bcenter(:,:,:,3).LT.(mesh%zmin+0.25*xyzlen)).OR. &
336 (mesh%bcenter(:,:,:,3).GT.(sim%Mesh%zmin+0.75*xyzlen)))
337 pvar%density%data3d(:,:,:) = rho0
338 pvar%velocity%data4d(:,:,:,1) = v0
339 pvar%velocity%data4d(:,:,:,2) = 0.0
340 pvar%pressure%data3d(:,:,:) = p0
341 ELSEWHERE
342 pvar%density%data3d(:,:,:) = rho1
343 pvar%velocity%data4d(:,:,:,1) = v1
344 pvar%velocity%data4d(:,:,:,2) = 0.0
345 pvar%pressure%data3d(:,:,:) = p1
346 END WHERE
347 ! set fargo background velocity to initial velocity along
348 ! x-direction
349 SELECT CASE(mesh%fargo%GetType())
350 CASE(1,2)
351 WHERE ((mesh%bcenter(mesh%IMIN,:,:,3).LT.(mesh%zmin+0.25*xyzlen)).OR. &
352 (mesh%bcenter(mesh%IMIN,:,:,3).GT.(sim%Mesh%zmin+0.75*xyzlen)))
353 timedisc%w(:,:) = v0
354 ELSEWHERE
355 timedisc%w(:,:) = v1
356 END WHERE
357 END SELECT
358 END SELECT
359 ! add perturbations
360 pvar%velocity%data4d(:,:,:,1:2) = pvar%velocity%data4d(:,:,:,1:2) &
361 + (dv(:,:,:,1:2)-0.5)*0.02
362 SELECT TYPE(cvar => timedisc%cvar)
363 TYPE IS(statevector_euler)
364 SELECT TYPE(phys => physics)
365 TYPE IS(physics_euler)
366 CALL phys%Convert2Conservative(pvar,cvar)
367 CALL mesh%Info(" DATA-----> initial condition: Kelvin-Helmholtz instability")
368 END SELECT
369 CLASS DEFAULT
370 CALL physics%Error("KHI2D::InitData","only non-isothermal HD supported")
371 END SELECT
372 CLASS DEFAULT
373 CALL physics%Error("KHI2D::InitData","only non-isothermal HD supported")
374 END SELECT
375 END SUBROUTINE initdata
376
377 SUBROUTINE transposedata(Mesh,pvar_in,pvar_out,dir)
378 IMPLICIT NONE
379 !------------------------------------------------------------------------!
380 CLASS(mesh_base), INTENT(IN) :: Mesh
381 CLASS(marray_compound), INTENT(INOUT) :: pvar_in,pvar_out
382 CHARACTER(LEN=2), INTENT(IN) :: dir
383! REAL,DIMENSION(1:SIZE(pvar_in%data3d,DIM=1)),TARGET :: data1d
384 REAL,ALLOCATABLE,TARGET :: data1d(:)
385! REAL,DIMENSION(1:SIZE(pvar_in%data4d,DIM=1),1:SIZE(pvar_in%data4d,DIM=3)),TARGET :: data2d
386 REAL,POINTER,CONTIGUOUS :: tpdata2d(:,:),tpdata3d(:,:,:)
387 !------------------------------------------------------------------------!
388 INTEGER :: j,k
389 !------------------------------------------------------------------------!
390 SELECT TYPE(pin => pvar_in)
391 TYPE IS(statevector_euler) ! non-isothermal HD
392 SELECT TYPE(pout => pvar_out)
393 TYPE IS(statevector_euler) ! non-isothermal HD
394 SELECT CASE(dir)
395 CASE("xy")
396 ! transpose x-y directions
397 DO k=mesh%KGMIN,mesh%KGMAX
398 pout%density%data3d(:,:,k) = reshape(transpose(pin%density%data3d(:,:,k)), &
399 shape(pout%density%data3d(:,:,k)))
400 pout%pressure%data3d(:,:,k) = reshape(transpose(pin%pressure%data3d(:,:,k)), &
401 shape(pout%pressure%data3d(:,:,k)))
402 ! ... and exchange x- and z-velocities
403 pout%velocity%data4d(:,:,k,1) = reshape(transpose(pin%velocity%data4d(:,:,k,2)), &
404 shape(pout%velocity%data4d(:,:,k,1)))
405 pout%velocity%data4d(:,:,k,2) = reshape(transpose(pin%velocity%data4d(:,:,k,1)), &
406 shape(pout%velocity%data4d(:,:,k,2)))
407 END DO
408 CASE("xz")
409 ! transpose x-z directions
410 DO j=mesh%JGMIN,mesh%JGMAX
411 ! simply transpose the 1st and 3rd indices ...
412 pout%density%data3d(:,j,:) = reshape(transpose(pin%density%data3d(:,j,:)), &
413 shape(pout%density%data3d(:,j,:)))
414 pout%pressure%data3d(:,j,:) = reshape(transpose(pin%pressure%data3d(:,j,:)), &
415 shape(pout%pressure%data3d(:,j,:)))
416 ! ... and exchange x- and z-velocities
417 pout%velocity%data4d(:,j,:,1) = reshape(transpose(pin%velocity%data4d(:,j,:,2)), &
418 shape(pout%velocity%data4d(:,j,:,1)))
419 pout%velocity%data4d(:,j,:,2) = reshape(transpose(pin%velocity%data4d(:,j,:,1)), &
420 shape(pout%velocity%data4d(:,j,:,2)))
421 END DO
422 CASE DEFAULT
423 CALL mesh%Error("KHI2D::TransposeData","only transposition of xy- and xz-direction supported")
424 END SELECT
425 END SELECT
426 END SELECT
427 END SUBROUTINE transposedata
428
429END PROGRAM khi
khi
program khi
Kelvin-Helmholtz instability.
Definition: KHI2D.f90:37
transposedata
subroutine transposedata(Mesh, pvar_in, pvar_out, dir)
Definition: KHI2D.f90:378
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
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
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
physics_eulerisotherm_mod::statevector_eulerisotherm
Definition: physics_eulerisotherm.f90:122
Generated by   doxygen 1.9.4