sblintheo3D.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: sblintheo.f90 #
5 !# #
6 !# Copyright (C) 2015 #
7 !# Jannes Klee <jklee@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 !#############################################################################
26 !----------------------------------------------------------------------------!
71 !----------------------------------------------------------------------------!
72 PROGRAM sblintheo
73  USE fosite_mod
74 #include "tap.h"
75  IMPLICIT NONE
76  !--------------------------------------------------------------------------!
77  ! general constants
78  REAL, PARAMETER :: gn = 1.0 ! grav. constant [GEOM] !
79  ! make a shearingsheet simulation and denote the direction of shearing
80  ! (applies boundaries, ficitious forces and fargo automatically)
81  INTEGER, PARAMETER :: shearsheet_direction = 1
82  ! simulation parameter
83  REAL, PARAMETER :: omega = 1.0 ! rotation at fid. point !
84  REAL, PARAMETER :: sigma0 = 1.0 ! mean surf.dens. !
85  REAL, PARAMETER :: delsigma = sigma0*5.0e-4 ! disturbance !
86  REAL, PARAMETER :: tsim = 10.0/omega ! simulation time !
87  REAL, PARAMETER :: gamma = 2.0 ! dep. on vert. struct. !
88  REAL, PARAMETER :: q = 1.5 ! shearing parameter !
89  ! set (initial) speed of sound using the Toomre criterion and
90  ! assuming marginal stabily; in non-isothermal simulations this is
91  ! used to set the initial pressure (see InitData) whereas in isothermal
92  ! simulations this is constant throughout the simulation
93  REAL, PARAMETER :: soundspeed = pi*gn*sigma0/omega ! det. by. Toomre !
94  ! mesh settings
95  INTEGER, PARAMETER :: mgeo = cartesian
96  INTEGER, PARAMETER :: xres = 30 ! amount of cells in x- !
97  INTEGER, PARAMETER :: yres = 30 ! y-direction (rho/phi) !
98  INTEGER, PARAMETER :: zres = 30
99  REAL :: domainx = 40.0 ! domain size [GEOM] !
100  REAL :: domainy = 40.0 ! domain size [GEOM] !
101  REAL :: domainz = 40.0 ! domain size [GEOM] !
102  ! number of output time steps
103 ! INTEGER, PARAMETER :: ONUM = 1
104  INTEGER, PARAMETER :: onum = 10
105  ! output directory and output name
106  CHARACTER(LEN=256), PARAMETER :: odir = "./"
107  CHARACTER(LEN=256), PARAMETER :: ofname = "sblintheo"
108  !--------------------------------------------------------------------------!
109  CLASS(fosite), ALLOCATABLE :: sim
110  !--------------------------------------------------------------------------!
111 
112 !TAP_PLAN(1)
113 
114 ALLOCATE(sim)
115 
116 CALL sim%InitFosite()
117 CALL makeconfig(sim, sim%config)
118 CALL sim%Setup()
119 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc%pvar, sim%Timedisc%cvar)
120 CALL sim%Run()
121 
122 ! search for the amplitude
123 !maximum = MAXVAL(Sim%Timedisc%pvar(:,:,:,Sim%Physics%DENSITY)) - SIGMA0
124 !TAP_CHECK_CLOSE(maximum, 0.0316463969505, 0.005, "Last max. < 0.005 deviation")
125 
126 CALL sim%Finalize()
127 DEALLOCATE(sim)
128 
129 !TAP_DONE
130 
131 CONTAINS
132 
134  SUBROUTINE makeconfig(Sim,config)
135  IMPLICIT NONE
136  !--------------------------------------------------------------------------!
137  CLASS(fosite) :: Sim
138  TYPE(Dict_TYP), POINTER :: config
139  !--------------------------------------------------------------------------!
140  TYPE(Dict_TYP), POINTER :: mesh,physics,fluxes,&
141  grav,sources,timedisc,shearingbox,&
142  datafile
143  REAL :: XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
144  !--------------------------------------------------------------------------!
145  domainx = domainx*gn*sigma0/(omega*omega)
146  domainy = domainy*gn*sigma0/(omega*omega)
147  xmin = -0.5*domainx
148  xmax = +0.5*domainx
149  ymin = -0.5*domainy
150  ymax = +0.5*domainy
151  zmin = -0.5*domainz
152  zmax = +0.5*domainz
153 
154  ! physics settings
155  physics => dict(&
156  "problem" / euler_isotherm, &
157  "cs" / soundspeed, &
158 ! "problem" / EULER3D, &
159 ! "gamma" / GAMMA, &
160  "units" / geometrical &
161  )
162 
163  ! mesh settings
164  mesh => dict(&
165  "meshtype" / midpoint, &
166  "geometry" / mgeo, &
167  "shear_dir" / shearsheet_direction, &
168  "inum" / xres, &
169  "jnum" / yres, &
170  "knum" / zres, &
171  "xmin" / xmin, &
172  "xmax" / xmax, &
173  "ymin" / ymin, &
174  "ymax" / ymax, &
175  "zmin" / zmin, &
176  "zmax" / zmax, &
177  "omega" / omega, &
178  "decomposition"/ (/ 1, -1, 1/), &
179  "output/rotation" / 0, &
180  "output/volume" / 0, &
181  "output/bh" / 0, &
182  "output/dl" / 0, &
183  "Q" / q &
184  )
185 
186  ! fluxes settings
187  fluxes => dict(&
188  "order" / linear, &
189  "fluxtype" / kt, &
190  "variables" / primitive, &
191  "limiter" / vanleer &
192  )
193 
194  ! gravity settings (source term)
195  grav => dict(&
196  "stype" / gravity, &
197  "self/gtype" / sboxspectral, &
198 ! "output/accel" / 1, &
199  "self/output/phi" / 0, &
200  "self/output/accel_x" / 0, &
201  "self/output/accel_y" / 0 &
202  )
203 
204  ! shearing box fictious forces
205  shearingbox => dict(&
206  "stype" / shearbox &
207  )
208 
209  ! sources settings (contains source terms)
210  sources => dict(&
211 ! "grav" / grav, &
212  "shearing" / shearingbox &
213  )
214 
215  ! time discretization settings
216  timedisc => dict( &
217  "method" / dormand_prince, &
218  "cfl" / 0.4, &
219  "stoptime" / tsim, &
220  "dtlimit" / 1.0e-40, &
221  "tol_rel" / 0.1, &
222  "maxiter" / 100000000)
223 
224  ! initialize data input/output
225  datafile => dict( &
226  "fileformat" / vtk, &
227  "filename" / (trim(odir) // trim(ofname)), &
228  "count" / onum)
229 
230  ! overall config settings
231  config => dict(&
232  "mesh" / mesh, &
233  "physics" / physics, &
234  "fluxes" / fluxes, &
235  "sources" / sources, &
236  "timedisc" / timedisc, &
237  "datafile" / datafile &
238  )
239  END SUBROUTINE makeconfig
240 
242  SUBROUTINE initdata(Mesh,Physics,pvar,cvar)
243  IMPLICIT NONE
244  !------------------------------------------------------------------------!
245  CLASS(mesh_base), INTENT(IN) :: Mesh
246  CLASS(physics_base), INTENT(IN) :: Physics
247  CLASS(marray_compound), POINTER, INTENT(INOUT) :: pvar,cvar
248  !------------------------------------------------------------------------!
249  ! local variable declaration
250  REAL :: kx, ky, kz
251  !------------------------------------------------------------------------!
252 
253  !---------------------- linear theory test ------------------------------!
254  kx = -2*(2*pi/domainx)
255  ky = 2*pi/domainy
256  kz = 2*pi/domainz
257 
258  ! initial condition
259  SELECT TYPE(p => pvar)
260  TYPE IS(statevector_eulerisotherm)
261  IF(mesh%shear_dir.EQ.2)THEN
262  p%density%data3d(:,:,:) = sigma0 &
263  + delsigma*cos(kx*mesh%bcenter(:,:,:,1) &
264  + ky*mesh%bcenter(:,:,:,2) &
265  + kz*mesh%bcenter(:,:,:,3))
266  p%velocity%data2d(:,1) = 0.0
267  p%velocity%data4d(:,:,:,2) = -q*omega*mesh%bcenter(:,:,:,1)
268  p%velocity%data2d(:,3) = 0.0
269  ELSE IF(mesh%shear_dir.EQ.1)THEN
270  p%density%data3d(:,:,:) = sigma0 &
271  + delsigma*cos(kx*mesh%bcenter(:,:,:,2) &
272  - ky*mesh%bcenter(:,:,:,1) &
273  + kz*mesh%bcenter(:,:,:,3))
274  p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
275  p%velocity%data2d(:,2) = 0.0
276  p%velocity%data2d(:,3) = 0.0
277  END IF
278  TYPE IS(statevector_euler) ! non-isothermal HD
279  IF(mesh%shear_dir.EQ.2)THEN
280  p%density%data3d(:,:,:) = sigma0 &
281  + delsigma*cos(kx*mesh%bcenter(:,:,:,1) &
282  + ky*mesh%bcenter(:,:,:,2) &
283  + kz*mesh%bcenter(:,:,:,3))
284  p%velocity%data2d(:,1) = 0.0
285  p%velocity%data4d(:,:,:,2) = -q*omega*mesh%bcenter(:,:,:,1)
286  p%velocity%data2d(:,3) = 0.0
287  ELSE IF(mesh%shear_dir.EQ.1)THEN
288  p%density%data3d(:,:,:) = sigma0 &
289  + delsigma*cos(kx*mesh%bcenter(:,:,:,2) &
290  - ky*mesh%bcenter(:,:,:,1) &
291  + kz*mesh%bcenter(:,:,:,3))
292  p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
293  p%velocity%data2d(:,2) = 0.0
294  p%velocity%data2d(:,3) = 0.0
295  END IF
296  CLASS DEFAULT
297  CALL physics%Error("shear::InitData","only (non-)isothermal HD supported")
298  END SELECT
299 
300  CALL physics%Convert2Conservative(pvar,cvar)
301  CALL mesh%Info(" DATA-----> initial condition: " // &
302  "Linear Theory Test - Shearingsheet")
303  END SUBROUTINE initdata
304 END PROGRAM sblintheo
omega
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
Definition: sources_viscosity.f90:436
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
makeconfig
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
sblintheo
program sblintheo
Definition: sblintheo.f90:72
fosite_mod
Definition: fosite.f90:41
Generated by   doxygen 1.8.14