Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
f94c65fc
Commit
f94c65fc
authored
Aug 10, 2016
by
Matthieu Schaller
Browse files
Revived the Gresho-Chan vortex test-case. Runs perfectly in 2D.
parent
a8fd4934
Changes
7
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
f94c65fc
...
...
@@ -29,6 +29,7 @@ examples/used_parameters.yml
examples/energy.txt
examples/*/*.xmf
examples/*/*.hdf5
examples/*/*.png
examples/*/*.txt
examples/*/used_parameters.yml
examples/*/*/*.xmf
...
...
examples/GreshoVortex/getGlass.sh
0 → 100755
View file @
f94c65fc
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
examples/GreshoVortex/gresho.yml
0 → 100644
View file @
f94c65fc
# Define the system of units to use internally.
InternalUnitSystem
:
UnitMass_in_cgs
:
1
# Grams
UnitLength_in_cgs
:
1
# Centimeters
UnitVelocity_in_cgs
:
1
# Centimeters per second
UnitCurrent_in_cgs
:
1
# Amperes
UnitTemp_in_cgs
:
1
# Kelvin
# Parameters governing the time integration
TimeIntegration
:
time_begin
:
0.
# The starting time of the simulation (in internal units).
time_end
:
1.
# The end time of the simulation (in internal units).
dt_min
:
1e-6
# The minimal time-step size of the simulation (in internal units).
dt_max
:
1e-2
# The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots
:
basename
:
gresho
# Common part of the name of output files
time_first
:
0.
# Time of the first output (in internal units)
delta_time
:
1e-2
# Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics
:
delta_time
:
1e-2
# Time between statistics output
# Parameters for the hydrodynamics scheme
SPH
:
resolution_eta
:
1.2348
# Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
delta_neighbours
:
0.1
# The tolerance for the targetted number of neighbours.
max_smoothing_length
:
0.1
# Maximal smoothing length allowed (in internal units).
CFL_condition
:
0.1
# Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions
:
file_name
:
./greshoVortex.hdf5
# The file to read
examples/GreshoVortex/makeIC.py
View file @
f94c65fc
...
...
@@ -23,32 +23,32 @@ import random
from
numpy
import
*
import
sys
# Generates a swift IC file for the Gresho
V
ortex in a periodic box
# Generates a swift IC file for the Gresho
-Chan v
ortex in a periodic box
# Parameters
inputFile
=
"glass_256.hdf5"
periodic
=
1
# 1 For periodic box
gamma
=
5.
/
3.
# Gas adiabatic index
rho
=
1
# Gas density
rho
0
=
1
# Gas density
P0
=
0.
# Constant additional pressure (should have no impact on the dynamics)
fileOutputName
=
"greshoVortex.hdf5"
fileGlass
=
"glassPlane_128.hdf5"
#---------------------------------------------------
fileInput
=
h5py
.
File
(
inputFile
,
"r"
)
header
=
fileInput
[
"/Header"
]
boxSize
=
header
.
attrs
[
"BoxSize"
][
0
]
# Get position and smoothing lengths from the glass
fileInput
=
h5py
.
File
(
fileGlass
,
'r'
)
coords
=
fileInput
[
"/PartType0/Coordinates"
][:,:]
h
=
fileInput
[
"/PartType0/SmoothingLength"
][:]
ids
=
fileInput
[
"/PartType0/ParticleIDs"
][:]
numPart
=
size
(
ids
)
m
=
ones
(
shape
(
ids
))
*
boxSize
**
2
/
numPart
u
=
zeros
(
shape
(
m
))
v
=
zeros
(
shape
(
coords
))
boxSize
=
fileInput
[
"/Header"
].
attrs
[
"BoxSize"
][
0
]
numPart
=
size
(
h
)
fileInput
.
close
()
# Now generate the rest
m
=
ones
(
numPart
)
*
rho0
*
boxSize
**
2
/
numPart
u
=
zeros
(
numPart
)
v
=
zeros
((
numPart
,
3
))
for
i
in
range
(
numPart
):
x
=
coords
[
i
,
0
]
y
=
coords
[
i
,
1
]
z
=
coords
[
i
,
2
]
...
...
@@ -74,7 +74,7 @@ for i in range(numPart):
P
=
P
+
9.
+
12.5
*
r2
-
20.
*
r
+
4.
*
log
(
r
/
0.2
)
else
:
P
=
P
+
3.
+
4.
*
log
(
2.
)
u
[
i
]
=
P
/
((
gamma
-
1.
)
*
rho
)
u
[
i
]
=
P
/
((
gamma
-
1.
)
*
rho
0
)
...
...
@@ -94,7 +94,7 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
#Runtime parameters
grp
=
fileOutput
.
create_group
(
"/RuntimePars"
)
grp
.
attrs
[
"PeriodicBoundariesOn"
]
=
periodic
grp
.
attrs
[
"PeriodicBoundariesOn"
]
=
1
#Units
grp
=
fileOutput
.
create_group
(
"/Units"
)
...
...
examples/GreshoVortex/plotSolution.py
0 → 100644
View file @
f94c65fc
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
# Computes the analytical solution of the Gresho-Chan vortex and plots the SPH answer
# Parameters
gas_gamma
=
5.
/
3.
# Gas adiabatic index
rho0
=
1
# Gas density
P0
=
0.
# Constant additional pressure (should have no impact on the dynamics)
# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------
import
matplotlib
matplotlib
.
use
(
"Agg"
)
from
pylab
import
*
import
h5py
# Plot parameters
params
=
{
'axes.labelsize'
:
10
,
'axes.titlesize'
:
10
,
'font.size'
:
12
,
'legend.fontsize'
:
12
,
'xtick.labelsize'
:
10
,
'ytick.labelsize'
:
10
,
'text.usetex'
:
True
,
'figure.figsize'
:
(
6.45
,
6.45
),
'figure.subplot.left'
:
0.07
,
'figure.subplot.right'
:
0.985
,
'figure.subplot.bottom'
:
0.05
,
'figure.subplot.top'
:
0.99
,
'figure.subplot.wspace'
:
0.15
,
'figure.subplot.hspace'
:
0.12
,
'lines.markersize'
:
6
,
'lines.linewidth'
:
3.
,
'text.latex.unicode'
:
True
}
rcParams
.
update
(
params
)
rc
(
'font'
,
**
{
'family'
:
'sans-serif'
,
'sans-serif'
:[
'Times'
]})
snap
=
int
(
sys
.
argv
[
1
])
# Generate the analytic solution at this time
N
=
200
R_max
=
0.8
solution_r
=
arange
(
0
,
R_max
,
R_max
/
N
)
solution_P
=
zeros
(
N
)
solution_v_phi
=
zeros
(
N
)
solution_v_r
=
zeros
(
N
)
for
i
in
range
(
N
):
if
solution_r
[
i
]
<
0.2
:
solution_P
[
i
]
=
P0
+
5.
+
12.5
*
solution_r
[
i
]
**
2
solution_v_phi
[
i
]
=
5.
*
solution_r
[
i
]
elif
solution_r
[
i
]
<
0.4
:
solution_P
[
i
]
=
P0
+
9.
+
12.5
*
solution_r
[
i
]
**
2
-
20.
*
solution_r
[
i
]
+
4.
*
log
(
solution_r
[
i
]
/
0.2
)
solution_v_phi
[
i
]
=
2.
-
5.
*
solution_r
[
i
]
else
:
solution_P
[
i
]
=
P0
+
3.
+
4.
*
log
(
2.
)
solution_v_phi
[
i
]
=
0.
solution_rho
=
ones
(
N
)
*
rho0
solution_s
=
solution_P
/
solution_rho
**
gas_gamma
solution_u
=
solution_P
/
((
gas_gamma
-
1.
)
*
solution_rho
)
# Read the simulation data
sim
=
h5py
.
File
(
"gresho_%03d.hdf5"
%
snap
,
"r"
)
boxSize
=
sim
[
"/Header"
].
attrs
[
"BoxSize"
][
0
]
time
=
sim
[
"/Header"
].
attrs
[
"Time"
][
0
]
scheme
=
sim
[
"/HydroScheme"
].
attrs
[
"Scheme"
]
kernel
=
sim
[
"/HydroScheme"
].
attrs
[
"Kernel function"
]
pos
=
sim
[
"/PartType0/Coordinates"
][:,:]
x
=
pos
[:,
0
]
-
boxSize
/
2
y
=
pos
[:,
1
]
-
boxSize
/
2
vel
=
sim
[
"/PartType0/Velocities"
][:,:]
r
=
sqrt
(
x
**
2
+
y
**
2
)
v_r
=
(
x
*
vel
[:,
0
]
+
y
*
vel
[:,
1
])
/
r
v_phi
=
(
-
y
*
vel
[:,
0
]
+
x
*
vel
[:,
1
])
/
r
v_norm
=
sqrt
(
vel
[:,
0
]
**
2
+
vel
[:,
1
]
**
2
)
u
=
sim
[
"/PartType0/InternalEnergy"
][:]
# Plot the interesting quantities
figure
()
# Internal energy profile --------------------------------
subplot
(
221
)
plot
(
r
,
u
,
'x'
,
color
=
'r'
)
plot
(
solution_r
,
solution_u
,
'--'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
':'
,
color
=
'k'
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
':'
,
color
=
'k'
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"${
\\
rm{Radius}}~r$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Internal~Energy}}~u$"
,
labelpad
=
0
)
text
(
0.02
,
8.97
,
"%s"
%
scheme
,
fontsize
=
9
,
backgroundcolor
=
"w"
)
text
(
0.02
,
8.82
,
"%s"
%
kernel
,
fontsize
=
9
,
backgroundcolor
=
"w"
)
xlim
(
0
,
R_max
)
ylim
(
7.3
,
9.1
)
# Azimuthal velocity profile -----------------------------
subplot
(
222
)
plot
(
r
,
v_phi
,
'x'
,
color
=
'r'
)
plot
(
solution_r
,
solution_v_phi
,
'--'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
':'
,
color
=
'k'
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
':'
,
color
=
'k'
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"${
\\
rm{Radius}}~r$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Azimuthal~velocity}}~v_
\\
phi$"
,
labelpad
=
0
)
xlim
(
0
,
R_max
)
ylim
(
-
0.1
,
1.2
)
text
(
0.75
,
1
,
"$t=%.3f$"
%
(
time
),
ha
=
"right"
,
va
=
"bottom"
)
# Radial velocity profile --------------------------------
subplot
(
223
)
plot
(
r
,
v_r
,
'x'
,
color
=
'r'
,
label
=
"$
\\
texttt{SWIFT}$"
)
plot
(
solution_r
,
solution_v_r
,
'--'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
,
label
=
"$
\\
rm{Solution}$"
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
':'
,
color
=
'k'
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
':'
,
color
=
'k'
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"${
\\
rm{Radius}}~r$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Radial~velocity}}~v_r$"
,
labelpad
=-
5
)
legend
(
loc
=
"upper right"
,
fontsize
=
10
,
frameon
=
False
,
numpoints
=
1
,
handletextpad
=
0.1
,
handlelength
=
1.2
)
xlim
(
0
,
R_max
)
ylim
(
-
0.2
,
0.2
)
yticks
([
-
0.2
,
-
0.1
,
0.
,
0.1
,
0.2
])
# Image --------------------------------------------------
subplot
(
224
)
scatter
(
pos
[:,
0
],
pos
[:,
1
],
c
=
v_norm
,
cmap
=
"PuBu"
,
edgecolors
=
'face'
,
s
=
4
,
vmin
=
0
,
vmax
=
1
)
text
(
0.95
,
0.95
,
"$|v|$"
,
ha
=
"right"
,
va
=
"top"
)
xlim
(
0
,
1
)
ylim
(
0
,
1
)
xlabel
(
"$x$"
,
labelpad
=
0
)
ylabel
(
"$y$"
,
labelpad
=
0
)
savefig
(
"greshoVortex_%03d.png"
%
snap
)
examples/GreshoVortex/run.sh
0 → 100755
View file @
f94c65fc
#!/bin/bash
# Generate the initial conditions if they are not present.
if
[
!
-e
glass_128.hdf5
]
then
echo
"Fetching initial glass file for the Gresho-Chan vortex example..."
./getGlass.sh
fi
if
[
!
-e
greshoVortex.hdf5
]
then
echo
"Generating initial conditions for the Gresho-Chan vortex example..."
python makeIC.py
fi
# Run SWIFT
../swift
-s
-t
1 gresho.yml
# Plot the solution
for
i
in
{
0..100
}
do
python plotSolution.py
$i
done
# Make a movie
mencoder mf://
*
.png
-mf
w
=
645:h
=
645:fps
=
12:type
=
png
-ovc
lavc
-lavcopts
vcodec
=
mpeg4:mbd
=
2:trell
-oac
copy
-o
gresho.avi
examples/GreshoVortex/solution.py
deleted
100644 → 0
View file @
a8fd4934
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import
random
from
numpy
import
*
# Computes the analytical solution of the Gresho-Chan vortex
# The script works for a given initial box and background pressure and computes the solution for any time t (The solution is constant over time).
# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
# entropic function on N radial points between r=0 and r=R_max.
# Parameters
rho0
=
1.
# Background Density
P0
=
0.
# Background Pressure
gamma
=
5.
/
3.
# Gas polytropic index
N
=
1000
# Number of radial points
R_max
=
1.
# Maximal radius
# ---------------------------------------------------------------
# Don't touch anything after this.
# ---------------------------------------------------------------
r
=
arange
(
0
,
R_max
,
R_max
/
N
)
rho
=
ones
(
N
)
P
=
zeros
(
N
)
v
=
zeros
(
N
)
u
=
zeros
(
N
)
s
=
zeros
(
N
)
for
i
in
range
(
N
):
if
r
[
i
]
<
0.2
:
P
[
i
]
=
P0
+
5.
+
12.5
*
r
[
i
]
**
2
v
[
i
]
=
5.
*
r
[
i
]
elif
r
[
i
]
<
0.4
:
P
[
i
]
=
P0
+
9.
+
12.5
*
r
[
i
]
**
2
-
20.
*
r
[
i
]
+
4.
*
log
(
r
[
i
]
/
0.2
)
v
[
i
]
=
2.
-
5.
*
r
[
i
]
else
:
P
[
i
]
=
P0
+
3.
+
4.
*
log
(
2.
)
v
[
i
]
=
0.
rho
[
i
]
=
rho0
s
[
i
]
=
P
[
i
]
/
rho
[
i
]
**
gamma
u
[
i
]
=
P
[
i
]
/
((
gamma
-
1.
)
*
rho
[
i
])
savetxt
(
"rho.dat"
,
column_stack
((
r
,
rho
)))
savetxt
(
"P.dat"
,
column_stack
((
r
,
P
)))
savetxt
(
"v.dat"
,
column_stack
((
r
,
v
)))
savetxt
(
"u.dat"
,
column_stack
((
r
,
u
)))
savetxt
(
"s.dat"
,
column_stack
((
r
,
s
)))
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment