43 REAL,
PARAMETER :: tsim = 30.0
44 REAL,
PARAMETER :: gamma = 1.4
45 REAL,
PARAMETER :: re = 1.0e+4
47 REAL,
PARAMETER :: rho0 = 1.0
48 REAL,
PARAMETER :: v0 =-0.5
49 REAL,
PARAMETER :: p0 = 2.5
50 REAL,
PARAMETER :: rho1 = 2.0
51 REAL,
PARAMETER :: v1 = 0.5
52 REAL,
PARAMETER :: p1 = p0
54 INTEGER,
PARAMETER :: mgeo = cartesian
55 INTEGER,
PARAMETER :: res = 8
56 REAL,
PARAMETER :: xyzlen= 1.0
58 INTEGER,
PARAMETER :: onum = 1
59 CHARACTER(LEN=256),
PARAMETER &
61 CHARACTER(LEN=256),
PARAMETER &
64 CLASS(
fosite),
ALLOCATABLE :: sim
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
75 CALL sim%Warning(
"KHI3d::main",
"in parallel mode the symmetry test currently works only for xy-direction")
81 CALL random_seed(size=n)
84 CALL system_clock(count=clock)
85 seed = clock + (sim%GetRank()+1) * (/(i-1, i=1,n)/)
86 CALL random_seed(put=seed)
93 CALL setattr(sim%config,
"/mesh/decomposition", (/1,1,-1/))
95 CALL setattr(sim%config,
"/datafile/filename", (trim(odir) // trim(ofname) //
"_xy"))
98 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,
'xy')
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.)
108 CALL sim%InitFosite()
112 CALL setattr(sim%config,
"/mesh/decomposition", (/1,1,-1/))
115 CALL setattr(sim%config,
"/datafile/filename", (trim(odir) // trim(ofname) //
"_yx"))
117 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,
'yx')
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")
123 CALL sim%Finalize(mpifinalize_=.false.)
129 CALL sim%InitFosite()
131 CALL setattr(sim%config,
"/datafile/filename", (trim(odir) // trim(ofname) //
"_xz"))
133 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,
'xz')
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.)
143 CALL sim%InitFosite()
145 CALL setattr(sim%config,
"/datafile/filename", (trim(odir) // trim(ofname) //
"_zx"))
147 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,
'zx')
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.)
158 CALL sim%InitFosite()
160 CALL setattr(sim%config,
"/datafile/filename", (trim(odir) // trim(ofname) //
"_zy"))
162 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,
'zy')
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.)
172 CALL sim%InitFosite()
174 CALL setattr(sim%config,
"/datafile/filename", (trim(odir) // trim(ofname) //
"_yz"))
176 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc,
'yz')
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")
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))
192 sigma_1 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_yx(:,:,:,:))**2)/
SIZE(pvar_temp))
193 DEALLOCATE(pvar_temp)
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))
205 sigma_2 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_zx(:,:,:,:))**2)/
SIZE(pvar_temp))
206 DEALLOCATE(pvar_temp)
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))
217 sigma_3 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_yz(:,:,:,:))**2)/
SIZE(pvar_temp))
218 DEALLOCATE(pvar_temp)
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))
229 sigma_4 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_zy(:,:,:,:))**2)/
SIZE(pvar_temp))
230 DEALLOCATE(pvar_temp)
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))
241 sigma_5 = sqrt(sum((pvar_temp(:,:,:,:) - pvar_yz(:,:,:,:))**2)/
SIZE(pvar_temp))
242 DEALLOCATE(pvar_temp)
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))
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)
258 DEALLOCATE(sim,pvar_xy,pvar_yx,dv,dv_temp,seed)
260 tap_check_small(sigma_1,4*epsilon(sigma_1),
"xy-yx symmetry test")
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")
276 TYPE(dict_typ),
POINTER :: config
279 TYPE(dict_typ),
POINTER :: mesh, physics, boundary, datafile, &
280 sources, timedisc, fluxes
285 "meshtype" / midpoint, &
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), &
298 "output/rotation" / 0 &
311 "variables" / conservative, &
312 "limiter" / monocent, &
318 "western" / periodic, &
319 "eastern" / periodic, &
320 "southern" / periodic, &
321 "northern" / periodic, &
322 "bottomer" / periodic, &
323 "topper" / periodic &
330 dynvis = abs(rho0 * xyzlen * (v0-v1) / re)
331 IF (dynvis.GT.tiny(1.0))
THEN
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 &
346 "method" / modified_euler, &
350 "dtlimit" / 1.0e-10, &
361 "fileformat" / vtk, &
364 "filename" / (trim(odir) // trim(ofname)), &
370 "physics" / physics, &
371 "boundary" / boundary, &
373 "timedisc" / timedisc, &
375 "datafile" / datafile &
378 IF (
ASSOCIATED(sources)) &
379 CALL setattr(config,
"sources", sources)
383 SUBROUTINE initdata(Mesh,Physics,Timedisc,flow_dir)
386 CLASS(physics_base) :: Physics
387 CLASS(mesh_base) :: Mesh
388 CLASS(timedisc_base) :: Timedisc
389 CHARACTER(2) :: flow_dir
394 INTENT(IN) :: mesh,physics,flow_dir
395 INTENT(INOUT) :: timedisc
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)
402 CALL random_number(dv)
404 IF (.NOT.
ALLOCATED(dv_temp))
ALLOCATE(dv_temp,mold=pvar%velocity%data4d,stat=sim%err)
407 SELECT CASE(flow_dir)
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
418 pvar%density%data3d(:,:,:) = rho1
419 pvar%velocity%data4d(:,:,:,1) = v1
420 pvar%pressure%data3d(:,:,:) = p1
423 pvar%velocity%data4d(:,:,:,:) = pvar%velocity%data4d(:,:,:,:) + (dv(:,:,:,:)-0.5)*0.2
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
434 pvar%density%data3d(:,:,:) = rho1
435 pvar%velocity%data4d(:,:,:,2) = v1
436 pvar%pressure%data3d(:,:,:) = p1
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
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
454 pvar%density%data3d(:,:,:) = rho1
455 pvar%velocity%data4d(:,:,:,1) = v1
456 pvar%pressure%data3d(:,:,:) = p1
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
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
474 pvar%density%data3d(:,:,:) = rho1
475 pvar%velocity%data4d(:,:,:,3) = v1
476 pvar%pressure%data3d(:,:,:) = p1
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))
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
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
500 pvar%density%data3d(:,:,:) = rho1
501 pvar%velocity%data4d(:,:,:,3) = v1
502 pvar%pressure%data3d(:,:,:) = p1
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
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
520 pvar%density%data3d(:,:,:) = rho1
521 pvar%velocity%data4d(:,:,:,2) = v1
522 pvar%pressure%data3d(:,:,:) = p1
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))
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
537 CALL mesh%Error(
"KHI INIT",
"Unknown flow direction")
541 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
542 CALL mesh%Info(
" DATA-----> initial condition: " // &
543 "Kelvin-Helmholtz instability " // flow_dir //
"-direction")
program khi
Kelvin-Helmholtz instability.
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
subroutine makeconfig(Sim, config)