sblintheo.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-2018 #
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 = 4.54845485/omega ! simulation time (first maximum) !
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
94  ! mesh settings
95  INTEGER, PARAMETER :: mgeo = cartesian
96  INTEGER, PARAMETER :: xres = 128 ! amount of cells in x- !
97  INTEGER, PARAMETER :: yres = 128 ! y-direction (rho/phi) !
98  INTEGER, PARAMETER :: zres = 1
99  REAL :: domainx = 40.0 ! domain size [GEOM] !
100  REAL :: domainy = 40.0 ! domain size [GEOM] !
101  ! number of output time steps
102  INTEGER, PARAMETER :: onum = 10
103  ! output directory and output name
104  CHARACTER(LEN=256), PARAMETER :: odir = "./"
105  CHARACTER(LEN=256), PARAMETER :: ofname = "sblintheo"
106  !--------------------------------------------------------------------------!
107  CLASS(fosite), ALLOCATABLE :: sim
108  REAL :: maximum
109  !--------------------------------------------------------------------------!
110 
111  tap_plan(1)
112 
113  ALLOCATE(sim)
114 
115  CALL sim%InitFosite()
116  CALL makeconfig(sim, sim%config)
117  CALL sim%Setup()
118  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc%pvar, sim%Timedisc%cvar)
119  CALL sim%Run()
120 
121  ! check amplitude
122  maximum = maxval(sim%Timedisc%pvar%data4d(:,:,:,sim%Physics%DENSITY)) - sigma0
123 
124  CALL sim%Finalize()
125  DEALLOCATE(sim)
126 
127  tap_check_close(maximum, 0.02395931, 0.0003, "First max. < 3e-4 deviation")
128  tap_done
129 
130 CONTAINS
131 
133  SUBROUTINE makeconfig(Sim,config)
134  IMPLICIT NONE
135  !--------------------------------------------------------------------------!
136  CLASS(fosite) :: Sim
137  TYPE(Dict_TYP), POINTER :: config
138  !--------------------------------------------------------------------------!
139  TYPE(Dict_TYP), POINTER :: mesh,physics,fluxes,&
140  grav,sources,timedisc,shearingbox,&
141  datafile
142  REAL :: XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
143  !--------------------------------------------------------------------------!
144  domainx = domainx*gn*sigma0/(omega*omega)
145  domainy = domainy*gn*sigma0/(omega*omega)
146  xmin = -0.5*domainx
147  xmax = +0.5*domainx
148  ymin = -0.5*domainy
149  ymax = +0.5*domainy
150  zmin = 0.0
151  zmax = 0.0
152 
153  ! physics settings
154  physics => dict(&
155  "problem" / euler_isotherm, &
156  "cs" / soundspeed, &
157  "units" / geometrical &
158  )
159 
160  ! mesh settings
161  mesh => dict(&
162  "meshtype" / midpoint, &
163  "geometry" / mgeo, &
164  "shear_dir" / shearsheet_direction, &
165  "inum" / xres, &
166  "jnum" / yres, &
167  "knum" / zres, &
168  "xmin" / xmin, &
169  "xmax" / xmax, &
170  "ymin" / ymin, &
171  "ymax" / ymax, &
172  "zmin" / zmin, &
173  "zmax" / zmax, &
174  "omega" / omega, &
175  "decomposition"/ (/ 1, -1, 1/), &
176  "output/rotation" / 0, &
177  "output/volume" / 0, &
178  "output/bh" / 0, &
179  "output/dl" / 0, &
180  "Q" / q &
181  )
182 
183  ! fluxes settings
184  fluxes => dict(&
185  "order" / linear, &
186  "fluxtype" / kt, &
187  "variables" / primitive, &
188  "limiter" / vanleer &
189  )
190 
191  ! gravity settings (source term)
192  grav => dict(&
193  "stype" / gravity, &
194  "self/gtype" / sboxspectral, &
195 ! "output/accel" / 1, &
196  "self/output/phi" / 0, &
197  "self/output/accel_x" / 0, &
198  "self/output/accel_y" / 0 &
199  )
200 
201  ! shearing box fictious forces
202  shearingbox => dict(&
203  "stype" / shearbox &
204  )
205 
206  ! sources settings (contains source terms)
207  sources => dict(&
208  "shearing" / shearingbox, &
209  "grav" / grav &
210  )
211 
212  ! time discretization settings
213  timedisc => dict( &
214  "method" / dormand_prince, &
215  "cfl" / 0.4, &
216  "stoptime" / tsim, &
217  "dtlimit" / 1.0e-40, &
218  "tol_rel" / 0.1, &
219  "maxiter" / 100000000)
220 
221  ! initialize data input/output
222  datafile => dict( &
223  "fileformat" / vtk, &
224  "filename" / (trim(odir) // trim(ofname)), &
225  "count" / onum)
226 
227  ! overall config settings
228  config => dict(&
229  "mesh" / mesh, &
230  "physics" / physics, &
231  "fluxes" / fluxes, &
232  "sources" / sources, &
233  "timedisc" / timedisc, &
234  "datafile" / datafile &
235  )
236  END SUBROUTINE makeconfig
237 
240  SUBROUTINE initdata(Mesh,Physics,pvar,cvar)
241  IMPLICIT NONE
242  !------------------------------------------------------------------------!
243  CLASS(mesh_base), INTENT(IN) :: Mesh
244  CLASS(physics_base), INTENT(IN) :: Physics
245  CLASS(marray_compound), POINTER, INTENT(INOUT) :: pvar,cvar
246  !------------------------------------------------------------------------!
247  ! local variable declaration
248  REAL :: kx, ky
249  !------------------------------------------------------------------------!
250 
251  !---------------------- linear theory test ------------------------------!
252  kx = -2*(2*pi/domainx)
253  ky = 2*pi/domainy
254 
255  ! initial condition
256  SELECT TYPE(p => pvar)
257  TYPE IS(statevector_eulerisotherm)
258  IF(mesh%shear_dir.EQ.2)THEN
259  p%density%data3d(:,:,:) = sigma0 &
260  + delsigma*cos(kx*mesh%bcenter(:,:,:,1) + ky*mesh%bcenter(:,:,:,2))
261  p%velocity%data2d(:,1) = 0.0
262  p%velocity%data4d(:,:,:,2) = -q*omega*mesh%bcenter(:,:,:,1)
263  ELSE IF(mesh%shear_dir.EQ.1)THEN
264  p%density%data3d(:,:,:) = sigma0 &
265  + delsigma*cos(kx*mesh%bcenter(:,:,:,2) - ky*mesh%bcenter(:,:,:,1))
266  p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
267  p%velocity%data2d(:,2) = 0.0
268  END IF
269  TYPE IS(statevector_euler) ! non-isothermal HD
270  IF(mesh%shear_dir.EQ.2)THEN
271  p%density%data3d(:,:,:) = sigma0 &
272  + delsigma*cos(kx*mesh%bcenter(:,:,:,1) + ky*mesh%bcenter(:,:,:,2))
273  p%velocity%data2d(:,1) = 0.0
274  p%velocity%data4d(:,:,:,2) = -q*omega*mesh%bcenter(:,:,:,1)
275  ELSE IF(mesh%shear_dir.EQ.1)THEN
276  p%density%data3d(:,:,:) = sigma0 &
277  + delsigma*cos(kx*mesh%bcenter(:,:,:,2) - ky*mesh%bcenter(:,:,:,1))
278  p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
279  p%velocity%data2d(:,2) = 0.0
280  END IF
281  p%pressure%data3d(:,:,:) = soundspeed**2 * sigma0 / gamma
282  CLASS DEFAULT
283  CALL physics%Error("shear::InitData","only (non-)isothermal HD supported")
284  END SELECT
285 
286  CALL physics%Convert2Conservative(pvar,cvar)
287 
288  !------------------------------------------------------------------------!
289 
290  CALL mesh%Info(" DATA-----> initial condition: " // &
291  "Linear Theory Test - Shearingsheet")
292  END SUBROUTINE initdata
293 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