RTI.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: RTI.f90 #
5!# #
6!# Copyright (C) 2008-2024 #
7!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
8!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
9!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
27!----------------------------------------------------------------------------!
48!----------------------------------------------------------------------------!
49PROGRAM rti
50 USE fosite_mod
51#include "tap.h"
52 IMPLICIT NONE
53 !--------------------------------------------------------------------------!
54 ! simulation parameters
55 REAL, PARAMETER :: tsim = 10.0 ! simulation time !
56! REAL, PARAMETER :: DYNVIS = 0.0 ! dynamic viscosity constant !
57! REAL, PARAMETER :: BULKVIS = 0.0 ! bulk viscosity constant !
58 REAL, PARAMETER :: dynvis = 1.0e-4
59 REAL, PARAMETER :: bulkvis = -6.67e-5
60 ! initial condition (SI units)
61 REAL, PARAMETER :: rho0 = 2.0 ! density: upper region !
62 REAL, PARAMETER :: rho1 = 1.0 ! density: lower region !
63 REAL, PARAMETER :: yacc = 0.2 ! grav. acceleration !
64 REAL, PARAMETER :: p0 = 1.2 ! pressure at the top !
65 REAL, PARAMETER :: a0 = 0.02 ! amplitude of init. disturbance !
66 ! mesh settings
67 INTEGER, PARAMETER :: xres = 50 ! resolution in x !
68 INTEGER, PARAMETER :: yres = 100 ! resolution in y !
69 INTEGER, PARAMETER :: zres = 1 ! z-resolution
70 REAL, PARAMETER :: width = 1.0 ! width of comp. domain !
71 REAL, PARAMETER :: height = 2.0 ! height of comp. domaina !
72 ! output parameters
73 INTEGER, PARAMETER :: onum = 10 ! number of output data sets !
74 CHARACTER(LEN=256), PARAMETER &
75 :: odir = './' ! output data dir !
76 CHARACTER(LEN=256), PARAMETER &
77 :: ofname = 'RTI' ! output data file name !
78 !--------------------------------------------------------------------------!
79 CLASS(fosite), ALLOCATABLE :: sim
80 !--------------------------------------------------------------------------!
81
82 tap_plan(1)
83
84 ALLOCATE(sim)
85
86 CALL sim%InitFosite()
87 CALL makeconfig(sim, sim%config)
88 CALL sim%Setup()
89 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc%pvar, sim%Timedisc%cvar)
90 CALL sim%Run()
91 CALL sim%Finalize()
92 DEALLOCATE(sim)
93
94 tap_check(.true.,"Simulation finished")
95 tap_done
96
97 CONTAINS
98
100 SUBROUTINE makeconfig(Sim, config)
101 IMPLICIT NONE
102 !------------------------------------------------------------------------!
103 CLASS(fosite) :: Sim
104 TYPE(dict_typ), POINTER :: config
105 !------------------------------------------------------------------------!
106 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
107 sources, timedisc, fluxes, caccel, vis
108 REAL :: zmax
109 !------------------------------------------------------------------------!
110 INTENT(INOUT) :: sim
111 !------------------------------------------------------------------------!
112 IF (zres.GT.1) THEN
113 zmax = 0.5*width
114 ELSE
115 zmax = 0.0
116 END IF
117 ! mesh settings
118 mesh => dict( &
119 "meshtype" / midpoint, & ! use midpoint rule !
120 "geometry" / cartesian, & ! cartesian grid !
121 "inum" / xres, & ! resolution in x-direction !
122 "jnum" / yres, & ! resolution in y-direction !
123 "knum" / zres, & ! resolution in z-direction !
124 "xmin" / (-0.5*width), & ! minimum value in x-dir. !
125 "xmax" / (0.5*width), & ! maximum value in x-dir. !
126 "ymin" / 0., & ! minimum value in y-dir. !
127 "ymax" / height, & ! maximum value in y-dir. !
128 "zmin" / (-zmax), & ! minimum value in z-dir. !
129 "zmax" / zmax) ! maximum value in z-dir. !
130
131 ! physics settings
132 physics => dict( &
133 "problem" / euler, & ! standard hydrodynamics !
134 "gamma" / 1.4) ! ratio of specific heats !
135
136 ! flux calculation and reconstruction method
137 fluxes => dict( &
138 "fluxtype" / kt, & ! use Kurganov-Tadmor flux !
139 "order" / linear, & ! linear reconstruction !
140 "variables" / conservative, & ! vars. for reconstruction !
141 "limiter" / monocent, & ! type of the limiter !
142 "theta" / 1.2) ! optional param. for MC !
143
144 ! boundary conditions
145 boundary => dict( &
146 "western" / reflecting, & ! reflecting boundary cond. !
147 "eastern" / reflecting, & ! reflecting boundary cond. !
148 "southern" / reflecting, & ! reflecting boundary cond. !
149 "northern" / reflecting, & ! reflecting boundary cond. !
150 "bottomer" / reflecting, & ! reflecting boundary cond. !
151 "topper" / reflecting) ! reflecting boundary cond. !
152
153 ! c-acceleration term
154 caccel => dict( &
155 "stype" / c_accel, & ! source 1: acceleration !
156 "yaccel" / (-yacc)) ! constant acc.in y-dir. !
157
158 ! add to sources
159 sources => dict("caccel" / caccel) ! incl. accel. from above !
160
161
162 ! viscosity source term
163 IF (dynvis.GT.tiny(dynvis).OR.bulkvis.GT.tiny(bulkvis)) THEN
164 vis => dict( &
165 "stype" / viscosity, & ! viscosity source !
166 "vismodel" / molecular, & ! visc. model: molecular !
167 "dynconst" / dynvis, & ! const. dynamic viscosity !
168 "bulkconst" / bulkvis) ! const. bulk viscosity !
169 ! add to sources
170 CALL setattr(sources, "vis", vis)
171 END IF
172
173 ! time discretization settings
174 timedisc => dict( &
175 "method" / modified_euler, &! use modified euler !
176 "order" / 3, & ! third order accuracy !
177 "cfl" / 0.4, & ! courant-number !
178 "stoptime" / tsim, & ! simulation stop-time !
179 "dtlimit" / 1.0e-4, & ! smallest allowed timestep !
180 "maxiter" / 100000) ! max. iters before abort !
181
182 ! initialize data input/output
183 datafile => dict( &
184 "fileformat" / vtk, &
185 "filename" / (trim(odir) // trim(ofname)), &
186 "count" / onum)
187
188 ! collect all above dicts in the configuration dict
189 config => dict( &
190 "mesh" / mesh, & ! all mesh-settings !
191 "physics" / physics, & ! all physics settings !
192 "boundary" / boundary, & ! all bounary settings !
193 "fluxes" / fluxes, & ! all fluxes settings !
194 "sources" / sources, & ! all sources !
195 "timedisc" / timedisc, & ! all timedisc settings !
196 "datafile" / datafile) ! all input/output settings !
197 END SUBROUTINE makeconfig
198
199
201 SUBROUTINE initdata(Mesh,Physics,pvar,cvar)
202 IMPLICIT NONE
203 !------------------------------------------------------------------------!
204 CLASS(mesh_base), INTENT(IN) :: Mesh
205 CLASS(physics_base), INTENT(IN) :: Physics
206 CLASS(marray_compound), POINTER, INTENT(INOUT) :: pvar,cvar
207 !------------------------------------------------------------------------!
208 ! Local variable declaration
209 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
210 :: y0
211 !------------------------------------------------------------------------!
212 ! y-coordiante of the contact surface between the two fluids:
213 ! half of the height with a little dip on the axis (x=0,y=0)
214 y0(:,:,:) = 0.5*mesh%ymax - a0*exp(-0.5/(0.4*width)**2 &
215 * (mesh%bcenter(:,:,:,1)**2+mesh%bcenter(:,:,:,3)**2))
216
217 ! initial hydrostatic stratification
218 SELECT TYPE(p => pvar)
219 TYPE IS(statevector_euler)
220 WHERE (mesh%bcenter(:,:,:,2).GT.y0(:,:,:))
221 ! upper fluid
222 p%density%data3d(:,:,:) = rho0
223 p%pressure%data3d(:,:,:) = p0 + yacc * rho0 * (mesh%ymax-mesh%bcenter(:,:,:,2))
224 ELSEWHERE
225 ! lower fluid
226 p%density%data3d(:,:,:) = rho1
227 p%pressure%data3d(:,:,:) = p0 + yacc * (rho1 * (y0(:,:,:)-mesh%bcenter(:,:,:,2)) &
228 + rho0 * (mesh%ymax-y0(:,:,:)))
229 END WHERE
230 ! velocity vanishes everywhere
231 p%velocity%data1d(:) = 0.
232 CLASS DEFAULT
233 CALL physics%Error("RTI::InitData","only non-isothermal HD supported")
234 END SELECT
235
236 CALL physics%Convert2Conservative(pvar,cvar)
237 CALL mesh%Info(" DATA-----> initial condition: " // &
238 –"RayleighTaylor instability")
239 END SUBROUTINE initdata
240END PROGRAM rti
rti
program rti
Definition: RTI.f90:49
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
makeconfig
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
fosite_mod
Definition: fosite.f90:41
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
Generated by   doxygen 1.9.4