KHI3D.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: KHI.f90 #
5 !# #
6 !# Copyright (C) 2006-2012 #
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 !----------------------------------------------------------------------------!
32 
33 !! References:
34 !! \cite thomson1871
35 !! \cite helmholtz1868
36 !----------------------------------------------------------------------------!
37 PROGRAM khi
38  USE fosite_mod
39 #include "tap.h"
40  IMPLICIT NONE
41  !--------------------------------------------------------------------------!
42  ! simulation parameters
43  REAL, PARAMETER :: tsim = 30.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  ! initial condition
47  REAL, PARAMETER :: rho0 = 1.0 ! outer region: density
48  REAL, PARAMETER :: v0 =-0.5 ! velocity
49  REAL, PARAMETER :: p0 = 2.5 ! pressure
50  REAL, PARAMETER :: rho1 = 2.0 ! inner region: density
51  REAL, PARAMETER :: v1 = 0.5 ! x-velocity
52  REAL, PARAMETER :: p1 = p0 ! pressure
53  ! mesh settings
54  INTEGER, PARAMETER :: mgeo = cartesian ! geometry of the mesh
55  INTEGER, PARAMETER :: res = 10 ! resolution
56  REAL, PARAMETER :: xyzlen= 1.0 ! spatial extend
57  ! output file parameter
58  INTEGER, PARAMETER :: onum = 10 ! 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 = 'KHI'
63  !--------------------------------------------------------------------------!
64 ! TYPE(fosite_TYP) :: Sim
65  INTEGER :: n
66  INTEGER :: i,j,k
67  REAL :: sigma_1,sigma_2,sigma_3,sigma_4,sigma_5,sigma_6
68  INTEGER, DIMENSION(:), ALLOCATABLE :: seed
69  REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: &
70  pvar_xy,pvar_xz,pvar_yx,pvar_yz,pvar_zx,pvar_zy, pvar_temp
71  !--------------------------------------------------------------------------!
72  CLASS(fosite), ALLOCATABLE :: sim
73  !--------------------------------------------------------------------------!
74  tap_plan(6)
75 
76  ! allocate memory for random seed variable
77  CALL random_seed(size=n)
78  ALLOCATE(seed(n))
79 
80  ! test flow along x-direction
81  ALLOCATE(sim)
82  CALL sim%InitFosite()
83  CALL makeconfig(sim, sim%config)
84  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_xy"))
85  CALL sim%Setup()
86  ! allocate memory to store the result of the first run
87  ALLOCATE(pvar_xy(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%VNUM))
88  ALLOCATE(pvar_xz(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%VNUM))
89  ALLOCATE(pvar_yx(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%VNUM))
90  ALLOCATE(pvar_yz(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%VNUM))
91  ALLOCATE(pvar_zx(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%VNUM))
92  ALLOCATE(pvar_zy(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%VNUM))
93  ALLOCATE(pvar_temp(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%VNUM))
94  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,'xy',.false.)
95  CALL sim%Run()
96  pvar_xy(:,:,:,:) = sim%Timedisc%pvar%data4d(:,:,:,:)
97  CALL sim%Finalize(.false.)
98  DEALLOCATE(sim)
99 
100  ! test flow along y-direction
101  ALLOCATE(sim)
102  CALL sim%InitFosite()
103  CALL makeconfig(sim, sim%config)
104  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_yx"))
105  CALL sim%Setup()
106  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,'yx',.true.)
107  CALL sim%Run()
108  pvar_yx(:,:,:,:) = sim%Timedisc%pvar%data4d(:,:,:,:)
109  CALL sim%Finalize(.false.)
110  DEALLOCATE(sim)
111 
112  ! test flow along z-direction
113  ALLOCATE(sim)
114  CALL sim%InitFosite()
115  CALL makeconfig(sim, sim%config)
116  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_xz"))
117  CALL sim%Setup()
118  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,'xz',.true.)
119  CALL sim%Run()
120  pvar_xz(:,:,:,:) = sim%Timedisc%pvar%data4d(:,:,:,:)
121  CALL sim%Finalize(.false.)
122  DEALLOCATE(sim)
123 
124  ALLOCATE(sim)
125  CALL sim%InitFosite()
126  CALL makeconfig(sim, sim%config)
127  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_zx"))
128  CALL sim%Setup()
129  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,'zx',.true.)
130  CALL sim%Run()
131  pvar_zx(:,:,:,:) = sim%Timedisc%pvar%data4d(:,:,:,:)
132  CALL sim%Finalize(.false.)
133  DEALLOCATE(sim)
134 
135  ALLOCATE(sim)
136  CALL sim%InitFosite()
137  CALL makeconfig(sim, sim%config)
138  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_zy"))
139  CALL sim%Setup()
140  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,'zy',.true.)
141  CALL sim%Run()
142  pvar_zy(:,:,:,:) = sim%Timedisc%pvar%data4d(:,:,:,:)
143  CALL sim%Finalize(.false.)
144  DEALLOCATE(sim)
145 
146  ALLOCATE(sim)
147  CALL sim%InitFosite()
148  CALL makeconfig(sim, sim%config)
149  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_yz"))
150  CALL sim%Setup()
151  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,'yz',.true.)
152  CALL sim%Run()
153  pvar_yz(:,:,:,:) = sim%Timedisc%pvar%data4d(:,:,:,:)
154 
155  DO k=sim%Mesh%KGMIN, sim%Mesh%KGMAX
156  pvar_temp(:,:,k,sim%Physics%DENSITY) = transpose(pvar_xy(:,:,k,sim%Physics%DENSITY))
157  pvar_temp(:,:,k,sim%Physics%XVELOCITY) = transpose(pvar_xy(:,:,k,sim%Physics%YVELOCITY))
158  pvar_temp(:,:,k,sim%Physics%YVELOCITY) = transpose(pvar_xy(:,:,k,sim%Physics%XVELOCITY))
159  pvar_temp(:,:,k,sim%Physics%ZVELOCITY) = transpose(pvar_xy(:,:,k,sim%Physics%ZVELOCITY))
160  pvar_temp(:,:,k,sim%Physics%ENERGY) = transpose(pvar_xy(:,:,k,sim%Physics%ENERGY))
161  END DO
162  sigma_1 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_yx(:,:,:,:))**2)/SIZE(pvar_temp))
163 
164  DO j=sim%Mesh%JGMIN, sim%Mesh%JGMAX
165  pvar_temp(:,j,:,sim%Physics%DENSITY) = transpose(pvar_xz(:,j,:,sim%Physics%DENSITY))
166  pvar_temp(:,j,:,sim%Physics%XVELOCITY) = transpose(pvar_xz(:,j,:,sim%Physics%ZVELOCITY))
167  pvar_temp(:,j,:,sim%Physics%YVELOCITY) = transpose(pvar_xz(:,j,:,sim%Physics%YVELOCITY))
168  pvar_temp(:,j,:,sim%Physics%ZVELOCITY) = transpose(pvar_xz(:,j,:,sim%Physics%XVELOCITY))
169  pvar_temp(:,j,:,sim%Physics%ENERGY) = transpose(pvar_xz(:,j,:,sim%Physics%ENERGY))
170  END DO
171  sigma_2 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_zx(:,:,:,:))**2)/SIZE(pvar_temp))
172 
173  DO i=sim%Mesh%IGMIN, sim%Mesh%IGMAX
174  pvar_temp(i,:,:,sim%Physics%DENSITY) = transpose(pvar_zy(i,:,:,sim%Physics%DENSITY))
175  pvar_temp(i,:,:,sim%Physics%XVELOCITY) = transpose(pvar_zy(i,:,:,sim%Physics%XVELOCITY))
176  pvar_temp(i,:,:,sim%Physics%YVELOCITY) = transpose(pvar_zy(i,:,:,sim%Physics%ZVELOCITY))
177  pvar_temp(i,:,:,sim%Physics%ZVELOCITY) = transpose(pvar_zy(i,:,:,sim%Physics%YVELOCITY))
178  pvar_temp(i,:,:,sim%Physics%ENERGY) = transpose(pvar_zy(i,:,:,sim%Physics%ENERGY))
179  END DO
180  sigma_3 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_yz(:,:,:,:))**2)/SIZE(pvar_temp))
181 
182  DO k=sim%Mesh%KGMIN, sim%Mesh%KGMAX
183  pvar_temp(:,:,k,sim%Physics%DENSITY) = transpose(pvar_zx(:,:,k,sim%Physics%DENSITY))
184  pvar_temp(:,:,k,sim%Physics%XVELOCITY) = transpose(pvar_zx(:,:,k,sim%Physics%YVELOCITY))
185  pvar_temp(:,:,k,sim%Physics%YVELOCITY) = transpose(pvar_zx(:,:,k,sim%Physics%XVELOCITY))
186  pvar_temp(:,:,k,sim%Physics%ZVELOCITY) = transpose(pvar_zx(:,:,k,sim%Physics%ZVELOCITY))
187  pvar_temp(:,:,k,sim%Physics%ENERGY) = transpose(pvar_zx(:,:,k,sim%Physics%ENERGY))
188  END DO
189  sigma_4 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_zy(:,:,:,:))**2)/SIZE(pvar_temp))
190 
191  DO j=sim%Mesh%JGMIN, sim%Mesh%JGMAX
192  pvar_temp(:,j,:,sim%Physics%DENSITY) = transpose(pvar_yx(:,j,:,sim%Physics%DENSITY))
193  pvar_temp(:,j,:,sim%Physics%XVELOCITY) = transpose(pvar_yx(:,j,:,sim%Physics%ZVELOCITY))
194  pvar_temp(:,j,:,sim%Physics%YVELOCITY) = transpose(pvar_yx(:,j,:,sim%Physics%YVELOCITY))
195  pvar_temp(:,j,:,sim%Physics%ZVELOCITY) = transpose(pvar_yx(:,j,:,sim%Physics%XVELOCITY))
196  pvar_temp(:,j,:,sim%Physics%ENERGY) = transpose(pvar_yx(:,j,:,sim%Physics%ENERGY))
197  END DO
198  sigma_5 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_yz(:,:,:,:))**2)/SIZE(pvar_temp))
199 
200  DO i=sim%Mesh%IGMIN, sim%Mesh%IGMAX
201  pvar_temp(i,:,:,sim%Physics%DENSITY) = transpose(pvar_xy(i,:,:,sim%Physics%DENSITY))
202  pvar_temp(i,:,:,sim%Physics%XVELOCITY) = transpose(pvar_xy(i,:,:,sim%Physics%XVELOCITY))
203  pvar_temp(i,:,:,sim%Physics%YVELOCITY) = transpose(pvar_xy(i,:,:,sim%Physics%ZVELOCITY))
204  pvar_temp(i,:,:,sim%Physics%ZVELOCITY) = transpose(pvar_xy(i,:,:,sim%Physics%YVELOCITY))
205  pvar_temp(i,:,:,sim%Physics%ENERGY) = transpose(pvar_xy(i,:,:,sim%Physics%ENERGY))
206  END DO
207  sigma_6 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_xz(:,:,:,:))**2)/SIZE(pvar_temp))
208 
209  CALL sim%Finalize(.true.)
210  DEALLOCATE(sim,pvar_xy,pvar_xz,pvar_yx,pvar_yz,pvar_zx,pvar_zy,pvar_temp,seed)
211  tap_check_small(sigma_1,3*epsilon(sigma_1),"xy-yx symmetry test")
212  tap_check_small(sigma_2,3*epsilon(sigma_2),"xz-zx symmetry test")
213  tap_check_small(sigma_3,3*epsilon(sigma_3),"zy-yz symmetry test")
214  tap_check_small(sigma_4,3*epsilon(sigma_4),"zx-zy symmetry test")
215  tap_check_small(sigma_5,3*epsilon(sigma_5),"yx-yz symmetry test")
216  tap_check_small(sigma_6,3*epsilon(sigma_6),"xy-xz symmetry test")
217 
218 
219  tap_done
220 
221 CONTAINS
222 
223  SUBROUTINE makeconfig(Sim, config)
224  IMPLICIT NONE
225  !------------------------------------------------------------------------!
226  CLASS(fosite) :: Sim
227  TYPE(Dict_TYP),POINTER :: config
228  !------------------------------------------------------------------------!
229  ! Local variable declaration
230  TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
231  sources, timedisc, fluxes
232  REAL :: dynvis
233  !------------------------------------------------------------------------!
234  ! mesh settings
235  mesh => dict( &
236  "meshtype" / midpoint, &
237  "geometry" / mgeo, &
238  "inum" / res, & ! resolution in x, !
239  "jnum" / res, & ! y and !
240  "knum" / res, & ! z direction !
241  "xmin" / (-0.5*xyzlen), &
242  "xmax" / (0.5*xyzlen), &
243  "ymin" / (-0.5*xyzlen), &
244  "ymax" / (0.5*xyzlen), &
245  "zmin" / (-0.5*xyzlen), &
246  "zmax" / (0.5*xyzlen), &
247  "output/dl" / 0, &
248  "output/bh" / 0, &
249  "output/rotation" / 0 &
250  )
251 
252  ! physics settings
253  physics => dict( &
254  "problem" / euler, &
255  "gamma" / gamma & ! ratio of specific heats !
256  )
257 
258  ! flux calculation and reconstruction method
259  fluxes => dict( &
260  "fluxtype" / kt, &
261  "order" / linear, &
262  "variables" / conservative, & ! vars. to use for reconstruction !
263  "limiter" / monocent, & ! one of: minmod, monocent,... !
264  "theta" / 1.2 & ! optional parameter for limiter !
265  )
266 
267  ! boundary conditions
268  boundary => dict( &
269  "western" / periodic, &
270  "eastern" / periodic, &
271  "southern" / periodic, &
272  "northern" / periodic, &
273  "bottomer" / periodic, &
274  "topper" / periodic &
275  )
276 
277 
278  NULLIFY(sources)
279  ! viscosity source term
280  ! compute dynamic viscosity constant using typical scales and Reynolds number
281  dynvis = abs(rho0 * xyzlen * (v0-v1) / re)
282  IF (dynvis.GT.tiny(1.0)) THEN
283  sources => dict( &
284  "vis/stype" / viscosity, &
285  "vis/vismodel" / molecular, &
286  "vis/dynconst" / dynvis, &
287  "vis/bulkconst" / (-2./3.*dynvis), &
288  "vis/output/dynvis" / 0, &
289  "vis/output/stress" / 0, &
290  "vis/output/kinvis" / 0, &
291  "vis/output/bulkvis" / 0 &
292  )
293  END IF
294 
295  ! time discretization settings
296  timedisc => dict( &
297  "method" / modified_euler, &
298  "order" / 3, &
299  "cfl" / 0.4, &
300  "stoptime" / tsim, &
301  "dtlimit" / 1.0e-10, &
302  "maxiter" / 100000 &
303 ! "output/energy" / 1, &
304 ! "output/xmomentum" / 1, &
305 ! "output/ymomentum" / 1, &
306 ! "output/zmomentum" / 1, &
307 ! "output/rhs" / 1 &
308  )
309 
310  ! initialize data input/output
311  datafile => dict(&
312 ! "fileformat" / HDF, &
313  "fileformat" / vtk, &
314 ! "fileformat" / XDMF, &
315 ! "fileformat" / GNUPLOT, "filecycles" / 0, &
316 ! "fileformat" / BINARY, &
317 ! "fileformat" / NETCDF, &
318  "filename" / (trim(odir) // trim(ofname)), &
319  "count" / onum &
320  )
321 
322  config => dict( &
323  "mesh" / mesh, &
324  "physics" / physics, &
325  "boundary" / boundary, &
326  "fluxes" / fluxes, &
327  "timedisc" / timedisc, &
328 ! "logfile" / logfile, &
329  "datafile" / datafile &
330  )
331 
332  IF (ASSOCIATED(sources)) &
333  CALL setattr(config, "sources", sources)
334  END SUBROUTINE makeconfig
335 
336 
337 
338 
339  SUBROUTINE initdata(Mesh,Physics,Timedisc,flow_dir,reuse_random_seed)
340  IMPLICIT NONE
341  !------------------------------------------------------------------------!
342  CLASS(physics_base) :: Physics
343  CLASS(mesh_base) :: Mesh
344  CLASS(timedisc_base) :: Timedisc
345  LOGICAL, OPTIONAL :: reuse_random_seed
346  CHARACTER(2), OPTIONAL :: flow_dir
347  !------------------------------------------------------------------------!
348  ! Local variable declaration
349  INTEGER :: i,j,k
350  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3) :: dv,dv_temp
351  INTEGER :: clock
352  !------------------------------------------------------------------------!
353  INTENT(IN) :: mesh,physics,flow_dir,reuse_random_seed
354  INTENT(INOUT) :: timedisc
355  !------------------------------------------------------------------------!
356  ! Seed the random number generator
357  IF (.NOT.(PRESENT(reuse_random_seed).AND.reuse_random_seed)) THEN
358  ! initialize with a mix from current time and mpi rank
359  CALL system_clock(count=clock)
360  seed = clock + (timedisc%getrank()+1) * (/(i-1, i=1,n)/)
361  END IF
362  CALL random_seed(put=seed)
363  CALL random_number(dv)
364  ! initial condition
365  IF (PRESENT(flow_dir).AND.flow_dir.EQ.'yx') THEN
366  ! flow along yx-direction:
367  print *,"yx-Direction"
368  ! x and z-velocity vanish everywhere
369  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.
370  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = 0.
371 
372  WHERE ((mesh%bcenter(:,:,:,1).LT.(mesh%xmin+0.25*xyzlen)).OR. &
373  (mesh%bcenter(:,:,:,1).GT.(mesh%xmin+0.75*xyzlen)))
374  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho0
375  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = v0
376  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p0
377  ELSEWHERE
378  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho1
379  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = v1
380  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p1
381  END WHERE
382  DO k = mesh%KGMIN, mesh%KGMAX
383  ! add perturbation to the velocity field
384  timedisc%pvar%data4d(:,:,k,physics%XVELOCITY) = timedisc%pvar%data4d(:,:,k,physics%XVELOCITY) &
385  + transpose(dv(:,:,k,2)-0.5)*0.2
386  timedisc%pvar%data4d(:,:,k,physics%YVELOCITY) = timedisc%pvar%data4d(:,:,k,physics%YVELOCITY) &
387  + transpose(dv(:,:,k,1)-0.5)*0.2
388  timedisc%pvar%data4d(:,:,k,physics%ZVELOCITY) = timedisc%pvar%data4d(:,:,k,physics%ZVELOCITY) &
389  + transpose(dv(:,:,k,3)-0.5)*0.2
390  END DO
391  ! initial condition
392  ELSEIF (PRESENT(flow_dir).AND.flow_dir.EQ.'zx') THEN
393  ! flow along zx-direction:
394  print *,"zx-Direction"
395  ! x and y-velocity vanish everywhere
396  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.
397  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.
398 
399  WHERE ((mesh%bcenter(:,:,:,1).LT.(mesh%xmin+0.25*xyzlen)).OR. &
400  (mesh%bcenter(:,:,:,1).GT.(mesh%xmin+0.75*xyzlen)))
401  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho0
402  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = v0
403  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p0
404  ELSEWHERE
405  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho1
406  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = v1
407  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p1
408  END WHERE
409  DO k = mesh%KGMIN,mesh%KGMAX
410  dv_temp(:,:,k,1) = transpose(dv(:,:,k,2))
411  dv_temp(:,:,k,2) = transpose(dv(:,:,k,1))
412  dv_temp(:,:,k,3) = transpose(dv(:,:,k,3))
413  END DO
414  DO i= mesh%IGMIN, mesh%IGMAX
415  ! add perturbation to the velocity field
416  timedisc%pvar%data4d(i,:,:,physics%XVELOCITY) = timedisc%pvar%data4d(i,:,:,physics%XVELOCITY) &
417  + transpose(dv_temp(i,:,:,1)-0.5)*0.2
418  timedisc%pvar%data4d(i,:,:,physics%YVELOCITY) = timedisc%pvar%data4d(i,:,:,physics%YVELOCITY) &
419  + transpose(dv_temp(i,:,:,3)-0.5)*0.2
420  timedisc%pvar%data4d(i,:,:,physics%ZVELOCITY) = timedisc%pvar%data4d(i,:,:,physics%ZVELOCITY) &
421  + transpose(dv_temp(i,:,:,2)-0.5)*0.2
422  END DO
423  ELSEIF (PRESENT(flow_dir).AND.flow_dir.EQ.'xy') THEN
424  ! flow along xy-direction:
425  print *,"xy-Direction"
426  ! y and z-velocity vanish everywhere
427  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.
428  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = 0.
429  WHERE ((mesh%bcenter(:,:,:,2).LT.(mesh%ymin+0.25*xyzlen)).OR. &
430  (mesh%bcenter(:,:,:,2).GT.(mesh%ymin+0.75*xyzlen)))
431  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho0
432  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = v0
433  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p0
434  ELSEWHERE
435  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho1
436  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = v1
437  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p1
438  END WHERE
439  ! add perturbation to the velocity field
440  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) &
441  + (dv(:,:,:,1)-0.5)*0.2
442  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) &
443  + (dv(:,:,:,2)-0.5)*0.2
444  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) &
445  + (dv(:,:,:,3)-0.5)*0.2
446 
447  ELSEIF (PRESENT(flow_dir).AND.flow_dir.EQ.'xz') THEN
448  ! flow along xz-direction:
449  print *,"xz-Direction"
450  ! y and z-velocity vanish everywhere
451  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.
452  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = 0.
453  WHERE ((mesh%bcenter(:,:,:,3).LT.(mesh%zmin+0.25*xyzlen)).OR. &
454  (mesh%bcenter(:,:,:,3).GT.(mesh%zmin+0.75*xyzlen)))
455  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho0
456  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = v0
457  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p0
458  ELSEWHERE
459  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho1
460  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = v1
461  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p1
462  END WHERE
463  DO i = mesh%IGMIN, mesh%IGMAX
464  ! add perturbation to the velocity field
465  timedisc%pvar%data4d(i,:,:,physics%XVELOCITY) = timedisc%pvar%data4d(i,:,:,physics%XVELOCITY) &
466  + transpose(dv(i,:,:,1)-0.5)*0.2
467  timedisc%pvar%data4d(i,:,:,physics%YVELOCITY) = timedisc%pvar%data4d(i,:,:,physics%YVELOCITY) &
468  + transpose(dv(i,:,:,3)-0.5)*0.2
469  timedisc%pvar%data4d(i,:,:,physics%ZVELOCITY) = timedisc%pvar%data4d(i,:,:,physics%ZVELOCITY) &
470  + transpose(dv(i,:,:,2)-0.5)*0.2
471  END DO
472  ELSEIF (PRESENT(flow_dir).AND.flow_dir.EQ.'zy') THEN
473  ! flow along zy-direction:
474  print *,"zy-Direction"
475  ! x and y-velocity vanish everywhere
476  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.
477  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.
478  WHERE ((mesh%bcenter(:,:,:,2).LT.(mesh%ymin+0.25*xyzlen)).OR. &
479  (mesh%bcenter(:,:,:,2).GT.(mesh%ymin+0.75*xyzlen)))
480  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho0
481  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = v0
482  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p0
483  ELSEWHERE
484  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho1
485  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = v1
486  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p1
487  END WHERE
488  DO j = mesh%JGMIN, mesh%JGMAX
489  ! add perturbation to the velocity field
490  timedisc%pvar%data4d(:,j,:,physics%XVELOCITY) = timedisc%pvar%data4d(:,j,:,physics%XVELOCITY) &
491  + transpose(dv(:,j,:,3)-0.5)*0.2
492  timedisc%pvar%data4d(:,j,:,physics%YVELOCITY) = timedisc%pvar%data4d(:,j,:,physics%YVELOCITY) &
493  + transpose(dv(:,j,:,2)-0.5)*0.2
494  timedisc%pvar%data4d(:,j,:,physics%ZVELOCITY) = timedisc%pvar%data4d(:,j,:,physics%ZVELOCITY) &
495  + transpose(dv(:,j,:,1)-0.5)*0.2
496  END DO
497  ELSEIF (PRESENT(flow_dir).AND.flow_dir.EQ.'yz') THEN
498  ! flow along yz-direction:
499  print *,"yz-Direction"
500  ! x and z-velocity vanish everywhere
501  timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.
502  timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = 0.
503  WHERE ((mesh%bcenter(:,:,:,3).LT.(mesh%zmin+0.25*xyzlen)).OR. &
504  (mesh%bcenter(:,:,:,3).GT.(mesh%zmin+0.75*xyzlen)))
505  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho0
506  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = v0
507  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p0
508  ELSEWHERE
509  timedisc%pvar%data4d(:,:,:,physics%DENSITY) = rho1
510  timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = v1
511  timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = p1
512  END WHERE
513  DO i = mesh%IGMIN, mesh%IGMAX
514  dv_temp(i,:,:,1) = transpose(dv(i,:,:,1))
515  dv_temp(i,:,:,2) = transpose(dv(i,:,:,3))
516  dv_temp(i,:,:,3) = transpose(dv(i,:,:,2))
517  END DO
518  DO k = mesh%KGMIN, mesh%KGMAX
519  ! add perturbation to the velocity field
520  timedisc%pvar%data4d(:,:,k,physics%XVELOCITY) = timedisc%pvar%data4d(:,:,k,physics%XVELOCITY) &
521  + transpose(dv_temp(:,:,k,2)-0.5)*0.2
522  timedisc%pvar%data4d(:,:,k,physics%YVELOCITY) = timedisc%pvar%data4d(:,:,k,physics%YVELOCITY) &
523  + transpose(dv_temp(:,:,k,1)-0.5)*0.2
524  timedisc%pvar%data4d(:,:,k,physics%ZVELOCITY) = timedisc%pvar%data4d(:,:,k,physics%ZVELOCITY) &
525  + transpose(dv_temp(:,:,k,3)-0.5)*0.2
526  END DO
527  ELSE
528  CALL mesh%Error("KHI INIT","Unknown flow direction")
529  END IF
530 
531  CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
532  CALL mesh%Info(" DATA-----> initial condition: " // &
533  "Kelvin-Helmholtz instability")
534 
535  END SUBROUTINE initdata
536 
537 END PROGRAM khi
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
makeconfig
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
fosite_mod
Definition: fosite.f90:41
khi
program khi
Kelvin-Helmholtz instability.
Definition: KHI2D.f90:37
Generated by   doxygen 1.8.14