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