gauss3d.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: gauss3d.f90 #
5 !# #
6 !# Copyright (C) 2006-2014 #
7 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8 !# #
9 !# This program is free software; you can redistribute it and/or modify #
10 !# it under the terms of the GNU General Public License as published by #
11 !# the Free Software Foundation; either version 2 of the License, or (at #
12 !# your option) any later version. #
13 !# #
14 !# This program is distributed in the hope that it will be useful, but #
15 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
16 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17 !# NON INFRINGEMENT. See the GNU General Public License for more #
18 !# details. #
19 !# #
20 !# You should have received a copy of the GNU General Public License #
21 !# along with this program; if not, write to the Free Software #
22 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23 !# #
24 !#############################################################################
25 
26 !----------------------------------------------------------------------------!
31 !----------------------------------------------------------------------------!
32 PROGRAM gauss3d
33  USE fosite_mod
34 #include "tap.h"
35  IMPLICIT NONE
36  !--------------------------------------------------------------------------!
37  ! simulation parameter
38  REAL, PARAMETER :: tsim = 0.5 ! simulation time
39  REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats
40  REAL, PARAMETER :: csiso = &
41  0.0 ! non-isothermal simulation
42 ! 1.127 ! isothermal simulation
43  ! with CSISO as sound speed
44  ! initial condition (dimensionless units)
45  REAL, PARAMETER :: rho0 = 1.0 ! ambient density
46  REAL, PARAMETER :: rho1 = 1.0 ! peak density above RHO0
47  REAL, PARAMETER :: rwidth = 0.06 ! half width of the Gaussian
48  REAL, PARAMETER :: p0 = 1.0 ! ambient pressure
49  REAL, PARAMETER :: p1 = 100.0 ! peak pressure above P0
50  REAL, PARAMETER :: pwidth = 0.06 ! half width of the Gaussian
51  REAL, PARAMETER :: omega0 = 0.0 ! angular velocity
52  REAL, PARAMETER :: eta = 0.0 ! dynamic viscosity (0.0 disables)
53  ! location of the initial pulse (cartesian coordinates)
54  REAL, PARAMETER :: x0 = 0.0 ! x-position
55  REAL, PARAMETER :: y0 = 0.0 ! y-position
56  REAL, PARAMETER :: z0 = 0.0 ! z-position
57  ! mesh settings
58  INTEGER, PARAMETER :: mgeo = cartesian! geometry
59 ! INTEGER, PARAMETER :: MGEO = CYLINDRICAL
60 ! INTEGER, PARAMETER :: MGEO = SPHERICAL
61  INTEGER, PARAMETER :: res = 30 ! x-/y-/z-resolution
62  REAL, PARAMETER :: rmax = 0.5 ! maximal radial extent
63  REAL, PARAMETER :: gpar = 0.8 ! geometry scaling parameter
64  ! output parameters
65  INTEGER, PARAMETER :: onum = 10 ! number of output data sets
66  CHARACTER(LEN=256), PARAMETER & ! output data dir
67  :: odir = './'
68  CHARACTER(LEN=256), PARAMETER & ! output data file name
69  :: ofname = 'gauss3D'
70  !--------------------------------------------------------------------------!
71  CLASS(fosite), ALLOCATABLE :: sim
72  LOGICAL :: ok
73  !--------------------------------------------------------------------------!
74 
75  tap_plan(1)
76 
77  ALLOCATE(sim)
78  CALL sim%InitFosite()
79  CALL makeconfig(sim, sim%config)
80  CALL sim%Setup()
81  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
82 
83  CALL sim%Run()
84  ok = .NOT.sim%aborted
85  CALL sim%Finalize()
86  DEALLOCATE(sim)
87 
88  tap_check(ok,"stoptime reached")
89  tap_done
90 
91 CONTAINS
92 
93  SUBROUTINE makeconfig(Sim, config)
94  USE functions, ONLY : asinh, acosh
95  IMPLICIT NONE
96  !------------------------------------------------------------------------!
97  CLASS(fosite) :: Sim
98  TYPE(Dict_TYP),POINTER :: config
99  !------------------------------------------------------------------------!
100  ! Local variable declaration
101  INTEGER :: bc(6),xres,yres,zres
102  TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, sources, &
103  timedisc, fluxes
104  REAL :: x1,x2,y1,y2,z1,z2
105  !------------------------------------------------------------------------!
106  INTENT(INOUT) :: sim
107  !------------------------------------------------------------------------!
108  ! mesh settings and boundary conditions
109  SELECT CASE(mgeo)
110  CASE(cartesian)
111  x1 = -rmax
112  x2 = rmax
113  y1 = -rmax
114  y2 = rmax
115  z1 = -rmax
116  z2 = rmax
117  xres = res
118  yres = res
119  zres = res
120 ! bc(WEST) = NO_GRADIENTS
121 ! bc(EAST) = NO_GRADIENTS
122 ! bc(SOUTH) = NO_GRADIENTS
123 ! bc(NORTH) = NO_GRADIENTS
124 ! bc(BOTTOM) = NO_GRADIENTS
125 ! bc(TOP) = NO_GRADIENTS
126  bc(west) = absorbing
127  bc(east) = absorbing
128  bc(south) = absorbing
129  bc(north) = absorbing
130  bc(bottom) = absorbing
131  bc(top) = absorbing
132  CASE(cylindrical)
133  x1 = 0.01
134  x2 = rmax
135  y1 = 0.0
136  y2 = 2*pi
137  z1 = -rmax
138  z2 = rmax
139  xres = res / 2
140  yres = res / 2
141  zres = res
142  bc(west) = axis
143  bc(east) = absorbing
144  bc(south) = periodic
145  bc(north) = periodic
146  bc(bottom) = absorbing
147  bc(top) = absorbing
148  CASE(spherical)
149  x1 = 0.01
150  x2 = rmax
151  y1 = 0.0
152  y2 = pi
153  z1 = 0.0
154  z2 = 2*pi
155  xres = res / 2
156  yres = res / 4
157  zres = res / 2
158  bc(west) = axis
159  bc(east) = absorbing
160  bc(south) = axis
161  bc(north) = axis
162  bc(bottom)= periodic
163  bc(top) = periodic
164  CASE DEFAULT
165  CALL error(sim%Physics,"InitProgram","geometry not supported for this test")
166  END SELECT
167 
168  !mesh settings
169  mesh => dict( &
170  "meshtype" / midpoint, &
171  "geometry" / mgeo, &
172  "inum" / xres, &
173  "jnum" / yres, &
174  "knum" / zres, &
175  "xmin" / x1, &
176  "xmax" / x2, &
177  "ymin" / y1, &
178  "ymax" / y2, &
179  "zmin" / z1, &
180  "zmax" / z2, &
181 ! "decomposition"/ (/ 1, 1, 2/), &
182  "output/radius" / 1, &
183  "gparam" / gpar)
184 
185  ! boundary conditions
186  boundary => dict( &
187  "western" / bc(west), &
188  "eastern" / bc(east), &
189  "southern" / bc(south), &
190  "northern" / bc(north), &
191  "bottomer" / bc(bottom), &
192  "topper" / bc(top))
193 
194  ! physics settings
195  IF (csiso.GT.tiny(csiso)) THEN
196  physics => dict("problem" / euler_isotherm, &
197  "cs" / csiso) ! isothermal speed of sound
198  ELSE
199  physics => dict("problem" / euler, &
200  "gamma" / gamma) ! ratio of specific heats
201  END IF
202 
203  ! flux calculation and reconstruction method
204  fluxes => dict( &
205  "order" / linear, &
206  "fluxtype" / kt, &
207  "variables" / primitive, &
208  "limiter" / vanleer, &
209  "output/slopes" / 0)
210 
211  NULLIFY(sources)
212 
213  ! time discretization settings
214  timedisc => dict( &
215  "method" / modified_euler, &
216  "order" / 3, &
217  "tol_rel" / 10.0, & ! > 1 disables adaptive step size
218  "cfl" / 0.4, &
219  "stoptime" / tsim, &
220  "dtlimit" / 1.0e-8, &
221  "maxiter" / 10000000)
222 
223  datafile => dict( &
224  "fileformat" / vtk, &
225  "filename" / (trim(odir) // trim(ofname)), &
226  "count" / onum)
227 
228  config => dict( &
229  "mesh" / mesh, &
230  "physics" / physics, &
231  "boundary" / boundary, &
232  "fluxes" / fluxes, &
233  "datafile" / datafile, &
234  "timedisc" / timedisc)
235 
236  IF (ASSOCIATED(sources)) &
237  CALL setattr(config, "sources", sources)
238 
239  END SUBROUTINE makeconfig
240 
241 
242  SUBROUTINE initdata(Mesh,Physics,Timedisc)
243  USE geometry_base_mod
244  IMPLICIT NONE
245  !------------------------------------------------------------------------!
246  CLASS(physics_base) :: Physics
247  CLASS(mesh_base) :: Mesh
248  CLASS(timedisc_base) :: Timedisc
249  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
250  :: radius
251  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3) &
252  :: posvec
253  !------------------------------------------------------------------------!
254  INTENT(IN) :: mesh,physics
255  INTENT(INOUT) :: timedisc
256  !------------------------------------------------------------------------!
257  IF (any((/abs(x0),abs(y0),abs(z0)/).GT.tiny(z0))) THEN
258  ! shifted point mass position:
259  ! compute curvilinear components of shift vector
260  posvec(:,:,:,1) = x0
261  posvec(:,:,:,2) = y0
262  posvec(:,:,:,3) = z0
263  CALL mesh%geometry%Convert2Curvilinear(mesh%bcenter,posvec,posvec)
264  ! subtract the result from the position vector:
265  ! this gives you the curvilinear components of all vectors pointing
266  ! from the point mass to the bary center of any cell on the mesh
267  posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - posvec(:,:,:,:)
268  ! compute its absolute value
269  radius(:,:,:) = sqrt(posvec(:,:,:,1)**2+posvec(:,:,:,2)**2+posvec(:,:,:,3)**2)
270  ELSE
271  ! no shift of point mass set radius and posvec to Mesh defaults
272  radius(:,:,:) = mesh%radius%bcenter(:,:,:)
273  posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
274  END IF
275 
276  ! initial condition
277  SELECT TYPE(pvar => timedisc%pvar)
278  TYPE IS(statevector_eulerisotherm) ! isothermal HD
279  ! Gaussian density pulse
280  pvar%density%data3d(:,:,:) = rho0 + rho1*exp(-log(2.0) &
281  * (radius(:,:,:)/rwidth)**2)
282  ! vanishing velocities
283  pvar%velocity%data1d(:) = 0.0
284  TYPE IS(statevector_euler) ! non-isothermal HD
285  ! constant density
286  pvar%density%data1d(:) = rho0
287  ! Gaussian pressure pulse
288  pvar%pressure%data3d(:,:,:) = p0 + p1*exp(-log(2.0) &
289  * (radius(:,:,:)/rwidth)**2)
290  ! vanishing velocities
291  pvar%velocity%data1d(:) = 0.0
292  CLASS DEFAULT
293  CALL physics%Error("gauss3d::InitData","unsupported physics")
294  END SELECT
295 
296  CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
297 
298  CALL mesh%Info(" DATA-----> initial condition: 3D Gaussian pulse")
299 
300  END SUBROUTINE initdata
301 
302 END PROGRAM gauss3d
functions
Mathematical functions.
Definition: functions.f90:48
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
gauss3d
program gauss3d
3D Gaussian pressure or density pulse with and without rotation
Definition: gauss3d.f90:32
functions::asinh
elemental real function, public asinh(x)
inverse hyperbolic sine function
Definition: functions.f90:100
geometry_base_mod
base class for geometrical properties
Definition: geometry_base.f90:36
makeconfig
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
fosite_mod
Definition: fosite.f90:41
functions::acosh
elemental real function, public acosh(x)
inverse hyperbolic cosine function
Definition: functions.f90:110
Generated by   doxygen 1.8.14