keplerianvortex.f90

Vortex transport in a keplerian disk
.

Vortex transport in a keplerian disk
References:

  • [34] 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.
  • [3] 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!----------------------------------------------------------------------------!
40PROGRAM 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
68ALLOCATE(sim)
69CALL sim%InitFosite()
70CALL makeconfig(sim%config)
71CALL sim%Setup()
72CALL initdata(sim%Timedisc, sim%Mesh, sim%Physics, sim%Fluxes)
73CALL sim%Run()
74CALL sim%Finalize()
75DEALLOCATE(sim)
76
77CONTAINS
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" / 0,&
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 "bottomer" / reflecting, &
124 "topper" / reflecting)
125
126 pmass => dict( &
127 "gtype" / pointmass, &
128 "potential" / newton, &
129 "mass" / mcentral, &
130 "cvis" / 0.9, &
131 "outbound" / 0)!, &
132! "output/accel" / 1 )
133
134! damping => Dict( &
135! "stype" / WAVE_DAMPING, &
136! "zone0/rmin" / RMIN, &
137! "zone0/rmax" / (0.5), &
138! "zone0/tau" / (0.1*2.0*PI / SQRT(GN*MCENTRAL/(0.4)**3)), &
139! "zone1/rmin" / (1.9), &
140! "zone1/rmax" / RMAX, &
141! "zone1/tau" / (0.1*2.0*PI / SQRT(GN*MCENTRAL/(2.0)**3)))
142
143 sources => dict( &
144! "damping" / damping, &
145 "grav/stype" / gravity, &
146! "grav/energy" / GravEnergySource, &
147 "grav/output/accel" / 1, &
148 "grav/pmass" / pmass)
149
150 timedisc => dict( &
151 "method" / ssprk, &
152! "fargo" / 1, &
153 "tol_rel" / 1.0e-3, &
154 "tol_abs" / (/ 0., 1.e-3, 1.e-3, 0. /), &
155 "cfl" / 0.4, &
156 "stoptime" / t, &
157 "dtlimit" / 1.0e-10, &
158 "maxiter" / 2000000000, &
159 "output/xmomentum" / 1, &
160 "output/ymomentum" / 1, &
161 "output/energy" / 1, &
162! "rhstype" / 1, &
163 "output/geometrical_sources" / 1, &
164 "output/external_sources" / 1, &
165 !"output/bflux" / 1, &
166 "output/fluxes" / 1, &
167 "output/rhs" / 1)
168
169
170 datafile => dict( &
171 "fileformat" / vtk, &
172 "filename" / "keplerianvortex", &
173 "count" / onum)
174
175 config => dict( &
176 "physics" / physics, &
177 "fluxes" / fluxes, &
178 "mesh" / mesh, &
179 "boundary" / boundary, &
180 "sources" / sources, &
181 "timedisc" / timedisc, &
182 "datafile" / datafile)
183
184 END SUBROUTINE makeconfig
185
186
187 SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes)
188 IMPLICIT NONE
189 !------------------------------------------------------------------------!
190 CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
191 CLASS(mesh_base), INTENT(IN) :: Mesh
192 CLASS(physics_base), INTENT(INOUT) :: Physics
193 CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
194 !------------------------------------------------------------------------!
195 ! Local variable declaration
196 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
197 :: dvr, dvphi, x, y, phi, r
198 REAL :: h, kappa
199 INTEGER :: dir,j,ig,i,k
200 !------------------------------------------------------------------------!
201 x = mesh%bccart(:,:,:,1)
202 y = mesh%bccart(:,:,:,2)
203
204 IF(gravenergysource.EQ.0) THEN
205 CALL timedisc%Error("keplerianvortex:InitData","GravEnergySource=0 is currently not supported")
206! CALL SetPotential_euler2Dia(Physics,Mesh,-GN * MCENTRAL / Mesh%radius%bcenter)
207! CALL SetPotential_euler2Dia(Physics,Mesh,-GN * MCENTRAL / Mesh%radius%faces)
208 END IF
209
210 SELECT TYPE (pvar => timedisc%pvar)
211 CLASS IS(statevector_euler)
212 SELECT TYPE (cvar => timedisc%cvar)
213 CLASS IS(statevector_euler)
214 cvar%density%data3d(:,:,:) = sigma0
215 cvar%momentum%data4d(:,:,:,1:physics%VDIM) = 0.
216
217 CALL physics%Convert2Primitive(cvar,pvar)
218
219 pvar%pressure%data3d(:,:,:) = 1./(gamma*m**2)
220
221 pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
222 pvar%velocity%data4d(:,:,:,1:physics%VDIM) &
223 + timedisc%GetCentrifugalVelocity(mesh,physics,sim%Fluxes,sim%Sources,(/0.,0.,1./))
224
225 END SELECT
226 END SELECT
227
228 ! boundary conditions
229 ! custom boundary conditions at western boundary if requested
230 SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
231 CLASS IS (boundary_custom)
232 CALL bwest%SetCustomBoundaries(mesh,physics, &
233 (/custom_reflect,custom_reflneg,custom_kepler,custom_reflect/))
234 END SELECT
235 SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
236 CLASS IS (boundary_custom)
237 CALL beast%SetCustomBoundaries(mesh,physics, &
238 (/custom_reflect,custom_reflneg,custom_kepler,custom_reflect/))
239 END SELECT
240
241 h = 1./(2.*m)
242 kappa = -1.
243 phi = mesh%center(:,:,:,2)
244 x = x - 1./sqrt(2.)
245 y = y - 1./sqrt(2.)
246 r = sqrt(x**2 + y**2)
247
248 dvr = kappa * exp(-r**2/h**2) * (-y*cos(phi) + x*sin(phi))
249 dvphi = kappa * exp(-r**2/h**2) * ( y*sin(phi) + x*cos(phi))
250
251 SELECT TYPE (pvar => timedisc%pvar)
252 CLASS IS(statevector_euler)
253 SELECT TYPE (cvar => timedisc%cvar)
254 CLASS IS(statevector_euler)
255 pvar%velocity%data4d(:,:,:,1) = pvar%velocity%data4d(:,:,:,1) + dvr
256 pvar%velocity%data4d(:,:,:,2) = pvar%velocity%data4d(:,:,:,2) + dvphi
257 END SELECT
258 END SELECT
259
260 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
261
262 END SUBROUTINE initdata
263END PROGRAM keplerianvortex
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
makeconfig
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
keplerianvortex
program keplerianvortex
Definition: keplerianvortex.f90:40
fosite_mod
Definition: fosite.f90:41
Generated by   doxygen 1.9.4