keplerianvortex.f90

Vortex transport in a keplerian disk
References:

  • [mignone2012] Mignone, A. et al.: "A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code" Astronomy & Astrophysics, Volume 545, id.A152, 16 pp.
  • [bodo2007] Bodo, G. et al.: "Stability amd nonlinear adjustment of vortices in Keplerian flows" A & A 475, 51-61 (2007)
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: keplerianvortex.f90 #
5 !# #
6 !# Copyright (C) 2014-2018 #
7 !# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
8 !# Jannes Klee <jklee@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 !----------------------------------------------------------------------------!
39 !----------------------------------------------------------------------------!
40 PROGRAM keplerianvortex
41  USE fosite_mod
42  IMPLICIT NONE
43  !--------------------------------------------------------------------------!
44  ! general constants
45  REAL, PARAMETER :: GN = 1.0
46  REAL, PARAMETER :: MCENTRAL = 1.0
47  REAL, PARAMETER :: Sigma0 = 1.0
48  REAL, PARAMETER :: M = 10.
49  REAL, PARAMETER :: omega = 0.
50 
51  REAL, PARAMETER :: P0 = 2.0*pi/1.0
52  REAL, PARAMETER :: P = 3
53  REAL, PARAMETER :: t = p*p0
54 
55  REAL, PARAMETER :: GAMMA = 5./3.
56 
57  REAL, PARAMETER :: RMIN = 0.5
58  REAL, PARAMETER :: RMAX = 1.5
59 
60  INTEGER, PARAMETER :: GravEnergySource = 1
61  INTEGER, PARAMETER :: XRES = 64
62  INTEGER, PARAMETER :: YRES = 4*xres
63  INTEGER, PARAMETER :: ONUM = p * 10
64  !--------------------------------------------------------------------------!
65  CLASS(fosite), ALLOCATABLE :: Sim
66  !--------------------------------------------------------------------------!
67 
68 ALLOCATE(sim)
69 CALL sim%InitFosite()
70 CALL makeconfig(sim%config)
71 CALL sim%Setup()
72 CALL initdata(sim%Timedisc, sim%Mesh, sim%Physics, sim%Fluxes)
73 CALL sim%Run()
74 CALL sim%Finalize()
75 DEALLOCATE(sim)
76 
77 CONTAINS
78 
79  SUBROUTINE makeconfig(config)
80  IMPLICIT NONE
81  !------------------------------------------------------------------------!
82  TYPE(Dict_TYP),POINTER :: config
83  !------------------------------------------------------------------------!
84  ! Local variable declaration
85  TYPE(Dict_TYP),POINTER :: mesh,boundary,timedisc,datafile,&
86  sources,fluxes,pmass,physics,damping,vis,pot
87  !------------------------------------------------------------------------!
88  physics => dict( &
89  "problem" / euler, &
90  "gamma" / gamma, &
91  "units" / geometrical)
92 
93  fluxes => dict( &
94  "fluxtype" / kt, &
95  "order" / linear, &
96  "variables" / primitive, &
97  "limiter" / vanleer)
98 
99  mesh => dict( &
100  "meshtype" / midpoint, &
101  "geometry" / cylindrical, &
102  "omega" / omega, &
103  "inum" / xres, &
104  "jnum" / yres, &
105  "knum" / 1, &
106  "xmin" / rmin, &
107  "xmax" / rmax, &
108  "ymin" / (-pi), &
109  "ymax" / (pi), &
110  "zmin" / 0.0, &
111  "zmax" / 0.0, &
112  "decomposition" / (/ 1, -1/), &
113  "output/volume" / 1,&
114  "output/bh" / 1,&
115  "output/rotation" / 1,&
116  "output/dl" / 1)
117 
118  boundary => dict( &
119  "western" / custom, &
120  "eastern" / custom, &
121  "southern" / periodic, &
122  "northern" / periodic)
123 
124  pmass => dict( &
125  "gtype" / pointmass, &
126  "potential" / newton, &
127  "mass" / mcentral, &
128  "cvis" / 0.9, &
129  "outbound" / 0)!, &
130 ! "output/accel" / 1 )
131 
132 ! damping => Dict( &
133 ! "stype" / WAVE_DAMPING, &
134 ! "zone0/rmin" / RMIN, &
135 ! "zone0/rmax" / (0.5), &
136 ! "zone0/tau" / (0.1*2.0*PI / SQRT(GN*MCENTRAL/(0.4)**3)), &
137 ! "zone1/rmin" / (1.9), &
138 ! "zone1/rmax" / RMAX, &
139 ! "zone1/tau" / (0.1*2.0*PI / SQRT(GN*MCENTRAL/(2.0)**3)))
140 
141  sources => dict( &
142 ! "damping" / damping, &
143  "grav/stype" / gravity, &
144 ! "grav/energy" / GravEnergySource, &
145  "grav/output/accel" / 1, &
146  "grav/pmass" / pmass)
147 
148  timedisc => dict( &
149  "method" / ssprk, &
150 ! "fargo" / 1, &
151  "tol_rel" / 1.0e-3, &
152  "tol_abs" / (/ 0., 1.e-3, 1.e-3, 0. /), &
153  "cfl" / 0.4, &
154  "stoptime" / t, &
155  "dtlimit" / 1.0e-10, &
156  "maxiter" / 2000000000, &
157  "output/xmomentum" / 1, &
158  "output/ymomentum" / 1, &
159  "output/energy" / 1, &
160 ! "rhstype" / 1, &
161  "output/geometrical_sources" / 1, &
162  "output/external_sources" / 1, &
163  !"output/bflux" / 1, &
164  "output/fluxes" / 1, &
165  "output/rhs" / 1)
166 
167 
168  datafile => dict( &
169  "fileformat" / xdmf, &
170  "filename" / "keplerianvortex", &
171  "count" / onum)
172 
173  config => dict( &
174  "physics" / physics, &
175  "fluxes" / fluxes, &
176  "mesh" / mesh, &
177  "boundary" / boundary, &
178  "sources" / sources, &
179  "timedisc" / timedisc, &
180  "datafile" / datafile)
181 
182  END SUBROUTINE makeconfig
183 
184 
185  SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes)
186  IMPLICIT NONE
187  !------------------------------------------------------------------------!
188  CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
189  CLASS(mesh_base), INTENT(IN) :: Mesh
190  CLASS(physics_base), INTENT(INOUT) :: Physics
191  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
192  !------------------------------------------------------------------------!
193  ! Local variable declaration
194  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
195  :: dvr, dvphi, x, y, phi, r
196  REAL :: h, kappa
197  INTEGER :: dir,j,ig,i,k
198  CLASS(sources_base), POINTER :: sp
199  CLASS(sources_gravity), POINTER :: gp
200  !------------------------------------------------------------------------!
201  x = mesh%bccart(:,:,:,1)
202  y = mesh%bccart(:,:,:,2)
203 
204  IF(gravenergysource.EQ.0) THEN
205  CALL setpotential_euler2dia(physics,mesh,-gn * mcentral / mesh%radius%bcenter)
206  CALL setpotential_euler2dia(physics,mesh,-gn * mcentral / mesh%radius%faces)
207  END IF
208 
209  ! get gravitational acceleration
210  sp => sim%Sources
211  DO
212  IF (ASSOCIATED(sp).EQV..false.) RETURN
213  SELECT TYPE(sp)
214  CLASS IS(sources_gravity)
215  gp => sp
216  EXIT
217  END SELECT
218  sp => sp%next
219  END DO
220 
221  SELECT TYPE (pvar => timedisc%pvar)
222  CLASS IS(statevector_euler)
223  SELECT TYPE (cvar => timedisc%cvar)
224  CLASS IS(statevector_euler)
225  cvar%density%data3d(:,:,:) = sigma0
226  cvar%momentum%data4d(:,:,:,1:physics%VDIM) = 0.
227 
228  CALL physics%Convert2Primitive(cvar,pvar)
229 
230  pvar%pressure%data3d(:,:,:) = 1./(gamma*m**2)
231 
232  pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
233  pvar%velocity%data4d(:,:,:,1:physics%VDIM) &
234  + timedisc%GetCentrifugalVelocity(mesh,physics,sim%Fluxes,sim%Sources,(/0.,0.,1./))
235 
236  END SELECT
237  END SELECT
238 
239  ! boundary conditions
240  ! custom boundary conditions at western boundary if requested
241  SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
242  CLASS IS (boundary_custom)
243  CALL bwest%SetCustomBoundaries(mesh,physics, &
244  (/custom_reflect,custom_reflneg,custom_kepler,custom_reflect/))
245  END SELECT
246  SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
247  CLASS IS (boundary_custom)
248  CALL beast%SetCustomBoundaries(mesh,physics, &
249  (/custom_reflect,custom_reflneg,custom_kepler,custom_reflect/))
250  END SELECT
251 
252  h = 1./(2.*m)
253  kappa = -1.
254  phi = mesh%center(:,:,:,2)
255  x = x - 1./sqrt(2.)
256  y = y - 1./sqrt(2.)
257  r = sqrt(x**2 + y**2)
258 
259  dvr = kappa * exp(-r**2/h**2) * (-y*cos(phi) + x*sin(phi))
260  dvphi = kappa * exp(-r**2/h**2) * ( y*sin(phi) + x*cos(phi))
261 
262  SELECT TYPE (pvar => timedisc%pvar)
263  CLASS IS(statevector_euler)
264  SELECT TYPE (cvar => timedisc%cvar)
265  CLASS IS(statevector_euler)
266  pvar%velocity%data4d(:,:,:,1) = pvar%velocity%data4d(:,:,:,1) + dvr
267  pvar%velocity%data4d(:,:,:,2) = pvar%velocity%data4d(:,:,:,2) + dvphi
268  END SELECT
269  END SELECT
270 
271  CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
272 
273  END SUBROUTINE initdata
274 END PROGRAM keplerianvortex
Generated by   doxygen 1.8.14