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
b8e7f90a
Commit
b8e7f90a
authored
Aug 14, 2016
by
Matthieu Schaller
Browse files
Added the Kelvin-Helmholtz instability test
parent
35a653b7
Changes
5
Hide whitespace changes
Inline
Side-by-side
examples/GreshoVortex/makeIC.py
View file @
b8e7f90a
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 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
...
...
@@ -19,7 +18,6 @@
##############################################################################
import
h5py
import
random
from
numpy
import
*
import
sys
...
...
@@ -51,7 +49,6 @@ 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
)
...
...
examples/KelvinHelmoltz/kelvinHelmholtz.yml
0 → 100644
View file @
b8e7f90a
# 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.5
# 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
:
kelvinHelmholtz
# Common part of the name of output files
time_first
:
0.
# Time of the first output (in internal units)
delta_time
:
0.25
# 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.01
# 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
:
./kelvinHelmholtz.hdf5
# The file to read
examples/KelvinHelmoltz/makeIC.py
0 → 100644
View file @
b8e7f90a
###############################################################################
# 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
from
numpy
import
*
import
sys
# Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
# Parameters
L2
=
128
# Particles along one edge in the low-density region
gamma
=
5.
/
3.
# Gas adiabatic index
P1
=
2.5
# Central region pressure
P2
=
2.5
# Outskirts pressure
v1
=
0.5
# Central region velocity
v2
=
-
0.5
# Outskirts vlocity
rho1
=
2
# Central density
rho2
=
1
# Outskirts density
omega0
=
0.1
sigma
=
0.05
/
sqrt
(
2
)
fileOutputName
=
"kelvinHelmholtz.hdf5"
#---------------------------------------------------
# Start by generating grids of particles at the two densities
numPart2
=
L2
*
L2
L1
=
int
(
sqrt
(
numPart2
/
rho2
*
rho1
))
numPart1
=
L1
*
L1
#print "N2 =", numPart2, "N1 =", numPart1
#print "L2 =", L2, "L1 = ", L1
#print "rho2 =", rho2, "rho1 =", (float(L1*L1)) / (float(L2*L2))
coords1
=
zeros
((
numPart1
,
3
))
coords2
=
zeros
((
numPart2
,
3
))
h1
=
ones
(
numPart1
)
*
1.2348
/
L1
h2
=
ones
(
numPart2
)
*
1.2348
/
L2
m1
=
zeros
(
numPart1
)
m2
=
zeros
(
numPart2
)
u1
=
zeros
(
numPart1
)
u2
=
zeros
(
numPart2
)
vel1
=
zeros
((
numPart1
,
3
))
vel2
=
zeros
((
numPart2
,
3
))
# Particles in the central region
for
i
in
range
(
L1
):
for
j
in
range
(
L1
):
index
=
i
*
L1
+
j
x
=
i
/
float
(
L1
)
+
1.
/
(
2.
*
L1
)
y
=
j
/
float
(
L1
)
+
1.
/
(
2.
*
L1
)
coords1
[
index
,
0
]
=
x
coords1
[
index
,
1
]
=
y
u1
[
index
]
=
P1
/
(
rho1
*
(
gamma
-
1.
))
vel1
[
index
,
0
]
=
v1
# Particles in the outskirts
for
i
in
range
(
L2
):
for
j
in
range
(
L2
):
index
=
i
*
L2
+
j
x
=
i
/
float
(
L2
)
+
1.
/
(
2.
*
L2
)
y
=
j
/
float
(
L2
)
+
1.
/
(
2.
*
L2
)
coords2
[
index
,
0
]
=
x
coords2
[
index
,
1
]
=
y
u2
[
index
]
=
P2
/
(
rho2
*
(
gamma
-
1.
))
vel2
[
index
,
0
]
=
v2
# Now concatenate arrays
where1
=
abs
(
coords1
[:,
1
]
-
0.5
)
<
0.25
where2
=
abs
(
coords2
[:,
1
]
-
0.5
)
>
0.25
coords
=
append
(
coords1
[
where1
,
:],
coords2
[
where2
,
:],
axis
=
0
)
#print L2*(L2/2), L1*(L1/2)
#print shape(coords), shape(coords1[where1,:]), shape(coords2[where2,:])
#print shape(coords), shape(logical_not(coords1[where1,:])), shape(logical_not(coords2[where2,:]))
vel
=
append
(
vel1
[
where1
,
:],
vel2
[
where2
,
:],
axis
=
0
)
h
=
append
(
h1
[
where1
],
h2
[
where2
],
axis
=
0
)
m
=
append
(
m1
[
where1
],
m2
[
where2
],
axis
=
0
)
u
=
append
(
u1
[
where1
],
u2
[
where2
],
axis
=
0
)
numPart
=
size
(
h
)
ids
=
linspace
(
1
,
numPart
,
numPart
)
m
[:]
=
(
0.5
*
rho1
+
0.5
*
rho2
)
/
float
(
numPart
)
# Velocity perturbation
vel
[:,
1
]
=
omega0
*
sin
(
4
*
pi
*
coords
[:,
0
])
*
(
exp
(
-
(
coords
[:,
1
]
-
0.25
)
**
2
/
(
2
*
sigma
**
2
))
+
exp
(
-
(
coords
[:,
1
]
-
0.75
)
**
2
/
(
2
*
sigma
**
2
)))
#File
fileOutput
=
h5py
.
File
(
fileOutputName
,
'w'
)
# Header
grp
=
fileOutput
.
create_group
(
"/Header"
)
grp
.
attrs
[
"BoxSize"
]
=
[
1.
,
1.
,
0.1
]
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
[
"NumFileOutputsPerSnapshot"
]
=
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
=
fileOutput
.
create_group
(
"/RuntimePars"
)
grp
.
attrs
[
"PeriodicBoundariesOn"
]
=
1
#Units
grp
=
fileOutput
.
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.
grp
.
attrs
[
"Unit current in cgs (U_I)"
]
=
1.
grp
.
attrs
[
"Unit temperature in cgs (U_T)"
]
=
1.
#Particle group
grp
=
fileOutput
.
create_group
(
"/PartType0"
)
ds
=
grp
.
create_dataset
(
'Coordinates'
,
(
numPart
,
3
),
'd'
)
ds
[()]
=
coords
ds
=
grp
.
create_dataset
(
'Velocities'
,
(
numPart
,
3
),
'f'
)
ds
[()]
=
vel
ds
=
grp
.
create_dataset
(
'Masses'
,
(
numPart
,
1
),
'f'
)
ds
[()]
=
m
.
reshape
((
numPart
,
1
))
ds
=
grp
.
create_dataset
(
'SmoothingLength'
,
(
numPart
,
1
),
'f'
)
ds
[()]
=
h
.
reshape
((
numPart
,
1
))
ds
=
grp
.
create_dataset
(
'InternalEnergy'
,
(
numPart
,
1
),
'f'
)
ds
[()]
=
u
.
reshape
((
numPart
,
1
))
ds
=
grp
.
create_dataset
(
'ParticleIDs'
,
(
numPart
,
1
),
'L'
)
ds
[()]
=
ids
.
reshape
((
numPart
,
1
))
fileOutput
.
close
()
examples/KelvinHelmoltz/plotSolution.py
0 → 100644
View file @
b8e7f90a
###############################################################################
# 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
P1
=
2.5
# Central region pressure
P2
=
2.5
# Outskirts pressure
v1
=
0.5
# Central region velocity
v2
=
-
0.5
# Outskirts vlocity
rho1
=
2
# Central density
rho2
=
1
# Outskirts density
# ---------------------------------------------------------------
# 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
])
# Read the simulation data
sim
=
h5py
.
File
(
"kelvinHelmholtz_%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"
][:,:]
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
)
scatter
(
pos
[:,
0
],
pos
[:,
1
],
c
=
vel
[:,
0
],
cmap
=
"PuBu"
,
edgecolors
=
'face'
,
s
=
4
,
vmin
=-
1.
,
vmax
=
1.
)
text
(
0.97
,
0.97
,
"${
\\
rm{Velocity~along}}~x$"
,
ha
=
"right"
,
va
=
"top"
,
backgroundcolor
=
"w"
)
xlabel
(
"${
\\
rm{Position}}~x$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Position}}~y$"
,
labelpad
=
0
)
xlim
(
0
,
1
)
ylim
(
0
,
1
)
# Radial density profile --------------------------------
subplot
(
232
)
scatter
(
pos
[:,
0
],
pos
[:,
1
],
c
=
rho
,
cmap
=
"PuBu"
,
edgecolors
=
'face'
,
s
=
4
,
vmin
=
0.8
,
vmax
=
2.2
)
text
(
0.97
,
0.97
,
"${
\\
rm{Density}}$"
,
ha
=
"right"
,
va
=
"top"
,
backgroundcolor
=
"w"
)
xlabel
(
"${
\\
rm{Position}}~x$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Position}}~y$"
,
labelpad
=
0
)
xlim
(
0
,
1
)
ylim
(
0
,
1
)
# Radial pressure profile --------------------------------
subplot
(
233
)
scatter
(
pos
[:,
0
],
pos
[:,
1
],
c
=
P
,
cmap
=
"PuBu"
,
edgecolors
=
'face'
,
s
=
4
,
vmin
=
1
,
vmax
=
4
)
text
(
0.97
,
0.97
,
"${
\\
rm{Pressure}}$"
,
ha
=
"right"
,
va
=
"top"
,
backgroundcolor
=
"w"
)
xlabel
(
"${
\\
rm{Position}}~x$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Position}}~y$"
,
labelpad
=
0
)
xlim
(
0
,
1
)
ylim
(
0
,
1
)
# Internal energy profile --------------------------------
subplot
(
234
)
scatter
(
pos
[:,
0
],
pos
[:,
1
],
c
=
u
,
cmap
=
"PuBu"
,
edgecolors
=
'face'
,
s
=
4
,
vmin
=
1.5
,
vmax
=
5.
)
text
(
0.97
,
0.97
,
"${
\\
rm{Internal~energy}}$"
,
ha
=
"right"
,
va
=
"top"
,
backgroundcolor
=
"w"
)
xlabel
(
"${
\\
rm{Position}}~x$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Position}}~y$"
,
labelpad
=
0
)
xlim
(
0
,
1
)
ylim
(
0
,
1
)
# Radial entropy profile --------------------------------
subplot
(
235
)
scatter
(
pos
[:,
0
],
pos
[:,
1
],
c
=
S
,
cmap
=
"PuBu"
,
edgecolors
=
'face'
,
s
=
4
,
vmin
=
0.5
,
vmax
=
3.
)
text
(
0.97
,
0.97
,
"${
\\
rm{Entropy}}$"
,
ha
=
"right"
,
va
=
"top"
,
backgroundcolor
=
"w"
)
xlabel
(
"${
\\
rm{Position}}~x$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Position}}~y$"
,
labelpad
=
0
)
xlim
(
0
,
1
)
ylim
(
0
,
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
,
"Kelvin-Helmholtz instability at $t=%.2f$"
%
(
time
),
fontsize
=
10
)
text
(
-
0.49
,
0.8
,
"Centre:~~~ $(P,
\\
rho, v) = (%.3f, %.3f, %.3f)$"
%
(
P1
,
rho1
,
v1
),
fontsize
=
10
)
text
(
-
0.49
,
0.7
,
"Outskirts: $(P,
\\
rho, v) = (%.3f, %.3f, %.3f)$"
%
(
P2
,
rho2
,
v2
),
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
(
"KelvinHelmholtz.png"
,
dpi
=
200
)
examples/KelvinHelmoltz/run.sh
0 → 100755
View file @
b8e7f90a
#!/bin/bash
# Generate the initial conditions if they are not present.
if
[
!
-e
kelvinHelmholtz.hdf5
]
then
echo
"Generating initial conditions for the Kelvin-Helmholtz example..."
python makeIC.py
fi
# Run SWIFT
../swift
-s
-t
1 kelvinHelmholtz.yml
# Plot the solution
python plotSolution.py 6
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