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
98b8c242
Commit
98b8c242
authored
Aug 12, 2016
by
Peter W. Draper
Browse files
Merge remote-tracking branch 'origin/master' into tasks_cleanup
Conflicts: src/task.c
parents
35b95f0a
8f24cc3f
Changes
90
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
98b8c242
...
...
@@ -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
...
...
CONTRIBUTING.md
View file @
98b8c242
The SWIFT source code is using a variation of the 'Google' formatting style.
The script 'format.sh' in the root directory applies the clang-format-3.8
tool with our style choices to all the SWIFT C source file. Please apply
the formatting script to the files before submitting a merge request.
\ No newline at end of file
the formatting script to the files before submitting a merge request.
The SWIFT code comes with a series of unit tests that are run automatically
when a push to the master branch occurs. The suite can be run by doing a
`make
check`
in the root directory. Please check that the test suite still
runs with your changes applied before submitting a merge request and add
relevant unit tests probing the correctness of new modules. An example of how
to add a test to the suite can be found by considering the tests/testGreeting
case.
\ No newline at end of file
configure.ac
View file @
98b8c242
...
...
@@ -469,7 +469,7 @@ if test "$enable_warn" != "no"; then
CFLAGS="$CFLAGS -Wall"
;;
intel)
CFLAGS="$CFLAGS -w2"
CFLAGS="$CFLAGS -w2
-Wunused-variable
"
;;
*)
AX_CFLAGS_WARN_ALL
...
...
examples/GreshoVortex/getGlass.sh
0 → 100755
View file @
98b8c242
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
examples/GreshoVortex/gresho.yml
0 → 100644
View file @
98b8c242
# 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-1
# 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 @
98b8c242
...
...
@@ -21,93 +21,83 @@
import
h5py
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
periodic
=
1
# 1 For periodic box
factor
=
3
boxSize
=
[
1.0
,
1.0
,
1.0
/
factor
]
L
=
120
# Number of particles along one axis
gamma
=
5.
/
3.
# Gas adiabatic index
eta
=
1.2349
# 48 ngbs with cubic spline kernel
rho
=
1
# Gas density
rho0
=
1
# Gas density
P0
=
0.
# Constant additional pressure (should have no impact on the dynamics)
fileName
=
"greshoVortex.hdf5"
vol
=
boxSize
[
0
]
*
boxSize
[
1
]
*
boxSize
[
2
]
fileOutputName
=
"greshoVortex.hdf5"
fileGlass
=
"glassPlane_128.hdf5"
#---------------------------------------------------
numPart
=
L
**
3
/
factor
mass
=
boxSize
[
0
]
*
boxSize
[
1
]
*
boxSize
[
2
]
*
rho
/
numPart
#Generate particles
coords
=
zeros
((
numPart
,
3
))
v
=
zeros
((
numPart
,
3
))
m
=
zeros
((
numPart
,
1
))
h
=
zeros
((
numPart
,
1
))
u
=
zeros
((
numPart
,
1
))
ids
=
zeros
((
numPart
,
1
),
dtype
=
'L'
)
partId
=
0
for
i
in
range
(
L
):
for
j
in
range
(
L
):
for
k
in
range
(
L
/
factor
):
index
=
i
*
L
*
L
/
factor
+
j
*
L
/
factor
+
k
x
=
i
*
boxSize
[
0
]
/
L
+
boxSize
[
0
]
/
(
2
*
L
)
y
=
j
*
boxSize
[
0
]
/
L
+
boxSize
[
0
]
/
(
2
*
L
)
z
=
k
*
boxSize
[
0
]
/
L
+
boxSize
[
0
]
/
(
2
*
L
)
r2
=
(
x
-
boxSize
[
0
]
/
2
)
**
2
+
(
y
-
boxSize
[
1
]
/
2
)
**
2
r
=
sqrt
(
r2
)
coords
[
index
,
0
]
=
x
coords
[
index
,
1
]
=
y
coords
[
index
,
2
]
=
z
v_phi
=
0.
if
r
<
0.2
:
v_phi
=
5.
*
r
elif
r
<
0.4
:
v_phi
=
2.
-
5.
*
r
else
:
v_phi
=
0.
v
[
index
,
0
]
=
-
v_phi
*
(
y
-
boxSize
[
0
]
/
2
)
/
r
v
[
index
,
1
]
=
v_phi
*
(
x
-
boxSize
[
0
]
/
2
)
/
r
v
[
index
,
2
]
=
0.
m
[
index
]
=
mass
h
[
index
]
=
eta
*
boxSize
[
0
]
/
L
P
=
P0
if
r
<
0.2
:
P
=
P
+
5.
+
12.5
*
r2
elif
r
<
0.4
:
P
=
P
+
9.
+
12.5
*
r2
-
20.
*
r
+
4.
*
log
(
r
/
0.2
)
else
:
P
=
P
+
3.
+
4.
*
log
(
2.
)
u
[
index
]
=
P
/
((
gamma
-
1.
)
*
rho
)
ids
[
index
]
=
partId
+
1
partId
=
partId
+
1
# 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"
][:]
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
]
r2
=
(
x
-
boxSize
/
2
)
**
2
+
(
y
-
boxSize
/
2
)
**
2
r
=
sqrt
(
r2
)
v_phi
=
0.
if
r
<
0.2
:
v_phi
=
5.
*
r
elif
r
<
0.4
:
v_phi
=
2.
-
5.
*
r
else
:
v_phi
=
0.
v
[
i
,
0
]
=
-
v_phi
*
(
y
-
boxSize
/
2
)
/
r
v
[
i
,
1
]
=
v_phi
*
(
x
-
boxSize
/
2
)
/
r
v
[
i
,
2
]
=
0.
P
=
P0
if
r
<
0.2
:
P
=
P
+
5.
+
12.5
*
r2
elif
r
<
0.4
:
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.
)
*
rho0
)
#File
file
=
h5py
.
File
(
fileName
,
'w'
)
file
Output
=
h5py
.
File
(
file
Output
Name
,
'w'
)
# Header
grp
=
file
.
create_group
(
"/Header"
)
grp
=
file
Output
.
create_group
(
"/Header"
)
grp
.
attrs
[
"BoxSize"
]
=
boxSize
grp
.
attrs
[
"NumPart_Total"
]
=
[
numPart
,
0
,
0
,
0
,
0
,
0
]
grp
.
attrs
[
"NumPart_Total_HighWord"
]
=
[
0
,
0
,
0
,
0
,
0
,
0
]
grp
.
attrs
[
"NumPart_ThisFile"
]
=
[
numPart
,
0
,
0
,
0
,
0
,
0
]
grp
.
attrs
[
"Time"
]
=
0.0
grp
.
attrs
[
"NumFilesPerSnapshot"
]
=
1
grp
.
attrs
[
"NumFile
Output
sPerSnapshot"
]
=
1
grp
.
attrs
[
"MassTable"
]
=
[
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
]
grp
.
attrs
[
"Flag_Entropy_ICs"
]
=
[
0
,
0
,
0
,
0
,
0
,
0
]
#Runtime parameters
grp
=
file
.
create_group
(
"/RuntimePars"
)
grp
.
attrs
[
"PeriodicBoundariesOn"
]
=
periodic
grp
=
file
Output
.
create_group
(
"/RuntimePars"
)
grp
.
attrs
[
"PeriodicBoundariesOn"
]
=
1
#Units
grp
=
file
.
create_group
(
"/Units"
)
grp
=
file
Output
.
create_group
(
"/Units"
)
grp
.
attrs
[
"Unit length in cgs (U_L)"
]
=
1.
grp
.
attrs
[
"Unit mass in cgs (U_M)"
]
=
1.
grp
.
attrs
[
"Unit time in cgs (U_t)"
]
=
1.
...
...
@@ -115,20 +105,20 @@ grp.attrs["Unit current in cgs (U_I)"] = 1.
grp
.
attrs
[
"Unit temperature in cgs (U_T)"
]
=
1.
#Particle group
grp
=
file
.
create_group
(
"/PartType0"
)
grp
=
file
Output
.
create_group
(
"/PartType0"
)
ds
=
grp
.
create_dataset
(
'Coordinates'
,
(
numPart
,
3
),
'd'
)
ds
[()]
=
coords
ds
=
grp
.
create_dataset
(
'Velocities'
,
(
numPart
,
3
),
'f'
)
ds
[()]
=
v
ds
=
grp
.
create_dataset
(
'Masses'
,
(
numPart
,
1
),
'f'
)
ds
[()]
=
m
ds
=
grp
.
create_dataset
(
'Masses'
,
(
numPart
,
1
),
'f'
)
ds
[()]
=
m
.
reshape
((
numPart
,
1
))
ds
=
grp
.
create_dataset
(
'SmoothingLength'
,
(
numPart
,
1
),
'f'
)
ds
[()]
=
h
ds
[()]
=
h
.
reshape
((
numPart
,
1
))
ds
=
grp
.
create_dataset
(
'InternalEnergy'
,
(
numPart
,
1
),
'f'
)
ds
[()]
=
u
ds
[()]
=
u
.
reshape
((
numPart
,
1
))
ds
=
grp
.
create_dataset
(
'ParticleIDs'
,
(
numPart
,
1
),
'L'
)
ds
[()]
=
ids
[:]
ds
[()]
=
ids
.
reshape
((
numPart
,
1
))
file
.
close
()
file
Output
.
close
()
examples/GreshoVortex/plotSolution.py
0 → 100644
View file @
98b8c242
###############################################################################
# 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'
:
(
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'
]})
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"
]
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
)
rho
=
sim
[
"/PartType0/Density"
][:]
u
=
sim
[
"/PartType0/InternalEnergy"
][:]
S
=
sim
[
"/PartType0/Entropy"
][:]
P
=
sim
[
"/PartType0/Pressure"
][:]
# Plot the interesting quantities
figure
()
# 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
)
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
)
# 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
)
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])
# 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
)
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
)
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
)
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
)
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
)
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
)
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
)
xlim
(
-
0.5
,
0.5
)
ylim
(
0
,
1
)
xticks
([])
yticks
([])
savefig
(
"GreshoVortex.png"
,
dpi
=
200
)
examples/GreshoVortex/run.sh
0 → 100755
View file @
98b8c242
#!/bin/bash
# Generate the initial conditions if they are not present.
if
[
!
-e
glassPlane_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
python plotSolution.py 11
examples/GreshoVortex/solution.py
deleted
100644 → 0
View file @
35b95f0a
###############################################################################
# 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
)))
examples/PerturbedBox_2D/makeIC.py
0 → 100644
View file @
98b8c242
###############################################################################
# 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/>.
#
##############################################################################
import
h5py
import
sys
import
random
from
numpy
import
*
# Generates a swift IC file containing a perturbed cartesian distribution of particles
# at a constant density and pressure in a cubic box
# Parameters
periodic
=
1
# 1 For periodic box
boxSize
=
1.
L
=
int
(
sys
.
argv
[
1
])
# Number of particles along one axis
rho
=
1.
# Density
P
=
1.
# Pressure
gamma
=
5.
/
3.
# Gas adiabatic index
pert
=
0.1
# Perturbation scale (in units of the interparticle separation)
fileName
=
"perturbedPlane.hdf5"
#---------------------------------------------------
numPart
=
L
**
2
mass
=
boxSize
**
2
*
rho
/
numPart
internalEnergy
=
P
/
((
gamma
-
1.
)
*
rho
)
#Generate particles
coords
=
zeros
((
numPart
,
3
))
v
=
zeros
((
numPart
,
3
))
m
=
zeros
((
numPart
,
1
))
h
=
zeros
((
numPart
,
1
))
u
=
zeros
((
numPart
,
1
))
ids
=
zeros
((
numPart
,
1
),
dtype
=
'L'
)
for
i
in
range
(
L
):
for
j
in
range
(
L
):
index
=
i
*
L
+
j