poiseuille.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: poiseuille.f90 #
5!# #
6!# Copyright (C) 2006-2023 #
7!# Bjoern Sperling <sperling@astrophysik.uni-kiel.de> #
8!# Tobias Illenseer <tillense@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!#############################################################################
32!----------------------------------------------------------------------------!
33PROGRAM poiseuille
34 USE fosite_mod
35#include "tap.h"
36 IMPLICIT NONE
37 !--------------------------------------------------------------------------!
38 ! simulation parameters
39 REAL, PARAMETER :: tsim = 2.0e+0 ! simulation time [TVIS] !
40 REAL, PARAMETER :: re = 1.0e+0 ! Reynolds number !
41 REAL, PARAMETER :: ma = 1.0e-2 ! Mach number << 1 !
42 REAL, PARAMETER :: rho0 = 1.0e-0 ! initial density !
43 REAL, PARAMETER :: eta = 1.0e-0 ! dynamic viscosity (const.) !
44 REAL, PARAMETER :: gam = 1.4 ! ratio of spec. heat coeff. !
45 ! mesh settings
46 INTEGER, PARAMETER :: mgeo = cylindrical
47! INTEGER, PARAMETER :: MGEO = CARTESIAN
48 INTEGER, PARAMETER :: xres = 30 ! radial resolution !
49 INTEGER, PARAMETER :: yres = 1 ! azimuthal resolution !
50 INTEGER, PARAMETER :: zres = 5 ! resolution along the tube !
51 REAL, PARAMETER :: ltube = 10.0 ! length of the tube !
52 REAL, PARAMETER :: rtube = 1.0 ! radius of the tube !
53 !--------------------------------------------------------------------------!
54 ! output file parameter
55 INTEGER, PARAMETER :: onum = 10 ! number of output data sets
56 CHARACTER(LEN=256), PARAMETER & ! output data dir
57 :: odir = './'
58 CHARACTER(LEN=256), PARAMETER & ! output data file name
59 :: ofname = 'poiseuille'
60 ! derived parameters
61 REAL :: tvis ! viscous timescale !
62 REAL :: umax ! max. velocity of lam. flow !
63 REAL :: pin,pout ! inlet/outlet pressure !
64 !--------------------------------------------------------------------------!
65 CLASS(fosite),ALLOCATABLE :: sim
66 LOGICAL :: ok
67 !--------------------------------------------------------------------------!
68
69tap_plan(1)
70
71 ALLOCATE(sim)
72 CALL sim%InitFosite()
73
74 CALL makeconfig(sim, sim%config)
75
76! CALL PrintDict(config)
77
78 CALL sim%Setup()
79
80 ! set initial condition
81 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
82
83 CALL sim%Run()
84 ok = .NOT.sim%aborted
85
86 CALL sim%Finalize()
87 DEALLOCATE(sim)
88
89tap_check(ok,"stoptime reached")
90
91tap_done
92
93CONTAINS
94
95 SUBROUTINE makeconfig(Sim, config)
96 IMPLICIT NONE
97 !------------------------------------------------------------------------!
98 CLASS(fosite) :: Sim
99 TYPE(dict_typ),POINTER :: config
100 !------------------------------------------------------------------------!
101 ! Local variable declaration
102 INTEGER :: bc(6)
103 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, sources, &
104 vis, timedisc, fluxes
105 REAL :: x1,x2,y1,y2,z1,z2
106 !------------------------------------------------------------------------!
107 INTENT(INOUT) :: sim
108 !------------------------------------------------------------------------!
109 ! compute derived parameters
110 ! set velocity scale according to the given Reynolds number
111 umax = eta/rho0 * re/rtube
112 ! viscous time scale
113 tvis = rho0/eta * rtube**2
114 ! inlet pressure determined by Mach number
115 pin = rho0/gam * (umax/ma)**2
116 ! outlet pressure (for given length and maximum velocity of laminar solution
117 pout = pin * (1.0 - 4*gam * (ma**2/re) * (ltube/rtube) )
118 IF (pout .LT. 0.0) &
119 CALL sim%Error("poiseuille::MakeConfig","negative outlet pressure")
120
121 ! geometry dependent setttings
122 SELECT CASE(mgeo)
123 CASE(cylindrical)
124 x1 = 0.0
125 x2 = rtube
126 y1 = 0.0
127 y2 = 2*pi
128 z1 = 0.0
129 z2 = ltube
130 bc(west) = axis
131! bc(EAST) = NOSLIP
132 bc(east) = custom
133 bc(south) = periodic
134 bc(north) = periodic
135 bc(bottom)= fixed
136 bc(top) = fixed
137 CASE(cartesian)
138 x1 = 0.0
139 x2 = rtube
140 y1 = 0.0
141 y2 = rtube
142 z1 = 0.0
143 z2 = ltube
144 bc(west) = noslip
145 bc(east) = noslip
146 bc(south) = noslip
147 bc(north) = noslip
148 bc(bottom)= fixed
149 bc(top) = fixed
150 END SELECT
151
152 ! mesh settings
153 mesh => dict("meshtype" / midpoint, &
154 "geometry" / mgeo, &
155! "decomposition" / (/1,-1,1/), &
156 "inum" / xres, &
157 "jnum" / yres, &
158 "knum" / zres, &
159 "xmin" / x1, &
160 "xmax" / x2, &
161 "ymin" / y1, &
162 "ymax" / y2, &
163 "zmin" / z1, &
164 "zmax" / z2)
165
166 ! boundary conditions
167 boundary => dict("western" / bc(west), &
168 "eastern" / bc(east), &
169 "southern" / bc(south), &
170 "northern" / bc(north), &
171 "bottomer" / bc(bottom), &
172 "topper" / bc(top))
173
174 ! physics settings
175 physics => dict("problem" / euler, &
176 "gamma" / gam) ! ratio of specific heats !
177
178 ! flux calculation and reconstruction method
179 fluxes => dict("order" / linear, &
180 "fluxtype" / kt, &
181 "variables" / primitive, & ! vars. to use for reconstruction!
182 "limiter" / vanleer, & ! one of: minmod, monocent,... !
183 "theta" / 1.2) ! optional parameter for limiter !
184
185 ! viscosity source term
186 vis => dict("stype" / viscosity, &
187 "vismodel" / molecular, &
188 "dynconst" / eta, &
189 "bulkconst" / (-2./3.*eta), &
190 "cvis" / 0.4) ! viscous Courant number
191
192 sources => dict("vis" / vis)
193
194 ! time discretization settings
195 timedisc => dict("method" / ssprk, &
196 "order" / 3, &
197 "cfl" / 0.4, &
198 "stoptime" / (tsim*tvis), &
199 "dtlimit" / (1.0e-7*tvis), &
200 "maxiter" / 100000000)
201
202 ! initialize data input/output
203 datafile => dict( &
204 "fileformat" / vtk, &
205 "filename" / (trim(odir) // trim(ofname)), &
206 "count" / onum)
207
208 config => dict("mesh" / mesh, &
209 "physics" / physics, &
210 "boundary" / boundary, &
211 "fluxes" / fluxes, &
212 "sources" / sources, &
213 "timedisc" / timedisc, &
214 "datafile" / datafile)
215 END SUBROUTINE makeconfig
216
217 SUBROUTINE initdata(Mesh,Physics,Timedisc)
218 IMPLICIT NONE
219 !------------------------------------------------------------------------!
220 CLASS(physics_base) :: Physics
221 CLASS(mesh_base) :: Mesh
222 CLASS(timedisc_base):: Timedisc
223 !------------------------------------------------------------------------!
224 TYPE(marray_base) :: dv
225 INTEGER, ALLOCATABLE :: seed(:)
226 INTEGER :: i,n,clock
227 CHARACTER(LEN=32) :: info_str
228 !------------------------------------------------------------------------!
229 INTENT(IN) :: mesh,physics
230 INTENT(INOUT) :: timedisc
231 !------------------------------------------------------------------------!
232 ! initialize velocity perturbations
233 dv = marray_base(3)
234 CALL random_seed(size=n)
235 ALLOCATE(seed(n))
236 ! seed the rng with a mix from current time and mpi rank
237 CALL system_clock(count=clock)
238 seed = clock + (timedisc%GetRank()+1) * (/(i-1, i=1,n)/)
239 CALL random_number(dv%data1d)
240 DEALLOCATE(seed)
241
242 ! initial condition
243 SELECT TYPE(pvar => timedisc%pvar)
244 CLASS IS(statevector_euler)
245 pvar%density%data1d(:) = rho0
246 pvar%pressure%data3d(:,:,:) = pin + (pout-pin)/ltube * mesh%bccart(:,:,:,3)
247 pvar%velocity%data1d(:) = (dv%data1d(:)-0.5) * 0.01
248 CLASS DEFAULT
249 CALL timedisc%Error("poiseuille::InitData","physics not supported")
250 END SELECT
251
252 ! fixed boundary conditions
253 ! inflow
254 SELECT TYPE(bbottom => timedisc%Boundary%boundary(bottom)%p)
255 CLASS IS (boundary_fixed)
256 bbottom%data(:,:,:,:) = 0.0
257 bbottom%data(:,:,:,physics%DENSITY) = rho0
258 bbottom%data(:,:,:,physics%PRESSURE) = pin
259 ! imposed density, pressure and tangential velocities;
260 ! extrapolated normal velocity
261 bbottom%fixed(:,:,:) = .true.
262 bbottom%fixed(:,:,physics%ZVELOCITY) = .false.
263 END SELECT
264 ! outflow
265 SELECT TYPE(btop => timedisc%Boundary%boundary(top)%p)
266 CLASS IS (boundary_fixed)
267 btop%data(:,:,:,:) = 0.0
268 btop%data(:,:,:,physics%DENSITY) = rho0
269 btop%data(:,:,:,physics%PRESSURE) = pout
270 ! imposed pressure;
271 ! extrapolated density, tangential and normal velocities
272 btop%fixed(:,:,:) = .false.
273 btop%fixed(:,:,physics%PRESSURE) = .true.
274 END SELECT
275 ! noslip boundary at outer radius
276 SELECT TYPE(b => timedisc%Boundary%boundary(east)%p)
277 CLASS IS(boundary_noslip)
278 b%data(:,:,:,physics%ZVELOCITY) = 0.0
279 CLASS IS(boundary_custom)
280 CALL b%SetCustomBoundaries(mesh,physics, &
281 (/custom_reflect,custom_reflneg,custom_fixed,custom_fixed,custom_reflect/))
282 b%data(:,:,:,:) = 0.0
283 b%data(:,:,:,physics%DENSITY) = rho0 ! actually not used
284 b%data(:,:,:,physics%PRESSURE) = pout ! actually not used
285 END SELECT
286
287 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
288 CALL mesh%Info(" DATA-----> initial condition: " // &
289 "tube with pressure gradient")
290 WRITE (info_str, '(ES10.2)') re
291 CALL physics%Info(" Reynolds number: " // trim(info_str))
292 WRITE (info_str, '(ES10.2)') ma
293 CALL physics%Info(" Mach number: " // trim(info_str))
294 WRITE (info_str, '(ES10.2)') umax
295 CALL physics%Info(" max. laminar velocity: " // trim(info_str) // " m/s")
296 WRITE (info_str, '(ES10.2)') pin
297 CALL physics%Info(" inlet pressure: " // trim(info_str) // " Pa")
298 WRITE (info_str, '(ES10.2)') 1.-pout/pin
299 CALL physics%Info(" rel. pressure gradient: " // trim(info_str))
300 WRITE (info_str, '(ES10.2)') tvis
301 CALL physics%Info(" viscous time scale: " // trim(info_str) // " s")
302
303 END SUBROUTINE initdata
304
305END PROGRAM poiseuille
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
poiseuille
program poiseuille
Program and data initialization for testing the law of Hagen-Poiseuille.
Definition: poiseuille.f90:33
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
Generated by   doxygen 1.9.4