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
7fb703af
Commit
7fb703af
authored
Mar 14, 2017
by
James Willis
Browse files
Merge branch 'master' into dopair-vectorisation
parents
d896e167
3b974424
Changes
63
Hide whitespace changes
Inline
Side-by-side
INSTALL.swift
View file @
7fb703af
...
...
@@ -143,6 +143,8 @@ before you can build it.
-
DOXYGEN
:
the
doxygen
library
is
required
to
create
the
SWIFT
API
documentation
.
-
python
:
Examples
and
solution
script
use
python
and
rely
on
the
numpy
library
version
1.8
.
2
or
higher
.
SWIFT
Coding
style
...
...
configure.ac
View file @
7fb703af
...
...
@@ -200,6 +200,20 @@ if test "$enable_debugging_checks" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
fi
# Check if gravity force checks are on for some particles.
AC_ARG_ENABLE([gravity-force-checks],
[AS_HELP_STRING([--enable-gravity-force-checks],
[Activate expensive brute-force gravity checks for a fraction 1/N of all particles @<:@N@:>@]
)],
[gravity_force_checks="$enableval"],
[gravity_force_checks="no"]
)
if test "$gravity_force_checks" == "yes"; then
AC_MSG_ERROR(Need to specify the fraction of particles to check when using --enable-gravity-force-checks!)
elif test "$gravity_force_checks" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
fi
# Define HAVE_POSIX_MEMALIGN if it works.
AX_FUNC_POSIX_MEMALIGN
...
...
@@ -637,13 +651,13 @@ AC_ARG_WITH([hydro-dimension],
)
case "$with_dimension" in
1)
AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D
analysis
])
AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D
solver
])
;;
2)
AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D
analysis
])
AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D
solver
])
;;
3)
AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D
analysis
])
AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D
solver
])
;;
*)
AC_MSG_ERROR([Dimensionality must be 1, 2 or 3])
...
...
@@ -751,7 +765,7 @@ esac
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
[external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch default: none@:>@]
[external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch
, sine-wave
default: none@:>@]
)],
[with_potential="$withval"],
[with_potential="none"]
...
...
@@ -769,6 +783,9 @@ case "$with_potential" in
disc-patch)
AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
;;
sine-wave)
AC_DEFINE([EXTERNAL_POTENTIAL_SINE_WAVE], [1], [Sine wave external potential in 1D])
;;
*)
AC_MSG_ERROR([Unknown external potential: $with_potential])
;;
...
...
@@ -821,8 +838,10 @@ AC_MSG_RESULT([
Riemann solver : $with_riemann
Cooling function : $with_cooling
External potential : $with_potential
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks
])
# Make sure the latest git revision string gets included
...
...
examples/BigCosmoVolume/makeIC.py
View file @
7fb703af
...
...
@@ -107,7 +107,7 @@ if n_copy > 1:
for
i
in
range
(
n_copy
):
for
j
in
range
(
n_copy
):
for
k
in
range
(
n_copy
):
coords
=
np
.
append
(
coords
,
coords_tile
+
np
.
array
([
i
*
boxSize
,
j
*
boxSize
,
k
*
boxSize
]),
axis
=
0
)
coords
=
np
.
append
(
coords
,
coords_tile
+
np
.
array
([
i
*
boxSize
[
0
]
,
j
*
boxSize
[
1
]
,
k
*
boxSize
[
2
]
]),
axis
=
0
)
v
=
np
.
append
(
v
,
v_tile
,
axis
=
0
)
m
=
np
.
append
(
m
,
m_tile
)
h
=
np
.
append
(
h
,
h_tile
)
...
...
examples/CoolingHalo/makeIC.py
View file @
7fb703af
...
...
@@ -227,7 +227,7 @@ ds[()] = u
u
=
np
.
zeros
(
1
)
# Particle IDs
ids
=
1
+
np
.
linspace
(
0
,
N
,
N
,
endpoint
=
False
,
dtype
=
'L'
)
ids
=
1
+
np
.
linspace
(
0
,
N
,
N
,
endpoint
=
False
)
ds
=
grp
.
create_dataset
(
'ParticleIDs'
,
(
N
,
),
'L'
)
ds
[()]
=
ids
...
...
examples/CoolingHaloWithSpin/makeIC.py
View file @
7fb703af
...
...
@@ -233,7 +233,7 @@ ds[()] = u
u
=
np
.
zeros
(
1
)
# Particle IDs
ids
=
1
+
np
.
linspace
(
0
,
N
,
N
,
endpoint
=
False
,
dtype
=
'L'
)
ids
=
1
+
np
.
linspace
(
0
,
N
,
N
,
endpoint
=
False
)
ds
=
grp
.
create_dataset
(
'ParticleIDs'
,
(
N
,
),
'L'
)
ds
[()]
=
ids
...
...
examples/DiscPatch/GravityOnly/makeIC.py
View file @
7fb703af
...
...
@@ -150,7 +150,7 @@ ds[()] = m
m
=
numpy
.
zeros
(
1
)
ids
=
1
+
numpy
.
linspace
(
0
,
numPart
,
numPart
,
endpoint
=
False
,
dtype
=
'L'
)
ids
=
1
+
numpy
.
linspace
(
0
,
numPart
,
numPart
,
endpoint
=
False
)
ds
=
grp1
.
create_dataset
(
'ParticleIDs'
,
(
numPart
,
),
'L'
)
ds
[()]
=
ids
...
...
examples/DiscPatch/HydroStatic/makeIC.py
View file @
7fb703af
...
...
@@ -205,7 +205,7 @@ if (entropy_flag == 1):
else
:
ds
[()]
=
u
ids
=
1
+
numpy
.
linspace
(
0
,
numGas
,
numGas
,
endpoint
=
False
,
dtype
=
'L'
)
ids
=
1
+
numpy
.
linspace
(
0
,
numGas
,
numGas
,
endpoint
=
False
)
ds
=
grp0
.
create_dataset
(
'ParticleIDs'
,
(
numGas
,
),
'L'
)
ds
[()]
=
ids
...
...
examples/HydrostaticHalo/density_profile.py
View file @
7fb703af
...
...
@@ -28,7 +28,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs
=
float
(
params
.
attrs
[
"InternalUnitSystem:UnitLength_in_cgs"
])
unit_velocity_cgs
=
float
(
params
.
attrs
[
"InternalUnitSystem:UnitVelocity_in_cgs"
])
unit_time_cgs
=
unit_length_cgs
/
unit_velocity_cgs
v_c
=
float
(
params
.
attrs
[
"
Softened
IsothermalPotential:vrot"
])
v_c
=
float
(
params
.
attrs
[
"IsothermalPotential:vrot"
])
v_c_cgs
=
v_c
*
unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
...
...
examples/HydrostaticHalo/internal_energy_profile.py
View file @
7fb703af
...
...
@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs
=
float
(
params
.
attrs
[
"InternalUnitSystem:UnitLength_in_cgs"
])
unit_velocity_cgs
=
float
(
params
.
attrs
[
"InternalUnitSystem:UnitVelocity_in_cgs"
])
unit_time_cgs
=
unit_length_cgs
/
unit_velocity_cgs
v_c
=
float
(
params
.
attrs
[
"
Softened
IsothermalPotential:vrot"
])
v_c
=
float
(
params
.
attrs
[
"IsothermalPotential:vrot"
])
v_c_cgs
=
v_c
*
unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
...
...
examples/HydrostaticHalo/makeIC.py
View file @
7fb703af
...
...
@@ -227,7 +227,7 @@ ds[()] = u
u
=
np
.
zeros
(
1
)
# Particle IDs
ids
=
1
+
np
.
linspace
(
0
,
N
,
N
,
endpoint
=
False
,
dtype
=
'L'
)
ids
=
1
+
np
.
linspace
(
0
,
N
,
N
,
endpoint
=
False
)
ds
=
grp
.
create_dataset
(
'ParticleIDs'
,
(
N
,
),
'L'
)
ds
[()]
=
ids
...
...
examples/HydrostaticHalo/velocity_profile.py
View file @
7fb703af
...
...
@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs
=
float
(
params
.
attrs
[
"InternalUnitSystem:UnitLength_in_cgs"
])
unit_velocity_cgs
=
float
(
params
.
attrs
[
"InternalUnitSystem:UnitVelocity_in_cgs"
])
unit_time_cgs
=
unit_length_cgs
/
unit_velocity_cgs
v_c
=
float
(
params
.
attrs
[
"
Softened
IsothermalPotential:vrot"
])
v_c
=
float
(
params
.
attrs
[
"IsothermalPotential:vrot"
])
v_c_cgs
=
v_c
*
unit_velocity_cgs
header
=
f
[
"Header"
]
N
=
header
.
attrs
[
"NumPart_Total"
][
0
]
...
...
examples/IsothermalPotential/makeIC.py
View file @
7fb703af
...
...
@@ -138,7 +138,7 @@ ds = grp1.create_dataset('Masses', (numPart,), 'f')
ds
[()]
=
m
m
=
numpy
.
zeros
(
1
)
ids
=
1
+
numpy
.
linspace
(
0
,
numPart
,
numPart
,
endpoint
=
False
,
dtype
=
'L'
)
ids
=
1
+
numpy
.
linspace
(
0
,
numPart
,
numPart
,
endpoint
=
False
)
ds
=
grp1
.
create_dataset
(
'ParticleIDs'
,
(
numPart
,
),
'L'
)
ds
[()]
=
ids
...
...
examples/Noh_3D/makeIC.py
View file @
7fb703af
...
...
@@ -35,8 +35,8 @@ glass = h5py.File("glassCube_64.hdf5", "r")
vol
=
8.
pos
=
glass
[
"/PartType0/Coordinates"
][:,:]
*
cbrt
(
vol
)
h
=
glass
[
"/PartType0/SmoothingLength"
][:]
*
cbrt
(
vol
)
pos
=
glass
[
"/PartType0/Coordinates"
][:,:]
*
vol
**
(
1.
/
3.
)
h
=
glass
[
"/PartType0/SmoothingLength"
][:]
*
vol
**
(
1.
/
3.
)
numPart
=
size
(
h
)
# Generate extra arrays
...
...
@@ -65,7 +65,7 @@ file = h5py.File(fileName, 'w')
# Header
grp
=
file
.
create_group
(
"/Header"
)
grp
.
attrs
[
"BoxSize"
]
=
[
cbrt
(
vol
),
cbrt
(
vol
),
cbrt
(
vol
)]
grp
.
attrs
[
"BoxSize"
]
=
[
vol
**
(
1.
/
3.
),
vol
**
(
1.
/
3.
),
vol
**
(
1.
/
3.
)]
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
]
...
...
examples/SineWavePotential_1D/makeIC.py
0 → 100644
View file @
7fb703af
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# 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/>.
#
##############################################################################
# Generates a distorted 1D grid with a density profile that balances out the
# external sine wave potential if run with an isothermal equation of state.
import
numpy
as
np
import
h5py
# constant thermal energy
# the definition below assumes the same thermal energy is defined in const.h,
# and also that the code was configured with an adiabatic index of 5./3.
uconst
=
20.2615290634
cs2
=
2.
*
uconst
/
3.
A
=
10.
fileName
=
"sineWavePotential.hdf5"
numPart
=
100
boxSize
=
1.
coords
=
np
.
zeros
((
numPart
,
3
))
v
=
np
.
zeros
((
numPart
,
3
))
m
=
np
.
zeros
(
numPart
)
+
1.
h
=
np
.
zeros
(
numPart
)
+
2.
/
numPart
u
=
np
.
zeros
(
numPart
)
+
uconst
ids
=
np
.
arange
(
numPart
,
dtype
=
'L'
)
rho
=
np
.
zeros
(
numPart
)
# first set the positions, as we try to do a reasonable volume estimate to
# set the masses
for
i
in
range
(
numPart
):
coords
[
i
,
0
]
=
(
i
+
np
.
random
.
random
())
/
numPart
for
i
in
range
(
numPart
):
# reasonable mass estimate (actually not really good, but better than assuming
# a constant volume)
if
i
==
0
:
V
=
0.5
*
(
coords
[
1
,
0
]
-
coords
[
-
1
,
0
]
+
1.
)
elif
i
==
numPart
-
1
:
V
=
0.5
*
(
coords
[
0
,
0
]
+
1.
-
coords
[
-
2
,
0
])
else
:
V
=
0.5
*
(
coords
[
i
+
1
,
0
]
-
coords
[
i
-
1
,
0
])
rho
[
i
]
=
1000.
*
np
.
exp
(
-
0.5
*
A
/
np
.
pi
/
cs2
*
np
.
cos
(
2.
*
np
.
pi
*
coords
[
i
,
0
]))
m
[
i
]
=
rho
[
i
]
*
V
#File
file
=
h5py
.
File
(
fileName
,
'w'
)
# Header
grp
=
file
.
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
[
"MassTable"
]
=
[
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
]
grp
.
attrs
[
"Flag_Entropy_ICs"
]
=
0
grp
.
attrs
[
"Dimension"
]
=
1
#Runtime parameters
grp
=
file
.
create_group
(
"/RuntimePars"
)
grp
.
attrs
[
"PeriodicBoundariesOn"
]
=
1
#Units
grp
=
file
.
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
=
file
.
create_group
(
"/PartType0"
)
grp
.
create_dataset
(
'Coordinates'
,
data
=
coords
,
dtype
=
'd'
)
grp
.
create_dataset
(
'Velocities'
,
data
=
v
,
dtype
=
'f'
)
grp
.
create_dataset
(
'Masses'
,
data
=
m
,
dtype
=
'f'
)
grp
.
create_dataset
(
'SmoothingLength'
,
data
=
h
,
dtype
=
'f'
)
grp
.
create_dataset
(
'InternalEnergy'
,
data
=
u
,
dtype
=
'f'
)
grp
.
create_dataset
(
'ParticleIDs'
,
data
=
ids
,
dtype
=
'L'
)
grp
.
create_dataset
(
'Density'
,
data
=
rho
,
dtype
=
'f'
)
file
.
close
()
examples/SineWavePotential_1D/plotSolution.py
0 → 100644
View file @
7fb703af
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# 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/>.
#
##############################################################################
# Plots some quantities for the snapshot file which is passed on as a command
# line argument (full name)
import
numpy
as
np
import
h5py
import
sys
import
pylab
as
pl
# these should be the same as in makeIC.py
uconst
=
20.2615290634
cs2
=
2.
*
uconst
/
3.
A
=
10.
if
len
(
sys
.
argv
)
<
2
:
print
"Need to provide a filename argument!"
exit
()
fileName
=
sys
.
argv
[
1
]
file
=
h5py
.
File
(
fileName
,
'r'
)
coords
=
np
.
array
(
file
[
"/PartType0/Coordinates"
])
rho
=
np
.
array
(
file
[
"/PartType0/Density"
])
u
=
np
.
array
(
file
[
"/PartType0/InternalEnergy"
])
agrav
=
np
.
array
(
file
[
"/PartType0/GravAcceleration"
])
m
=
np
.
array
(
file
[
"/PartType0/Masses"
])
ids
=
np
.
array
(
file
[
"/PartType0/ParticleIDs"
])
x
=
np
.
linspace
(
0.
,
1.
,
1000
)
rho_x
=
1000.
*
np
.
exp
(
-
0.5
*
A
/
np
.
pi
/
cs2
*
np
.
cos
(
2.
*
np
.
pi
*
x
))
P
=
cs2
*
rho
ids_reverse
=
np
.
argsort
(
ids
)
gradP
=
np
.
zeros
(
P
.
shape
)
for
i
in
range
(
len
(
P
)):
iself
=
int
(
ids
[
i
])
corr
=
0.
im1
=
iself
-
1
if
im1
<
0
:
im1
=
len
(
P
)
-
1
corr
=
1.
ip1
=
iself
+
1
if
ip1
==
len
(
P
):
ip1
=
0
corr
=
1.
idxp1
=
ids_reverse
[
ip1
]
idxm1
=
ids_reverse
[
im1
]
gradP
[
i
]
=
(
P
[
idxp1
]
-
P
[
idxm1
])
/
(
coords
[
idxp1
,
0
]
-
coords
[
idxm1
,
0
])
fig
,
ax
=
pl
.
subplots
(
2
,
2
)
ax
[
0
][
0
].
plot
(
coords
[:,
0
],
rho
,
"r."
,
markersize
=
0.5
)
ax
[
0
][
0
].
plot
(
x
,
rho_x
,
"g-"
)
ax
[
0
][
1
].
plot
(
coords
[:,
0
],
gradP
/
rho
,
"b."
)
ax
[
1
][
0
].
plot
(
coords
[:,
0
],
agrav
[:,
0
],
"g."
,
markersize
=
0.5
)
ax
[
1
][
1
].
plot
(
coords
[:,
0
],
m
,
"y."
)
pl
.
savefig
(
"{fileName}.png"
.
format
(
fileName
=
fileName
[:
-
5
]))
examples/SineWavePotential_1D/run.sh
0 → 100755
View file @
7fb703af
#!/bin/bash
if
[
!
-e
sineWavePotential.hdf5
]
then
echo
"Generating initial conditions for the 1D SineWavePotential example..."
python makeIC.py
fi
../swift
-g
-s
-t
2 sineWavePotential.yml 2>&1 |
tee
output.log
for
f
in
sineWavePotential_
*
.hdf5
do
python plotSolution.py
$f
done
examples/SineWavePotential_1D/sineWavePotential.yml
0 → 100644
View file @
7fb703af
# Define the system of units to use internally.
InternalUnitSystem
:
UnitMass_in_cgs
:
1.
UnitLength_in_cgs
:
1.
UnitVelocity_in_cgs
:
1.
UnitCurrent_in_cgs
:
1.
UnitTemp_in_cgs
:
1.
# Parameters governing the time integration
TimeIntegration
:
time_begin
:
0.
# The starting time of the simulation (in internal units).
time_end
:
10.
# 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 conserved quantities statistics
Statistics
:
delta_time
:
1e-2
# Time between statistics output
# Parameters governing the snapshots
Snapshots
:
basename
:
sineWavePotential
# Common part of the name of output files
time_first
:
0.
# Time of the first output (in internal units)
delta_time
:
1.
# Time difference between consecutive outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH
:
resolution_eta
:
1.2349
# Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
delta_neighbours
:
0.1
# The tolerance for the targetted number of neighbours.
CFL_condition
:
0.1
# Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions
:
file_name
:
sineWavePotential.hdf5
# The file to read
# External potential parameters
SineWavePotential
:
amplitude
:
10.
timestep_limit
:
1.
examples/SineWavePotential_2D/makeIC.py
0 → 100644
View file @
7fb703af
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# 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/>.
#
##############################################################################
# Generates a distorted 1D grid with a density profile that balances out the
# external sine wave potential if run with an isothermal equation of state.
import
numpy
as
np
import
h5py
# constant thermal energy
# the definition below assumes the same thermal energy is defined in const.h,
# and also that the code was configured with an adiabatic index of 5./3.
uconst
=
20.2615290634
cs2
=
2.
*
uconst
/
3.
A
=
10.
fileName
=
"sineWavePotential.hdf5"
numPart_1D
=
50
boxSize
=
[
1.
,
1.
]
numPart
=
numPart_1D
**
2
coords
=
np
.
zeros
((
numPart
,
3
))
v
=
np
.
zeros
((
numPart
,
3
))
m
=
np
.
zeros
(
numPart
)
+
1.
h
=
np
.
zeros
(
numPart
)
+
2.
/
numPart
u
=
np
.
zeros
(
numPart
)
+
uconst
ids
=
np
.
arange
(
numPart
,
dtype
=
'L'
)
rho
=
np
.
zeros
(
numPart
)
# first set the positions, as we try to do a reasonable volume estimate to
# set the masses
for
i
in
range
(
numPart_1D
):
for
j
in
range
(
numPart_1D
):
coords
[
numPart_1D
*
i
+
j
,
0
]
=
(
i
+
0.5
)
/
numPart_1D
coords
[
numPart_1D
*
i
+
j
,
1
]
=
(
j
+
0.5
)
/
numPart_1D
V
=
1.
/
numPart
for
i
in
range
(
numPart
):
rho
[
i
]
=
1000.
*
np
.
exp
(
-
0.5
*
A
/
np
.
pi
/
cs2
*
np
.
cos
(
2.
*
np
.
pi
*
coords
[
i
,
0
]))
m
[
i
]
=
rho
[
i
]
*
V
#File
file
=
h5py
.
File
(
fileName
,
'w'
)
# Header
grp
=
file
.
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
[
"MassTable"
]
=
[
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
]
grp
.
attrs
[
"Flag_Entropy_ICs"
]
=
0
grp
.
attrs
[
"Dimension"
]
=
2
#Runtime parameters
grp
=
file
.
create_group
(
"/RuntimePars"
)
grp
.
attrs
[
"PeriodicBoundariesOn"
]
=
1
#Units
grp
=
file
.
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
=
file
.
create_group
(
"/PartType0"
)
grp
.
create_dataset
(
'Coordinates'
,
data
=
coords
,
dtype
=
'd'
)
grp
.
create_dataset
(
'Velocities'
,
data
=
v
,
dtype
=
'f'
)
grp
.
create_dataset
(
'Masses'
,
data
=
m
,
dtype
=
'f'
)
grp
.
create_dataset
(
'SmoothingLength'
,
data
=
h
,
dtype
=
'f'
)
grp
.
create_dataset
(
'InternalEnergy'
,
data
=
u
,
dtype
=
'f'
)
grp
.
create_dataset
(
'ParticleIDs'
,
data
=
ids
,
dtype
=
'L'
)
grp
.
create_dataset
(
'Density'
,
data
=
rho
,
dtype
=
'f'
)
file
.
close
()
examples/SineWavePotential_2D/plotSolution.py
0 → 100644
View file @
7fb703af
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# 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/>.
#
##############################################################################
# Plots some quantities for the snapshot file which is passed on as a command
# line argument (full name)
import
numpy
as
np
import
h5py
import
sys
import
pylab
as
pl
# these should be the same as in makeIC.py
uconst
=
20.2615290634
cs2
=
2.
*
uconst
/
3.
A
=
10.
if
len
(
sys
.
argv
)
<
2
:
print
"Need to provide a filename argument!"
exit
()
fileName
=
sys
.
argv
[
1
]
file
=
h5py
.
File
(
fileName
,
'r'
)
coords
=
np
.
array
(
file
[
"/PartType0/Coordinates"
])
rho
=
np
.
array
(
file
[
"/PartType0/Density"
])
u
=
np
.
array
(
file
[
"/PartType0/InternalEnergy"
])
agrav
=
np
.
array
(
file
[
"/PartType0/GravAcceleration"
])
m
=
np
.
array
(
file
[
"/PartType0/Masses"
])
ids
=
np
.
array
(
file
[
"/PartType0/ParticleIDs"
])
# ids_reverse gives the index original particle 0 now has in the particle arrays
# and so on
# note that this will only work if the particles do not move away too much from
# there initial positions
ids_reverse
=
np
.
argsort
(
ids
)
x
=
np
.
linspace
(
0.
,
1.
,
1000
)
rho_x
=
1000.
*
np
.
exp
(
-
0.5
*
A
/
np
.
pi
/
cs2
*
np
.
cos
(
2.
*
np
.
pi
*
x
))
P
=
cs2
*
rho
n1D
=
int
(
np
.
sqrt
(
len
(
P
)))
gradP
=
np
.
zeros
(
P
.
shape
)
for
i
in
range
(
len
(
P
)):
iself
=
int
(
ids
[
i
]
/
n1D
)
jself
=
int
(
ids
[
i
]
-
n1D
*
iself
)
corr
=
0.
im1
=
iself
-
1
if
im1
<
0
: