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
cdad038f
Commit
cdad038f
authored
Mar 03, 2017
by
Bert Vandenbroucke
Browse files
Added a sine wave potential and corresponding 1D hydro + external gravity test case.
parent
ab0491e9
Changes
7
Hide whitespace changes
Inline
Side-by-side
configure.ac
View file @
cdad038f
...
...
@@ -751,7 +751,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 +769,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])
;;
...
...
examples/SineWavePotential_1D/makeIC.py
0 → 100644
View file @
cdad038f
###############################################################################
# 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
)
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 @
cdad038f
###############################################################################
# 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"
])
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
gradP
=
np
.
zeros
(
P
.
shape
)
for
i
in
range
(
len
(
P
)):
im1
=
i
-
1
if
im1
<
0
:
im1
=
len
(
P
)
-
1
ip1
=
i
+
1
if
ip1
==
len
(
P
):
ip1
=
0
gradP
[
i
]
=
(
P
[
ip1
]
-
P
[
im1
])
/
(
coords
[
2
,
0
]
-
coords
[
0
,
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 @
cdad038f
#!/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 @
cdad038f
# 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.
src/potential.h
View file @
cdad038f
...
...
@@ -36,6 +36,8 @@
#include
"./potential/isothermal/potential.h"
#elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
#include
"./potential/disc_patch/potential.h"
#elif defined(EXTERNAL_POTENTIAL_SINE_WAVE)
#include
"./potential/sine_wave/potential.h"
#else
#error "Invalid choice of external potential"
#endif
...
...
src/potential/sine_wave/potential.h
0 → 100644
View file @
cdad038f
/*******************************************************************************
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_SINE_WAVE_H
#define SWIFT_SINE_WAVE_H
/* Config parameters. */
#include
"../config.h"
/* Some standard headers. */
#include
<float.h>
#include
<math.h>
/* Local includes. */
#include
"const.h"
#include
"error.h"
#include
"parser.h"
#include
"part.h"
#include
"physical_constants.h"
#include
"space.h"
#include
"units.h"
/**
* @brief External Potential Properties - Sine wave case
*/
struct
external_potential
{
/*! Amplitude of the sine wave. */
double
amplitude
;
/*! Time-step limiting factor. */
double
timestep_limit
;
};
/**
* @brief Computes the time-step from the acceleration due to a sine wave.
*
* @param time The current time.
* @param potential The properties of the potential.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__
((
always_inline
))
INLINE
static
float
external_gravity_timestep
(
double
time
,
const
struct
external_potential
*
restrict
potential
,
const
struct
phys_const
*
restrict
phys_const
,
const
struct
gpart
*
restrict
g
)
{
return
potential
->
timestep_limit
;
}
/**
* @brief Computes the gravitational acceleration along x given by the sine
* wave.
*
* @param time The current time in internal units.
* @param potential The properties of the potential.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__
((
always_inline
))
INLINE
static
void
external_gravity_acceleration
(
double
time
,
const
struct
external_potential
*
restrict
potential
,
const
struct
phys_const
*
restrict
phys_const
,
struct
gpart
*
restrict
g
)
{
g
->
a_grav
[
0
]
=
potential
->
amplitude
*
sin
(
2
.
*
M_PI
*
g
->
x
[
0
])
/
phys_const
->
const_newton_G
;
}
/**
* @brief Computes the gravitational potential energy of a particle in the
* sine wave.
*
* @param time The current time.
* @param potential The #external_potential used in the run.
* @param phys_const Physical constants in internal units.
* @param gp Pointer to the particle data.
*/
__attribute__
((
always_inline
))
INLINE
static
float
external_gravity_get_potential_energy
(
double
time
,
const
struct
external_potential
*
potential
,
const
struct
phys_const
*
const
phys_const
,
const
struct
gpart
*
gp
)
{
/* this potential does not really have a potential energy */
return
0
.;
}
/**
* @brief Initialises the external potential properties in the internal system
* of units.
*
* @param parameter_file The parsed parameter file
* @param phys_const Physical constants in internal units
* @param us The current internal system of units
* @param potential The external potential properties to initialize
*/
static
INLINE
void
potential_init_backend
(
const
struct
swift_params
*
parameter_file
,
const
struct
phys_const
*
phys_const
,
const
struct
unit_system
*
us
,
const
struct
space
*
s
,
struct
external_potential
*
potential
)
{
potential
->
amplitude
=
parser_get_param_double
(
parameter_file
,
"SineWavePotential:amplitude"
);
potential
->
timestep_limit
=
parser_get_param_double
(
parameter_file
,
"SineWavePotential:timestep_limit"
);
}
/**
* @brief Prints the properties of the external potential to stdout.
*
* @param potential The external potential properties.
*/
static
INLINE
void
potential_print_backend
(
const
struct
external_potential
*
potential
)
{
message
(
"External potential is a sine wave with amplitude %g"
,
potential
->
amplitude
);
}
#endif
/* SWIFT_SINE_WAVE_H */
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