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-2018 #
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 !----------------------------------------------------------------------------!
37 PROGRAM 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  ! 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 = 50 ! 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 = 'KHI2D'
63  !--------------------------------------------------------------------------!
64  CLASS(fosite), ALLOCATABLE :: sim
65  CLASS(marray_compound), POINTER :: pvar,pvar_init
66  INTEGER :: n
67  REAL :: sigma
68  INTEGER, DIMENSION(:), ALLOCATABLE :: seed
69  !--------------------------------------------------------------------------!
70  tap_plan(1)
71 
72  ! allocate memory for dynamic types and arrays
73  CALL random_seed(size=n)
74  ALLOCATE(sim,pvar,seed(n))
75 
76  ! simulate KHI with initial x-velocity along x-direction
77  CALL sim%InitFosite()
78  CALL makeconfig(sim,sim%config)
79  CALL sim%Setup()
80  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
81  ! store transposed initial data
82  CALL sim%Physics%new_statevector(pvar_init,primitive)
83  CALL transposedata(sim%Mesh,sim%Timedisc%pvar,pvar_init,"xy")
84  CALL sim%Run()
85  ! store transposed result of the first run
86  CALL sim%Physics%new_statevector(pvar,primitive)
87  CALL transposedata(sim%Mesh,sim%Timedisc%pvar,pvar,"xy")
88  ! finish the simulation
89  CALL sim%Finalize()
90  DEALLOCATE(sim)
91 
92  ! simulate KHI with initial y-velocity along y-direction
93  ALLOCATE(sim)
94  CALL sim%InitFosite()
95  CALL makeconfig(sim, sim%config)
96  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_rotate"))
97  CALL sim%Setup()
98  sim%Timedisc%pvar%data1d(:) = pvar_init%data1d(:)
99  CALL sim%Physics%Convert2Conservative(sim%Timedisc%pvar,sim%Timedisc%cvar)
100  CALL sim%Run()
101  ! compare results
102  sigma = sqrt(sum((sim%Timedisc%pvar%data4d(:,:,:,:)-pvar%data4d(:,:,:,:))**2)/SIZE(pvar%data4d))
103 
104  CALL pvar%Destroy()
105  CALL sim%Finalize()
106  DEALLOCATE(pvar,seed,sim)
107  tap_check_small(sigma,tiny(sigma),"x-y symmetry test")
108  tap_done
109 
110 CONTAINS
111 
112  SUBROUTINE makeconfig(Sim, config)
113  IMPLICIT NONE
114  !------------------------------------------------------------------------!
115  CLASS(fosite),INTENT(INOUT) :: Sim
116  TYPE(Dict_TYP),POINTER :: config
117  !------------------------------------------------------------------------!
118  ! Local variable declaration
119  TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
120  sources, timedisc, fluxes
121  REAL :: dynvis
122  !------------------------------------------------------------------------!
123  ! mesh settings
124  mesh => dict( &
125  "meshtype" / midpoint, &
126  "geometry" / mgeo, &
127  "inum" / res, & ! resolution in x, !
128  "jnum" / res, & ! y and !
129  "knum" / 1, & ! z direction !
130  "xmin" / (-0.5*xyzlen), &
131  "xmax" / (0.5*xyzlen), &
132  "ymin" / (-0.5*xyzlen), &
133  "ymax" / (0.5*xyzlen), &
134 ! "zmin" / (-0.5*XYZLEN), &
135 ! "zmax" / (0.5*XYZLEN) &
136  "zmin" / (-0.0), &
137  "zmax" / (0.0) &
138  )
139 
140  ! physics settings
141  physics => dict( &
142  "problem" / euler, &
143  "gamma" / gamma & ! ratio of specific heats !
144  )
145 
146  ! flux calculation and reconstruction method
147  fluxes => dict( &
148  "fluxtype" / kt, &
149  "order" / linear, &
150  "variables" / conservative, & ! vars. to use for reconstruction !
151 ! "variables" / PRIMITIVE, & ! vars. to use for reconstruction !
152  "limiter" / monocent, & ! one of: minmod, monocent,... !
153 ! "output/pfluxes" / 1, &
154 ! "output/pstates" / 1, &
155 ! "output/cstates" / 1, &
156  "theta" / 1.2 & ! optional parameter for limiter !
157  )
158 
159  ! boundary conditions
160  boundary => dict( &
161  "western" / periodic, &
162  "eastern" / periodic, &
163  "southern" / periodic, &
164  "northern" / periodic, &
165  "bottomer" / periodic, &
166  "topper" / periodic &
167  )
168 
169  NULLIFY(sources)
170  ! viscosity source term
171  ! compute dynamic viscosity constant using typical scales and Reynolds number
172  dynvis = abs(rho0 * xyzlen * (v0-v1) / re)
173  IF (dynvis.GT.tiny(1.0)) THEN
174  sources => dict( &
175  "vis/stype" / viscosity, &
176  "vis/vismodel" / molecular, &
177  "vis/dynconst" / dynvis, &
178  "vis/bulkconst" / (-2./3.*dynvis), &
179  "vis/output/dynvis" / 0, &
180  "vis/output/stress" / 0, &
181  "vis/output/kinvis" / 0, &
182  "vis/output/bulkvis" / 0 &
183  )
184  END IF
185 
186  ! time discretization settings
187  timedisc => dict( &
188  "method" / modified_euler, &
189  "order" / 3, &
190  "cfl" / 0.4, &
191  "stoptime" / tsim, &
192  "dtlimit" / 1.0e-10, &
193  "maxiter" / 100000, &
194  "output/pressure" / 1, &
195  "output/density" / 1, &
196  "output/xvelocity" / 1, &
197  "output/yvelocity" / 1 &
198  )
199 
200  ! initialize data input/output
201  datafile => dict(&
202  "fileformat" / vtk, &
203 ! "fileformat" / XDMF, &
204  "filename" / (trim(odir) // trim(ofname)), &
205  "count" / onum &
206  )
207 
208  config => dict( &
209  "mesh" / mesh, &
210  "physics" / physics, &
211  "boundary" / boundary, &
212  "fluxes" / fluxes, &
213  "timedisc" / timedisc, &
214  "datafile" / datafile &
215  )
216  IF (ASSOCIATED(sources)) &
217  CALL setattr(config, "sources", sources)
218  END SUBROUTINE makeconfig
219 
220  SUBROUTINE initdata(Mesh,Physics,Timedisc)
221  IMPLICIT NONE
222  !------------------------------------------------------------------------!
223  CLASS(physics_base) :: Physics
224  CLASS(mesh_base) :: Mesh
225  CLASS(timedisc_base) :: Timedisc
226  !------------------------------------------------------------------------!
227  ! Local variable declaration
228  INTEGER :: k,clock
229  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,2) :: dv
230  !------------------------------------------------------------------------!
231  INTENT(IN) :: mesh,physics
232  INTENT(INOUT) :: timedisc
233  !------------------------------------------------------------------------!
234  ! Seed the random number generator
235  ! initialize with a mix from current time and mpi rank
236  CALL system_clock(count=clock)
237  seed = clock + (timedisc%getrank()+1) * (/(k-1, k=1,n)/)
238  CALL random_seed(put=seed)
239  CALL random_number(dv)
240 
241  ! initial condition
242  SELECT TYPE(pvar => timedisc%pvar)
243  TYPE IS(statevector_euler) ! non-isothermal HD
244  ! flow along x-direction:
245  WHERE ((mesh%bcenter(:,:,:,2).LT.(mesh%ymin+0.25*xyzlen)).OR. &
246  (mesh%bcenter(:,:,:,2).GT.(sim%Mesh%ymin+0.75*xyzlen)))
247  pvar%density%data3d(:,:,:) = rho0
248  pvar%velocity%data4d(:,:,:,1) = v0
249  pvar%velocity%data4d(:,:,:,2) = 0.0
250  pvar%pressure%data3d(:,:,:) = p0
251  ELSEWHERE
252  pvar%density%data3d(:,:,:) = rho1
253  pvar%velocity%data4d(:,:,:,1) = v1
254  pvar%velocity%data4d(:,:,:,2) = 0.0
255  pvar%pressure%data3d(:,:,:) = p1
256  END WHERE
257  ! add perturbations
258  DO k=mesh%KGMIN,mesh%KGMAX
259  pvar%velocity%data4d(:,:,k,1) = pvar%velocity%data4d(:,:,k,1) &
260  + (dv(:,:,k,1)-0.5)*0.02
261  pvar%velocity%data4d(:,:,k,2) = pvar%velocity%data4d(:,:,k,2) &
262  + (dv(:,:,k,2)-0.5)*0.02
263  END DO
264  CLASS DEFAULT
265  CALL physics%Error("KHI2D::InitData","only non-isothermal HD supported")
266  END SELECT
267 
268  ! initial condition
269  CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
270  CALL mesh%Info(" DATA-----> initial condition: " // &
271  "Kelvin-Helmholtz instability")
272  END SUBROUTINE initdata
273 
274  SUBROUTINE transposedata(Mesh,pvar_in,pvar_out,dir)
275  IMPLICIT NONE
276  !------------------------------------------------------------------------!
277  CLASS(mesh_base), INTENT(IN) :: Mesh
278  CLASS(marray_compound), INTENT(INOUT) :: pvar_in,pvar_out
279  CHARACTER(LEN=2), INTENT(IN) :: dir
280  !------------------------------------------------------------------------!
281  INTEGER :: k
282  !------------------------------------------------------------------------!
283  SELECT TYPE(pin => pvar_in)
284  TYPE IS(statevector_euler) ! non-isothermal HD
285  SELECT TYPE(pout => pvar_out)
286  TYPE IS(statevector_euler) ! non-isothermal HD
287  SELECT CASE(dir)
288  CASE("xy")
289  ! transpose x-y directions
290  DO k=mesh%KGMIN,mesh%KGMAX
291  ! simply transpose the 1st and 2nd indices ...
292  pout%density%data3d(:,:,k) = transpose(pin%density%data3d(:,:,k))
293  pout%pressure%data3d(:,:,k) = transpose(pin%pressure%data3d(:,:,k))
294  ! ... and exchange x- and y-velocities
295  pout%velocity%data4d(:,:,k,1) = transpose(pin%velocity%data4d(:,:,k,2))
296  pout%velocity%data4d(:,:,k,2) = transpose(pin%velocity%data4d(:,:,k,1))
297  END DO
298  CASE DEFAULT
299  CALL mesh%Error("KHI2D::TransposeData","directions must be one of 'xy','xz' or 'yz'")
300  END SELECT
301  END SELECT
302  END SELECT
303  END SUBROUTINE transposedata
304 
305 END PROGRAM khi
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
transposedata
subroutine transposedata(Mesh, pvar_in, pvar_out, dir)
Definition: KHI2D.f90:275
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