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
ce3b8404
Commit
ce3b8404
authored
Feb 07, 2020
by
Matthieu Schaller
Browse files
Merge branch 'gresho_vortex_plotSolution_update' into 'master'
Added improved plotting for the Gresho Vortex See merge request
!1013
parents
3fc96373
b0aee79f
Changes
2
Hide whitespace changes
Inline
Side-by-side
examples/HydroTests/GreshoVortex_2D/plotSolution.py
View file @
ce3b8404
###############################################################################
# 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/>.
#
##############################################################################
# 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)
gas_gamma
=
5.
0
/
3.0
# Gas adiabatic index
rho0
=
1
# Gas density
P0
=
0.
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
*
from
scipy
import
stats
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'
:
(
9.90
,
6.45
),
'figure.subplot.left'
:
0.045
,
'figure.subplot.right'
:
0.99
,
'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'
]})
style
.
use
(
"../../../tools/stylesheets/mnras.mplstyle"
)
snap
=
int
(
sys
.
argv
[
1
])
...
...
@@ -69,21 +49,27 @@ 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
]
solution_P
[
i
]
=
P0
+
5.
0
+
12.5
*
solution_r
[
i
]
**
2
solution_v_phi
[
i
]
=
5.
0
*
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
]
solution_P
[
i
]
=
(
P0
+
9.0
+
12.5
*
solution_r
[
i
]
**
2
-
20.0
*
solution_r
[
i
]
+
4.0
*
log
(
solution_r
[
i
]
/
0.2
)
)
solution_v_phi
[
i
]
=
2.0
-
5.0
*
solution_r
[
i
]
else
:
solution_P
[
i
]
=
P0
+
3.
+
4.
*
log
(
2.
)
solution_v_phi
[
i
]
=
0.
solution_P
[
i
]
=
P0
+
3.
0
+
4.
0
*
log
(
2.
0
)
solution_v_phi
[
i
]
=
0.
0
solution_rho
=
ones
(
N
)
*
rho0
solution_s
=
solution_P
/
solution_rho
**
gas_gamma
solution_u
=
solution_P
/
((
gas_gamma
-
1.
)
*
solution_rho
)
solution_s
=
solution_P
/
solution_rho
**
gas_gamma
solution_u
=
solution_P
/
((
gas_gamma
-
1.
0
)
*
solution_rho
)
# Read the simulation data
sim
=
h5py
.
File
(
"gresho_%04d.hdf5"
%
snap
,
"r"
)
sim
=
h5py
.
File
(
"gresho_%04d.hdf5"
%
snap
,
"r"
)
boxSize
=
sim
[
"/Header"
].
attrs
[
"BoxSize"
][
0
]
time
=
sim
[
"/Header"
].
attrs
[
"Time"
][
0
]
scheme
=
sim
[
"/HydroScheme"
].
attrs
[
"Scheme"
]
...
...
@@ -92,133 +78,153 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta
=
sim
[
"/HydroScheme"
].
attrs
[
"Kernel eta"
]
git
=
sim
[
"Code"
].
attrs
[
"Git Revision"
]
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
)
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
)
rho
=
sim
[
"/PartType0/Densities"
][:]
u
=
sim
[
"/PartType0/InternalEnergies"
][:]
S
=
sim
[
"/PartType0/Entropies"
][:]
P
=
sim
[
"/PartType0/Pressures"
][:]
# Bin te data
r_bin_edge
=
np
.
arange
(
0.
,
1.
,
0.02
)
r_bin
=
0.5
*
(
r_bin_edge
[
1
:]
+
r_bin_edge
[:
-
1
])
rho_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
v_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
P_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
S_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
u_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
rho2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
v2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
P2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
S2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
u2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
rho_sigma_bin
=
np
.
sqrt
(
rho2_bin
-
rho_bin
**
2
)
v_sigma_bin
=
np
.
sqrt
(
v2_bin
-
v_bin
**
2
)
P_sigma_bin
=
np
.
sqrt
(
P2_bin
-
P_bin
**
2
)
S_sigma_bin
=
np
.
sqrt
(
S2_bin
-
S_bin
**
2
)
u_sigma_bin
=
np
.
sqrt
(
u2_bin
-
u_bin
**
2
)
r_bin_edge
=
np
.
arange
(
0.
0
,
1.
0
,
0.02
)
r_bin
=
0.5
*
(
r_bin_edge
[
1
:]
+
r_bin_edge
[:
-
1
])
rho_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
v_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
P_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
S_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
u_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
rho2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
v2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
P2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
S2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
u2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
rho_sigma_bin
=
np
.
sqrt
(
rho2_bin
-
rho_bin
**
2
)
v_sigma_bin
=
np
.
sqrt
(
v2_bin
-
v_bin
**
2
)
P_sigma_bin
=
np
.
sqrt
(
P2_bin
-
P_bin
**
2
)
S_sigma_bin
=
np
.
sqrt
(
S2_bin
-
S_bin
**
2
)
u_sigma_bin
=
np
.
sqrt
(
u2_bin
-
u_bin
**
2
)
# Plot the interesting quantities
figure
()
figure
(
figsize
=
(
7
,
7
/
1.6
))
line_color
=
"C4"
binned_color
=
"C2"
binned_marker_size
=
4
scatter_props
=
dict
(
marker
=
"."
,
ms
=
1
,
markeredgecolor
=
"none"
,
alpha
=
0.5
,
zorder
=-
1
,
rasterized
=
True
,
linestyle
=
"none"
,
)
errorbar_props
=
dict
(
color
=
binned_color
,
ms
=
binned_marker_size
,
fmt
=
"."
,
lw
=
1.2
)
# Azimuthal velocity profile -----------------------------
subplot
(
231
)
plot
(
r
,
v_phi
,
'.'
,
color
=
'r'
,
ms
=
0.5
)
plot
(
solution_r
,
solution_v_phi
,
'
--
'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
v_bin
,
yerr
=
v_sigma_bin
,
fmt
=
'.'
,
ms
=
8.0
,
color
=
'b'
,
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
)
plot
(
r
,
v_phi
,
**
scatter_props
)
plot
(
solution_r
,
solution_v_phi
,
"
--
"
,
color
=
line_color
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
v_bin
,
yerr
=
v_sigma_bin
,
**
errorbar_props
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"
Radius $r$"
)
ylabel
(
"Azimuthal
velocity
$
v_
\\
phi$"
)
xlim
(
0
,
R_max
)
ylim
(
-
0.1
,
1.2
)
# Radial density profile --------------------------------
subplot
(
232
)
plot
(
r
,
rho
,
'.'
,
color
=
'r'
,
ms
=
0.5
)
plot
(
solution_r
,
solution_rho
,
'
--
'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
rho_bin
,
yerr
=
rho_sigma_bin
,
fmt
=
'.'
,
ms
=
8.0
,
color
=
'b'
,
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{
Density
}}~
\\
rho$"
,
labelpad
=
0
)
xlim
(
0
,
R_max
)
ylim
(
rho0
-
0.3
,
rho0
+
0.3
)
#yticks([-0.2, -0.1, 0., 0.1, 0.2])
plot
(
r
,
rho
,
**
scatter_props
)
plot
(
solution_r
,
solution_rho
,
"
--
"
,
color
=
line_color
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
rho_bin
,
yerr
=
rho_sigma_bin
,
**
errorbar_props
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"
Radius $r$"
)
ylabel
(
"Density
$
\\
rho$"
)
xlim
(
0
,
R_max
)
ylim
(
rho0
-
0.3
,
rho0
+
0.3
)
#
yticks([-0.2, -0.1, 0., 0.1, 0.2])
# Radial pressure profile --------------------------------
subplot
(
233
)
plot
(
r
,
P
,
'.'
,
color
=
'r'
,
ms
=
0.5
)
plot
(
solution_r
,
solution_P
,
'
--
'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
P_bin
,
yerr
=
P_sigma_bin
,
fmt
=
'.'
,
ms
=
8.0
,
color
=
'b'
,
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{
Pressure
}}~P$"
,
labelpad
=
0
)
plot
(
r
,
P
,
**
scatter_props
)
plot
(
solution_r
,
solution_P
,
"
--
"
,
color
=
line_color
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
P_bin
,
yerr
=
P_sigma_bin
,
**
errorbar_props
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"
Radius $r$"
)
ylabel
(
"Pressure
$P$"
)
xlim
(
0
,
R_max
)
ylim
(
4.9
+
P0
,
P0
+
6.1
)
# Internal energy profile --------------------------------
subplot
(
234
)
plot
(
r
,
u
,
'.'
,
color
=
'r'
,
ms
=
0.5
)
plot
(
solution_r
,
solution_u
,
'
--
'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
u_bin
,
yerr
=
u_sigma_bin
,
fmt
=
'.'
,
ms
=
8.0
,
color
=
'b'
,
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
)
xlim
(
0
,
R_max
)
plot
(
r
,
u
,
**
scatter_props
)
plot
(
solution_r
,
solution_u
,
"
--
"
,
color
=
line_color
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
u_bin
,
yerr
=
u_sigma_bin
,
**
errorbar_props
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"$
Radius $r$"
)
ylabel
(
"Internal
Energy
$u$"
)
xlim
(
0
,
R_max
)
ylim
(
7.3
,
9.1
)
# Radial entropy profile --------------------------------
subplot
(
235
)
plot
(
r
,
S
,
'.'
,
color
=
'r'
,
ms
=
0.5
)
plot
(
solution_r
,
solution_s
,
'
--
'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
S_bin
,
yerr
=
S_sigma_bin
,
fmt
=
'.'
,
ms
=
8.0
,
color
=
'b'
,
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{Entropy}}~S$"
,
labelpad
=
0
)
plot
(
r
,
S
,
**
scatter_props
)
plot
(
solution_r
,
solution_s
,
"
--
"
,
color
=
line_color
,
alpha
=
0.8
,
lw
=
1.2
)
errorbar
(
r_bin
,
S_bin
,
yerr
=
S_sigma_bin
,
**
errorbar_props
)
plot
([
0.2
,
0.2
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
plot
([
0.4
,
0.4
],
[
-
100
,
100
],
":"
,
color
=
line_color
,
alpha
=
0.4
,
lw
=
1.2
)
xlabel
(
"
Radius $r$"
)
ylabel
(
"
Entropy $S$"
)
xlim
(
0
,
R_max
)
ylim
(
4.9
+
P0
,
P0
+
6.1
)
# Image --------------------------------------------------
#subplot(234)
#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)
# Information -------------------------------------
subplot
(
236
,
frameon
=
False
)
text
(
-
0.49
,
0.9
,
"Gresho-Chan vortex with $
\\
gamma=%.3f$ at $t=%.2f$"
%
(
gas_gamma
,
time
),
fontsize
=
10
)
text
(
-
0.49
,
0.8
,
"Background $
\\
rho_0=%.3f$"
%
rho0
,
fontsize
=
10
)
text
(
-
0.49
,
0.7
,
"Background $P_0=%.3f$"
%
P0
,
fontsize
=
10
)
plot
([
-
0.49
,
0.1
],
[
0.62
,
0.62
],
'k-'
,
lw
=
1
)
text
(
-
0.49
,
0.5
,
"$
\\
textsc{Swift}$ %s"
%
git
,
fontsize
=
10
)
text
(
-
0.49
,
0.4
,
scheme
,
fontsize
=
10
)
text
(
-
0.49
,
0.3
,
kernel
,
fontsize
=
10
)
text
(
-
0.49
,
0.2
,
"$%.2f$ neighbours ($
\\
eta=%.3f$)"
%
(
neighbours
,
eta
),
fontsize
=
10
)
text_fontsize
=
5
text
(
-
0.49
,
0.9
,
"Gresho-Chan vortex (2D) with $
\\
gamma=%.3f$ at $t=%.2f$"
%
(
gas_gamma
,
time
),
fontsize
=
text_fontsize
,
)
text
(
-
0.49
,
0.8
,
"Background $
\\
rho_0=%.3f$"
%
rho0
,
fontsize
=
text_fontsize
)
text
(
-
0.49
,
0.7
,
"Background $P_0=%.3f$"
%
P0
,
fontsize
=
text_fontsize
)
plot
([
-
0.49
,
0.1
],
[
0.62
,
0.62
],
"k-"
,
lw
=
1
)
text
(
-
0.49
,
0.5
,
"SWIFT %s"
%
git
.
decode
(
"utf-8"
),
fontsize
=
text_fontsize
)
text
(
-
0.49
,
0.4
,
scheme
.
decode
(
"utf-8"
),
fontsize
=
text_fontsize
)
text
(
-
0.49
,
0.3
,
kernel
.
decode
(
"utf-8"
),
fontsize
=
text_fontsize
)
text
(
-
0.49
,
0.2
,
"$%.2f$ neighbours ($
\\
eta=%.3f$)"
%
(
neighbours
,
eta
),
fontsize
=
text_fontsize
,
)
xlim
(
-
0.5
,
0.5
)
ylim
(
0
,
1
)
xticks
([])
yticks
([])
savefig
(
"GreshoVortex.png"
,
dpi
=
200
)
tight_layout
()
savefig
(
"GreshoVortex.png"
)
examples/HydroTests/GreshoVortex_3D/plotSolution.py
View file @
ce3b8404
...
...
@@ -22,43 +22,23 @@
# 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)
gas_gamma
=
5.
0
/
3.0
# Gas adiabatic index
rho0
=
1
# Gas density
P0
=
0.
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
*
from
scipy
import
stats
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'
:
(
9.90
,
6.45
),
'figure.subplot.left'
:
0.045
,
'figure.subplot.right'
:
0.99
,
'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'
]})
style
.
use
(
"../../../tools/stylesheets/mnras.mplstyle"
)
snap
=
int
(
sys
.
argv
[
1
])
...
...
@@ -72,21 +52,27 @@ 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
]
solution_P
[
i
]
=
P0
+
5.
0
+
12.5
*
solution_r
[
i
]
**
2
solution_v_phi
[
i
]
=
5.
0
*
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
]
solution_P
[
i
]
=
(
P0
+
9.0
+
12.5
*
solution_r
[
i
]
**
2
-
20.0
*
solution_r
[
i
]
+
4.0
*
log
(
solution_r
[
i
]
/
0.2
)
)
solution_v_phi
[
i
]
=
2.0
-
5.0
*
solution_r
[
i
]
else
:
solution_P
[
i
]
=
P0
+
3.
+
4.
*
log
(
2.
)
solution_v_phi
[
i
]
=
0.
solution_P
[
i
]
=
P0
+
3.
0
+
4.
0
*
log
(
2.
0
)
solution_v_phi
[
i
]
=
0.
0
solution_rho
=
ones
(
N
)
*
rho0
solution_s
=
solution_P
/
solution_rho
**
gas_gamma
solution_u
=
solution_P
/
((
gas_gamma
-
1.
)
*
solution_rho
)
solution_s
=
solution_P
/
solution_rho
**
gas_gamma
solution_u
=
solution_P
/
((
gas_gamma
-
1.
0
)
*
solution_rho
)
# Read the simulation data
sim
=
h5py
.
File
(
"gresho_%04d.hdf5"
%
snap
,
"r"
)
sim
=
h5py
.
File
(
"gresho_%04d.hdf5"
%
snap
,
"r"
)
boxSize
=
sim
[
"/Header"
].
attrs
[
"BoxSize"
][
0
]
time
=
sim
[
"/Header"
].
attrs
[
"Time"
][
0
]
scheme
=
sim
[
"/HydroScheme"
].
attrs
[
"Scheme"
]
...
...
@@ -95,14 +81,14 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta
=
sim
[
"/HydroScheme"
].
attrs
[
"Kernel eta"
]
git
=
sim
[
"Code"
].
attrs
[
"Git Revision"
]
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
)
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
)
rho
=
sim
[
"/PartType0/Densities"
][:]
u
=
sim
[
"/PartType0/InternalEnergies"
][:]
S
=
sim
[
"/PartType0/Entropies"
][:]
...
...
@@ -121,144 +107,183 @@ except:
plot_viscosity
=
False
# Bin te data
r_bin_edge
=
np
.
arange
(
0.
,
1.
,
0.02
)
r_bin
=
0.5
*
(
r_bin_edge
[
1
:]
+
r_bin_edge
[:
-
1
])
rho_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
v_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
P_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
S_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
u_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
rho2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
v2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
P2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
S2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
u2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
**
2
,
statistic
=
'
mean
'
,
bins
=
r_bin_edge
)
rho_sigma_bin
=
np
.
sqrt
(
rho2_bin
-
rho_bin
**
2
)
v_sigma_bin
=
np
.
sqrt
(
v2_bin
-
v_bin
**
2
)
P_sigma_bin
=
np
.
sqrt
(
P2_bin
-
P_bin
**
2
)
S_sigma_bin
=
np
.
sqrt
(
S2_bin
-
S_bin
**
2
)
u_sigma_bin
=
np
.
sqrt
(
u2_bin
-
u_bin
**
2
)
r_bin_edge
=
np
.
arange
(
0.
0
,
1.
0
,
0.02
)
r_bin
=
0.5
*
(
r_bin_edge
[
1
:]
+
r_bin_edge
[:
-
1
])
rho_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
v_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
P_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
S_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
u_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
rho2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
rho
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
v2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
v_phi
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
P2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
P
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
S2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
S
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
u2_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
u
**
2
,
statistic
=
"
mean
"
,
bins
=
r_bin_edge
)
rho_sigma_bin
=
np
.
sqrt
(
rho2_bin
-
rho_bin
**
2
)
v_sigma_bin
=
np
.
sqrt
(
v2_bin
-
v_bin
**
2
)
P_sigma_bin
=
np
.
sqrt
(
P2_bin
-
P_bin
**
2
)
S_sigma_bin
=
np
.
sqrt
(
S2_bin
-
S_bin
**
2
)
u_sigma_bin
=
np
.
sqrt
(
u2_bin
-
u_bin
**
2
)
if
plot_diffusion
:
alpha_diff_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
diffusion
,
statistic
=
'mean'
,
bins
=
r_bin_edge
)
alpha2_diff_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
diffusion
**
2
,
statistic
=
'mean'
,
bins
=
r_bin_edge
)
alpha_diff_sigma_bin
=
np
.
sqrt
(
alpha2_diff_bin
-
alpha_diff_bin
**
2
)
alpha_diff_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
diffusion
,
statistic
=
"mean"
,
bins
=
r_bin_edge
)
alpha2_diff_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
diffusion
**
2
,
statistic
=
"mean"
,
bins
=
r_bin_edge
)
alpha_diff_sigma_bin
=
np
.
sqrt
(
alpha2_diff_bin
-
alpha_diff_bin
**
2
)
if
plot_viscosity
:
alpha_visc_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
viscosity
,
statistic
=
'mean'
,
bins
=
r_bin_edge
)
alpha2_visc_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
viscosity
**
2
,
statistic
=
'mean'
,
bins
=
r_bin_edge
)
alpha_visc_sigma_bin
=
np
.
sqrt
(
alpha2_visc_bin
-
alpha_visc_bin
**
2
)
alpha_visc_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
viscosity
,
statistic
=
"mean"
,
bins
=
r_bin_edge
)
alpha2_visc_bin
,
_
,
_
=
stats
.
binned_statistic
(
r
,
viscosity
**
2
,
statistic
=
"mean"
,
bins
=
r_bin_edge
)
alpha_visc_sigma_bin
=
np
.
sqrt
(
alpha2_visc_bin
-
alpha_visc_bin
**
2
)
# Plot the interesting quantities
figure
()
figure
(
figsize
=
(
7
,
7
/
1.6
))
line_color
=
"C4"
binned_color
=
"C2"
binned_marker_size
=
4
scatter_props
=
dict
(
marker
=
"."
,
ms
=
1
,
markeredgecolor
=
"none"
,
alpha
=
0.1
,
zorder
=-
1
,
rasterized
=
True
,
linestyle
=
"none"
,