Commit 1957693f authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge remote-tracking branch 'origin/master' into back_of_the_queue

parents c7479fe8 3288efec
......@@ -407,6 +407,24 @@ if test "$enable_memuse_reports" = "yes"; then
AC_DEFINE([SWIFT_MEMUSE_REPORTS],1,[Enable memory usage reports])
fi
# The system memory report depends on the presence of /proc/self/statm, i.e.
# a linux OS, check for that.
if test -f "/proc/self/statm"; then
AC_DEFINE([SWIFT_MEMUSE_STATM],1,[Have /proc/self/statm capability])
fi
# Check if we want to make the dumper thread active.
AC_ARG_ENABLE([dumper],
[AS_HELP_STRING([--enable-dumper],
[Dump active tasks and memory use (if configured)@<:@yes/no@:>@]
)],
[enable_dumper="$enableval"],
[enable_dumper="no"]
)
if test "$enable_dumper" = "yes"; then
AC_DEFINE([SWIFT_DUMPER_THREAD],1,[Enable dumper thread])
fi
# See if we want mpi reporting.
AC_ARG_ENABLE([mpiuse-reports],
[AS_HELP_STRING([--enable-mpiuse-reports],
......@@ -1208,6 +1226,17 @@ fi
AC_SUBST([VELOCIRAPTOR_LIBS])
AM_CONDITIONAL([HAVEVELOCIRAPTOR],[test -n "$VELOCIRAPTOR_LIBS"])
# Now that we found VELOCIraptor, let's check how it was compiled.
if test "$have_velociraptor" == "yes"; then
AC_CHECK_LIB(
[velociraptor],
[VR_NOMASS],
[AC_DEFINE([HAVE_VELOCIRAPTOR_WITH_NOMASS],1,[The VELOCIraptor library has been compiled with the NOMASS option. Only useful if running a uniform box.])],
[AC_MSG_RESULT(VELOCIraptor not compiled to so as to *not* store masses per particle.)],
[$VELOCIRAPTOR_LIBS $HDF5_LDFLAGS $HDF5_LIBS $GSL_LIBS]
)
fi
# Check for dummy VELOCIraptor.
AC_ARG_ENABLE([dummy-velociraptor],
[AS_HELP_STRING([--enable-dummy-velociraptor],
......@@ -1899,7 +1928,7 @@ esac
# Feedback model
AC_ARG_WITH([feedback],
[AS_HELP_STRING([--with-feedback=<model>],
[Feedback model to use @<:@none, EAGLE, debug default: none@:>@]
[Feedback model to use @<:@none, EAGLE, GEAR default: none@:>@]
)],
[with_feedback="$withval"],
[with_feedback="none"]
......@@ -1932,7 +1961,7 @@ esac
# Black hole model.
AC_ARG_WITH([black-holes],
[AS_HELP_STRING([--with-black-holes=<model>],
[Black holes model to use @<:@none, default: none@:>@]
[Black holes model to use @<:@none, EAGLE default: none@:>@]
)],
[with_black_holes="$withval"],
[with_black_holes="none"]
......@@ -1961,14 +1990,14 @@ esac
# Task order
AC_ARG_WITH([task-order],
[AS_HELP_STRING([--with-task-order=<model>],
[Task order to use @<:@none, EAGLE, GEAR default: none@:>@]
[Task order to use @<:@default, EAGLE, GEAR default: default@:>@]
)],
[with_task_order="$withval"],
[with_task_order="none"]
[with_task_order="default"]
)
if test "$with_subgrid" != "none"; then
if test "$with_task_order" != "none"; then
if test "$with_task_order" != "default"; then
AC_MSG_ERROR([Cannot provide with-subgrid and with-task-order together])
else
with_task_order="$with_subgrid_task_order"
......@@ -1979,8 +2008,8 @@ case "$with_task_order" in
EAGLE)
AC_DEFINE([TASK_ORDER_EAGLE], [1], [EAGLE task order])
;;
none)
AC_DEFINE([TASK_ORDER_NONE], [1], [Default task order])
default)
AC_DEFINE([TASK_ORDER_NONE], [1], [Default (i.e. EAGLE/OWLS) task order])
;;
GEAR)
AC_DEFINE([TASK_ORDER_GEAR], [1], [GEAR task order])
......
......@@ -7,18 +7,24 @@
Equations of State
==================
Currently (if the documentation was well updated), we have two different gas
equations of state (EoS) implemented: ideal and isothermal; as well as a variety
of EoS for "planetary" materials.
The EoS describe the relations between our main thermodynamical variables:
the internal energy (\\(u\\)), the density (\\(\\rho\\)), the entropy (\\(A\\))
and the pressure (\\(P\\)).
Currently, SWIFT offers two different gas equations of state (EoS)
implemented: ``ideal`` and ``isothermal``; as well as a variety of EoS for
"planetary" materials. The EoS describe the relations between our
main thermodynamical variables: the internal energy per unit mass
(\\(u\\)), the mass density (\\(\\rho\\)), the entropy (\\(A\\)) and
the pressure (\\(P\\)).
Gas EoS
-------
In the following section, the variables not yet defined are: \\(\\gamma\\) for
the adiabatic index and \\( c_s \\) for the speed of sound.
We write the adiabatic index as \\(\\gamma \\) and \\( c_s \\) denotes
the speed of sound. The adiabatic index can be changed at configure
time by choosing one of the allowed values of the option
``--with-adiabatic-index``. The default value is \\(\\gamma = 5/3 \\).
The tables below give the expression for the thermodynamic quantities
on each row entry as a function of the gas density and the
thermodynamical quantity given in the header of each column.
.. csv-table:: Ideal Gas
:header: "Variable", "A", "u", "P"
......@@ -38,7 +44,11 @@ the adiabatic index and \\( c_s \\) for the speed of sound.
"P", "", "\\(\\left( \\gamma - 1\\right) u \\rho \\)", ""
"\\( c_s\\)", "", "\\(\\sqrt{ u \\gamma \\left( \\gamma - 1 \\right) } \\)", ""
Note that when running with an isothermal equation of state, the value
of the tracked thermodynamic variable (e.g. the entropy in a
density-entropy scheme or the internal enegy in a density-energy SPH
formulation) written to the snapshots is meaningless. The pressure,
however, is always correct in all scheme.
Planetary EoS
-------------
......
......@@ -986,26 +986,34 @@ should contain. The type of partitioning attempted is controlled by the::
DomainDecomposition:
initial_type:
parameter. Which can have the values *memory*, *region*, *grid* or
parameter. Which can have the values *memory*, *edgememory*, *region*, *grid* or
*vectorized*:
* *edgememory*
This is the default if METIS or ParMETIS is available. It performs a
partition based on the memory use of all the particles in each cell.
The total memory per cell is used to weight the cell vertex and all the
associated edges. This attempts to equalize the memory used by all the
ranks but with some consideration given to the need to not cut dense
regions (by also minimizing the edge cut). How successful this
attempt is depends on the granularity of cells and particles and the
number of ranks, clearly if most of the particles are in one cell, or a
small region of the volume, balance is impossible or difficult. Having
more top-level cells makes it easier to calculate a good distribution
(but this comes at the cost of greater overheads).
* *memory*
This is the default if METIS or ParMETIS is available. It performs a
partition based on the memory use of all the particles in each cell,
attempting to equalize the memory used by all the ranks.
How successful this attempt is depends on the granularity of cells and particles
and the number of ranks, clearly if most of the particles are in one cell,
or a small region of the volume, balance is impossible or
difficult. Having more top-level cells makes it easier to calculate a
good distribution (but this comes at the cost of greater overheads).
This is like *edgememory*, but doesn't include any edge weights, it should
balance the particle memory use per rank more exactly (but note effects
like the numbers of cells per rank will also have an effect, as that
changes the need for foreign cells).
* *region*
The one other METIS/ParMETIS option is "region". This attempts to assign equal
numbers of cells to each rank, with the surface area of the regions minimised
(so we get blobs, rather than rectangular volumes of cells).
numbers of cells to each rank, with the surface area of the regions minimised.
If ParMETIS and METIS are not available two other options are possible, but
will give a poorer partition:
......
......@@ -62,11 +62,9 @@ def get_analytic_pressure(x):
def get_data(name):
file = h5py.File(name, "r")
coords = np.array(file["/PartType0/Coordinates"])
rho = np.array(file["/PartType0/Density"])
u = np.array(file["/PartType0/InternalEnergy"])
rho = np.array(file["/PartType0/Densities"])
v = np.array(file["/PartType0/Velocities"])
P = (gamma - 1.) * rho * u
P = np.array(file["/PartType0/Pressures"])
vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )
......@@ -79,7 +77,7 @@ for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
if num < start or num > stop:
continue
print "processing", f, "..."
print("processing", f, "...")
xrange = np.linspace(0., 2. * x_disc, 1000)
time, x, rho, P, v = get_data(f)
......
###############################################################################
# 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")
......@@ -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