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
e240f03b
Commit
e240f03b
authored
Jul 25, 2017
by
Bert Vandenbroucke
Browse files
Switched disc patch potential direction from z to x to be able to run it in 1D and 2D.
parent
ffbf1f31
Changes
5
Hide whitespace changes
Inline
Side-by-side
examples/DiscPatch/HydroStatic/disc-patch-icc.yml
View file @
e240f03b
...
...
@@ -39,8 +39,8 @@ InitialConditions:
DiscPatchPotential
:
surface_density
:
10.
scale_height
:
100.
z
_disc
:
400.
z
_trunc
:
300.
z
_max
:
350.
x
_disc
:
400.
x
_trunc
:
300.
x
_max
:
350.
timestep_mult
:
0.03
growth_time
:
5.
examples/DiscPatch/HydroStatic/disc-patch.yml
View file @
e240f03b
...
...
@@ -39,7 +39,7 @@ InitialConditions:
DiscPatchPotential
:
surface_density
:
10.
scale_height
:
100.
z
_disc
:
400.
z
_trunc
:
300.
z
_max
:
380.
x
_disc
:
400.
x
_trunc
:
300.
x
_max
:
380.
timestep_mult
:
0.03
examples/DiscPatch/HydroStatic/makeIC.py
View file @
e240f03b
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
# Tom Theuns (tom.theuns@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 John A. Regan (john.a.regan@durham.ac.uk)
# Tom Theuns (tom.theuns@durham.ac.uk)
# 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# 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/>.
#
##############################################################################
import
h5py
import
sys
...
...
@@ -37,16 +39,16 @@ import random
# Size of the patch -- side_length
# Parameters of the gas disc
surface_density
=
10.
surface_density
=
10.
scale_height
=
100.
gas_gamma
=
5.
/
3.
# Parameters of the problem
z
_factor
=
2
x
_factor
=
2
side_length
=
400.
# File
fileName
=
"Disc-Patch.hdf5"
fileName
=
"Disc-Patch.hdf5"
####################################################################
...
...
@@ -64,30 +66,33 @@ unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
unit_velocity_in_cgs
=
(
1e5
)
unit_time_in_cgs
=
unit_length_in_cgs
/
unit_velocity_in_cgs
print
"UnitMass_in_cgs: %.5e"
%
unit_mass_in_cgs
print
"UnitMass_in_cgs: %.5e"
%
unit_mass_in_cgs
print
"UnitLength_in_cgs: %.5e"
%
unit_length_in_cgs
print
"UnitVelocity_in_cgs: %.5e"
%
unit_velocity_in_cgs
print
"UnitTime_in_cgs: %.5e"
%
unit_time_in_cgs
print
""
# Derived units
const_G
=
NEWTON_GRAVITY_CGS
*
unit_mass_in_cgs
*
unit_time_in_cgs
**
2
*
unit_length_in_cgs
**-
3
const_G
=
NEWTON_GRAVITY_CGS
*
unit_mass_in_cgs
*
unit_time_in_cgs
**
2
*
\
unit_length_in_cgs
**-
3
const_mp
=
PROTON_MASS_IN_CGS
*
unit_mass_in_cgs
**-
1
const_kb
=
BOLTZMANN_IN_CGS
*
unit_mass_in_cgs
**-
1
*
unit_length_in_cgs
**-
2
*
unit_time_in_cgs
**
2
const_kb
=
BOLTZMANN_IN_CGS
*
unit_mass_in_cgs
**-
1
*
unit_length_in_cgs
**-
2
*
\
unit_time_in_cgs
**
2
print
"--- Some constants [internal units] ---"
print
"G_Newton:
%.5e"
%
const_G
print
"m_proton:
%.5e"
%
const_mp
print
"k_boltzmann:
%.5e"
%
const_kb
print
"G_Newton: %.5e"
%
const_G
print
"m_proton: %.5e"
%
const_mp
print
"k_boltzmann: %.5e"
%
const_kb
print
""
# derived quantities
temp
=
math
.
pi
*
const_G
*
surface_density
*
scale_height
*
const_mp
/
const_kb
u_therm
=
const_kb
*
temp
/
((
gas_gamma
-
1
)
*
const_mp
)
v_disp
=
math
.
sqrt
(
2
*
u_therm
)
soundspeed
=
math
.
sqrt
(
u_therm
/
(
gas_gamma
*
(
gas_gamma
-
1.
)))
t_dyn
=
math
.
sqrt
(
scale_height
/
(
const_G
*
surface_density
))
t_cross
=
scale_height
/
soundspeed
temp
=
math
.
pi
*
const_G
*
surface_density
*
scale_height
*
const_mp
/
\
const_kb
u_therm
=
const_kb
*
temp
/
((
gas_gamma
-
1
)
*
const_mp
)
v_disp
=
math
.
sqrt
(
2
*
u_therm
)
soundspeed
=
math
.
sqrt
(
u_therm
/
(
gas_gamma
*
(
gas_gamma
-
1.
)))
t_dyn
=
math
.
sqrt
(
scale_height
/
(
const_G
*
surface_density
))
t_cross
=
scale_height
/
soundspeed
print
"--- Properties of the gas [internal units] ---"
print
"Gas temperature: %.5e"
%
temp
...
...
@@ -101,18 +106,20 @@ print ""
# Problem properties
boxSize_x
=
side_length
boxSize_y
=
boxSize_x
boxSize_z
=
boxSize_x
*
z_factor
boxSize_z
=
boxSize_x
boxSize_x
*=
x_factor
volume
=
boxSize_x
*
boxSize_y
*
boxSize_z
M_tot
=
boxSize_x
*
boxSize_y
*
surface_density
*
math
.
tanh
(
boxSize_z
/
(
2.
*
scale_height
))
M_tot
=
boxSize_y
*
boxSize_z
*
surface_density
*
\
math
.
tanh
(
boxSize_x
/
(
2.
*
scale_height
))
density
=
M_tot
/
volume
entropy
=
(
gas_gamma
-
1.
)
*
u_therm
/
density
**
(
gas_gamma
-
1.
)
print
"--- Problem properties [internal units] ---"
print
"Box:
[%.1f, %.1f, %.1f]"
%
(
boxSize_x
,
boxSize_y
,
boxSize_z
)
print
"Volume:
%.5e"
%
volume
print
"Total mass:
%.5e"
%
M_tot
print
"Density:
%.5e"
%
density
print
"Entropy:
%.5e"
%
entropy
print
"Box: [%.1f, %.1f, %.1f]"
%
(
boxSize_x
,
boxSize_y
,
boxSize_z
)
print
"Volume: %.5e"
%
volume
print
"Total mass: %.5e"
%
M_tot
print
"Density: %.5e"
%
density
print
"Entropy: %.5e"
%
entropy
print
""
####################################################################
...
...
@@ -123,34 +130,24 @@ one_glass_pos = infile["/PartType0/Coordinates"][:,:]
one_glass_h
=
infile
[
"/PartType0/SmoothingLength"
][:]
# Rescale to the problem size
one_glass_pos
*=
boxSize_x
one_glass_h
*=
boxSize_x
#print min(one_glass_p[:,0]), max(one_glass_p[:,0])
#print min(one_glass_p[:,1]), max(one_glass_p[:,1])
#print min(one_glass_p[:,2]), max(one_glass_p[:,2])
one_glass_pos
*=
side_length
one_glass_h
*=
side_length
# Now create enough copies to fill the volume in
z
# Now create enough copies to fill the volume in
x
pos
=
np
.
copy
(
one_glass_pos
)
h
=
np
.
copy
(
one_glass_h
)
for
i
in
range
(
1
,
z_factor
):
one_glass_pos
[:,
2
]
+=
boxSize_x
for
i
in
range
(
1
,
x_factor
):
one_glass_pos
[:,
0
]
+=
side_length
pos
=
np
.
append
(
pos
,
one_glass_pos
,
axis
=
0
)
h
=
np
.
append
(
h
,
one_glass_h
,
axis
=
0
)
#print min(pos[:,0]), max(pos[:,0])
#print min(pos[:,1]), max(pos[:,1])
#print min(pos[:,2]), max(pos[:,2])
# Compute further properties of ICs
numPart
=
np
.
size
(
h
)
mass
=
M_tot
/
numPart
print
"--- Particle properties [internal units] ---"
print
"Number part.: "
,
numPart
print
"Part. mass:
%.5e"
%
mass
print
"Part. mass: %.5e"
%
mass
print
""
# Create additional arrays
...
...
@@ -159,7 +156,6 @@ mass = np.ones(numPart) * mass
vel
=
np
.
zeros
((
numPart
,
3
))
ids
=
1
+
np
.
linspace
(
0
,
numPart
,
numPart
,
endpoint
=
False
)
####################################################################
# Create and write output file
...
...
@@ -169,7 +165,7 @@ file = h5py.File(fileName, 'w')
#Units
grp
=
file
.
create_group
(
"/Units"
)
grp
.
attrs
[
"Unit length in cgs (U_L)"
]
=
unit_length_in_cgs
grp
.
attrs
[
"Unit mass in cgs (U_M)"
]
=
unit_mass_in_cgs
grp
.
attrs
[
"Unit mass in cgs (U_M)"
]
=
unit_mass_in_cgs
grp
.
attrs
[
"Unit time in cgs (U_t)"
]
=
unit_time_in_cgs
grp
.
attrs
[
"Unit current in cgs (U_I)"
]
=
1.
grp
.
attrs
[
"Unit temperature in cgs (U_T)"
]
=
1.
...
...
@@ -200,13 +196,12 @@ ds = grp0.create_dataset('SmoothingLength', (numPart,), 'f', data=h)
ds
=
grp0
.
create_dataset
(
'InternalEnergy'
,
(
numPart
,),
'f'
,
data
=
u
)
ds
=
grp0
.
create_dataset
(
'ParticleIDs'
,
(
numPart
,
),
'L'
,
data
=
ids
)
####################################################################
print
"--- Runtime parameters (YAML file): ---"
print
"DiscPatchPotential:surface_density: "
,
surface_density
print
"DiscPatchPotential:scale_height: "
,
scale_height
print
"DiscPatchPotential:
z
_disc: "
,
boxSize_
z
/
2.
print
"DiscPatchPotential:
x
_disc: "
,
0.5
*
boxSize_
x
print
""
print
"--- Constant parameters: ---"
...
...
examples/DiscPatch/HydroStatic/plotSolution.py
View file @
e240f03b
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
# 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
...
...
@@ -34,9 +35,9 @@ import sys
# Parameters
surface_density
=
10.
scale_height
=
100.
z
_disc
=
400.
z
_trunc
=
300.
z
_max
=
350.
x
_disc
=
400.
x
_trunc
=
300.
x
_max
=
350.
utherm
=
20.2678457288
gamma
=
5.
/
3.
...
...
@@ -50,14 +51,14 @@ if len(sys.argv) > 2:
# Get the analytic solution for the density
def
get_analytic_density
(
x
):
return
0.5
*
surface_density
/
scale_height
/
\
np
.
cosh
(
(
x
-
z
_disc
)
/
scale_height
)
**
2
np
.
cosh
(
(
x
-
x
_disc
)
/
scale_height
)
**
2
# Get the analytic solution for the (isothermal) pressure
def
get_analytic_pressure
(
x
):
return
(
gamma
-
1.
)
*
utherm
*
get_analytic_density
(
x
)
# Get the data fields to plot from the snapshot file with the given name:
# snapshot time,
z
-coord, density, pressure, velocity norm
# snapshot time,
x
-coord, density, pressure, velocity norm
def
get_data
(
name
):
file
=
h5py
.
File
(
name
,
"r"
)
coords
=
np
.
array
(
file
[
"/PartType0/Coordinates"
])
...
...
@@ -69,7 +70,7 @@ def get_data(name):
vtot
=
np
.
sqrt
(
v
[:,
0
]
**
2
+
v
[:,
1
]
**
2
+
v
[:,
2
]
**
2
)
return
float
(
file
[
"/Header"
].
attrs
[
"Time"
]),
coords
[:,
2
],
rho
,
P
,
vtot
return
float
(
file
[
"/Header"
].
attrs
[
"Time"
]),
coords
[:,
0
],
rho
,
P
,
vtot
# scan the folder for snapshot files and plot all of them (within the requested
# range)
...
...
@@ -80,38 +81,38 @@ for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
print
"processing"
,
f
,
"..."
z
range
=
np
.
linspace
(
0.
,
2.
*
z
_disc
,
1000
)
time
,
z
,
rho
,
P
,
v
=
get_data
(
f
)
x
range
=
np
.
linspace
(
0.
,
2.
*
x
_disc
,
1000
)
time
,
x
,
rho
,
P
,
v
=
get_data
(
f
)
fig
,
ax
=
pl
.
subplots
(
3
,
1
,
sharex
=
True
)
ax
[
0
].
plot
(
z
,
rho
,
"r."
)
ax
[
0
].
plot
(
z
range
,
get_analytic_density
(
z
range
),
"k-"
)
ax
[
0
].
plot
([
z
_disc
-
z
_max
,
z
_disc
-
z
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
plot
([
z
_disc
+
z
_max
,
z
_disc
+
z
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
plot
([
z
_disc
-
z
_trunc
,
z
_disc
-
z
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
plot
([
z
_disc
+
z
_trunc
,
z
_disc
+
z
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
set_ylim
(
0.
,
1.2
*
get_analytic_density
(
z
_disc
))
ax
[
0
].
plot
(
x
,
rho
,
"r."
)
ax
[
0
].
plot
(
x
range
,
get_analytic_density
(
x
range
),
"k-"
)
ax
[
0
].
plot
([
x
_disc
-
x
_max
,
x
_disc
-
x
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
plot
([
x
_disc
+
x
_max
,
x
_disc
+
x
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
plot
([
x
_disc
-
x
_trunc
,
x
_disc
-
x
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
plot
([
x
_disc
+
x
_trunc
,
x
_disc
+
x
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
0
].
set_ylim
(
0.
,
1.2
*
get_analytic_density
(
x
_disc
))
ax
[
0
].
set_ylabel
(
"density"
)
ax
[
1
].
plot
(
z
,
v
,
"r."
)
ax
[
1
].
plot
(
z
range
,
np
.
zeros
(
len
(
z
range
)),
"k-"
)
ax
[
1
].
plot
([
z
_disc
-
z
_max
,
z
_disc
-
z
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
plot
([
z
_disc
+
z
_max
,
z
_disc
+
z
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
plot
([
z
_disc
-
z
_trunc
,
z
_disc
-
z
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
plot
([
z
_disc
+
z
_trunc
,
z
_disc
+
z
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
plot
(
x
,
v
,
"r."
)
ax
[
1
].
plot
(
x
range
,
np
.
zeros
(
len
(
x
range
)),
"k-"
)
ax
[
1
].
plot
([
x
_disc
-
x
_max
,
x
_disc
-
x
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
plot
([
x
_disc
+
x
_max
,
x
_disc
+
x
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
plot
([
x
_disc
-
x
_trunc
,
x
_disc
-
x
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
plot
([
x
_disc
+
x
_trunc
,
x
_disc
+
x
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
1
].
set_ylim
(
-
0.5
,
10.
)
ax
[
1
].
set_ylabel
(
"velocity norm"
)
ax
[
2
].
plot
(
z
,
P
,
"r."
)
ax
[
2
].
plot
(
z
range
,
get_analytic_pressure
(
z
range
),
"k-"
)
ax
[
2
].
plot
([
z
_disc
-
z
_max
,
z
_disc
-
z
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
plot
([
z
_disc
+
z
_max
,
z
_disc
+
z
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
plot
([
z
_disc
-
z
_trunc
,
z
_disc
-
z
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
plot
([
z
_disc
+
z
_trunc
,
z
_disc
+
z
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
set_xlim
(
0.
,
2.
*
z
_disc
)
ax
[
2
].
set_ylim
(
0.
,
1.2
*
get_analytic_pressure
(
z
_disc
))
ax
[
2
].
set_xlabel
(
"
z
"
)
ax
[
2
].
plot
(
x
,
P
,
"r."
)
ax
[
2
].
plot
(
x
range
,
get_analytic_pressure
(
x
range
),
"k-"
)
ax
[
2
].
plot
([
x
_disc
-
x
_max
,
x
_disc
-
x
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
plot
([
x
_disc
+
x
_max
,
x
_disc
+
x
_max
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
plot
([
x
_disc
-
x
_trunc
,
x
_disc
-
x
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
plot
([
x
_disc
+
x
_trunc
,
x
_disc
+
x
_trunc
],
[
0
,
10
],
"k--"
,
alpha
=
0.5
)
ax
[
2
].
set_xlim
(
0.
,
2.
*
x
_disc
)
ax
[
2
].
set_ylim
(
0.
,
1.2
*
get_analytic_pressure
(
x
_disc
))
ax
[
2
].
set_xlabel
(
"
x
"
)
ax
[
2
].
set_ylabel
(
"pressure"
)
pl
.
suptitle
(
"t = {0:.2f}"
.
format
(
time
))
...
...
src/potential/disc_patch/potential.h
View file @
e240f03b
...
...
@@ -56,20 +56,20 @@ struct external_potential {
/*! Inverse of disc scale-height (1/b) */
float
scale_height_inv
;
/*! Position of the disc along the
z
-axis */
float
z
_disc
;
/*! Position of the disc along the
x
-axis */
float
x
_disc
;
/*! Position above which the accelerations get truncated */
float
z
_trunc
;
float
x
_trunc
;
/*! Position above which the accelerations are zero */
float
z
_max
;
float
x
_max
;
/*! The truncated transition regime */
float
z
_trans
;
float
x
_trans
;
/*! Inverse of the truncated transition regime */
float
z
_trans_inv
;
float
x
_trans_inv
;
/*! Dynamical time of the system */
float
dynamical_time
;
...
...
@@ -115,36 +115,36 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const
float
norm
=
potential
->
norm
;
/* absolute value of height above disc */
const
float
d
z
=
fabsf
(
g
->
x
[
2
]
-
potential
->
z
_disc
);
const
float
d
x
=
fabsf
(
g
->
x
[
0
]
-
potential
->
x
_disc
);
/* vertical acceleration */
const
float
z
_accel
=
norm
*
tanhf
(
d
z
*
b_inv
);
const
float
x
_accel
=
norm
*
tanhf
(
d
x
*
b_inv
);
float
dt
=
dt_dyn
;
/* demand that dt * velocity < fraction of scale height of disc */
if
(
g
->
v_full
[
2
]
!=
0
.
f
)
{
if
(
g
->
v_full
[
0
]
!=
0
.
f
)
{
const
float
dt1
=
b
/
fabsf
(
g
->
v_full
[
2
]);
const
float
dt1
=
b
/
fabsf
(
g
->
v_full
[
0
]);
dt
=
min
(
dt1
,
dt
);
}
/* demand that dt^2 * acceleration < fraction of scale height of disc */
if
(
z
_accel
!=
0
.
f
)
{
if
(
x
_accel
!=
0
.
f
)
{
const
float
dt2
=
b
/
fabsf
(
z
_accel
);
const
float
dt2
=
b
/
fabsf
(
x
_accel
);
if
(
dt2
<
dt
*
dt
)
dt
=
sqrtf
(
dt2
);
}
/* demand that dt^3 * jerk < fraction of scale height of disc */
if
(
g
->
v_full
[
2
]
!=
0
.
f
)
{
if
(
g
->
v_full
[
0
]
!=
0
.
f
)
{
const
float
cosh_d
z
_inv
=
1
.
f
/
coshf
(
d
z
*
b_inv
);
const
float
cosh_d
z
_inv2
=
cosh_d
z
_inv
*
cosh_d
z
_inv
;
const
float
d
z
_accel_over_dt
=
norm
*
cosh_d
z
_inv2
*
b_inv
*
fabsf
(
g
->
v_full
[
2
]);
const
float
cosh_d
x
_inv
=
1
.
f
/
coshf
(
d
x
*
b_inv
);
const
float
cosh_d
x
_inv2
=
cosh_d
x
_inv
*
cosh_d
x
_inv
;
const
float
d
x
_accel_over_dt
=
norm
*
cosh_d
x
_inv2
*
b_inv
*
fabsf
(
g
->
v_full
[
0
]);
const
float
dt3
=
b
/
fabsf
(
d
z
_accel_over_dt
);
const
float
dt3
=
b
/
fabsf
(
d
x
_accel_over_dt
);
if
(
dt3
<
dt
*
dt
*
dt
)
dt
=
cbrtf
(
dt3
);
}
...
...
@@ -152,13 +152,13 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
}
/**
* @brief Computes the gravitational acceleration along
z
due to a hydrostatic
* @brief Computes the gravitational acceleration along
x
due to a hydrostatic
* disc
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 17.
* We truncate the accelerations beyond
z
_trunc using a 1-cos(
z
) function
* that smoothly brings the accelerations to 0 at
z
_max.
* We truncate the accelerations beyond
x
_trunc using a 1-cos(
x
) function
* that smoothly brings the accelerations to 0 at
x
_max.
*
* @param time The current time in internal units.
* @param potential The properties of the potential.
...
...
@@ -169,40 +169,40 @@ __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
)
{
const
float
d
z
=
g
->
x
[
2
]
-
potential
->
z
_disc
;
const
float
abs_d
z
=
fabsf
(
d
z
);
const
float
d
x
=
g
->
x
[
0
]
-
potential
->
x
_disc
;
const
float
abs_d
x
=
fabsf
(
d
x
);
const
float
t_growth
=
potential
->
growth_time
;
const
float
t_growth_inv
=
potential
->
growth_time_inv
;
const
float
b_inv
=
potential
->
scale_height_inv
;
const
float
z
_trunc
=
potential
->
z
_trunc
;
const
float
z
_max
=
potential
->
z
_max
;
const
float
z
_trans_inv
=
potential
->
z
_trans_inv
;
const
float
x
_trunc
=
potential
->
x
_trunc
;
const
float
x
_max
=
potential
->
x
_max
;
const
float
x
_trans_inv
=
potential
->
x
_trans_inv
;
const
float
norm_over_G
=
potential
->
norm_over_G
;
/* Are we still growing the disc ? */
const
float
reduction_factor
=
time
<
t_growth
?
time
*
t_growth_inv
:
1
.
f
;
/* Truncated or not ? */
float
a_
z
;
if
(
abs_d
z
<
z
_trunc
)
{
float
a_
x
;
if
(
abs_d
x
<
x
_trunc
)
{
/* Acc. 2 pi sigma tanh(
z
/b) */
a_
z
=
reduction_factor
*
norm_over_G
*
tanhf
(
abs_d
z
*
b_inv
);
}
else
if
(
abs_d
z
<
z
_max
)
{
/* Acc. 2 pi sigma tanh(
x
/b) */
a_
x
=
reduction_factor
*
norm_over_G
*
tanhf
(
abs_d
x
*
b_inv
);
}
else
if
(
abs_d
x
<
x
_max
)
{
/* Acc. 2 pi sigma tanh(
z
/b) [1/2 + 1/2cos((
z-z
max)/(pi
z
_trans))] */
a_
z
=
reduction_factor
*
norm_over_G
*
tanhf
(
abs_d
z
*
b_inv
)
*
(
0
.
5
f
+
0
.
5
f
*
cosf
((
float
)(
M_PI
)
*
(
abs_d
z
-
z
_trunc
)
*
z
_trans_inv
));
/* Acc. 2 pi sigma tanh(
x
/b) [1/2 + 1/2cos((
x-x
max)/(pi
x
_trans))] */
a_
x
=
reduction_factor
*
norm_over_G
*
tanhf
(
abs_d
x
*
b_inv
)
*
(
0
.
5
f
+
0
.
5
f
*
cosf
((
float
)(
M_PI
)
*
(
abs_d
x
-
x
_trunc
)
*
x
_trans_inv
));
}
else
{
/* Acc. 0 */
a_
z
=
0
.
f
;
a_
x
=
0
.
f
;
}
/* Get the correct sign. Recall G is multipiled in later on */
if
(
d
z
>
0
)
g
->
a_grav
[
2
]
-=
a_
z
;
if
(
d
z
<
0
)
g
->
a_grav
[
2
]
+=
a_
z
;
if
(
d
x
>
0
)
g
->
a_grav
[
0
]
-=
a_
x
;
if
(
d
x
<
0
)
g
->
a_grav
[
0
]
+=
a_
x
;
}
/**
...
...
@@ -211,8 +211,8 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 22.
* We truncate the accelerations beyond
z
_trunc using a 1-cos(
z
) function
* that smoothly brings the accelerations to 0 at
z
_max.
* We truncate the accelerations beyond
x
_trunc using a 1-cos(
x
) function
* that smoothly brings the accelerations to 0 at
x
_max.
*
* @param time The current time.
* @param potential The #external_potential used in the run.
...
...
@@ -224,28 +224,28 @@ external_gravity_get_potential_energy(
double
time
,
const
struct
external_potential
*
potential
,
const
struct
phys_const
*
const
phys_const
,
const
struct
gpart
*
gp
)
{
const
float
d
z
=
gp
->
x
[
2
]
-
potential
->
z
_disc
;
const
float
abs_d
z
=
fabsf
(
d
z
);
const
float
d
x
=
gp
->
x
[
0
]
-
potential
->
x
_disc
;
const
float
abs_d
x
=
fabsf
(
d
x
);
const
float
t_growth
=
potential
->
growth_time
;
const
float
t_growth_inv
=
potential
->
growth_time_inv
;
const
float
b
=
potential
->
scale_height
;
const
float
b_inv
=
potential
->
scale_height_inv
;
const
float
norm
=
potential
->
norm
;
const
float
z
_trunc
=
potential
->
z
_trunc
;
const
float
z
_max
=
potential
->
z
_max
;
const
float
x
_trunc
=
potential
->
x
_trunc
;
const
float
x
_max
=
potential
->
x
_max
;
/* Are we still growing the disc ? */
const
float
reduction_factor
=
time
<
t_growth
?
time
*
t_growth_inv
:
1
.
f
;
/* Truncated or not ? */
float
pot
;
if
(
abs_d
z
<
z
_trunc
)
{
if
(
abs_d
x
<
x
_trunc
)
{
/* Potential (2 pi G sigma b ln(cosh(
z
/b)) */
pot
=
b
*
logf
(
coshf
(
d
z
*
b_inv
));
}
else
if
(
abs_d
z
<
z
_max
)
{
/* Potential (2 pi G sigma b ln(cosh(
x
/b)) */
pot
=
b
*
logf
(
coshf
(
d
x
*
b_inv
));
}
else
if
(
abs_d
x
<
x
_max
)
{
/* Potential. At
z
>>b, phi(
z
) = norm *
z
/ b */
/* Potential. At
x
>>b, phi(
x
) = norm *
x
/ b */
pot
=
0
.
f
;
}
else
{
...
...
@@ -277,14 +277,14 @@ static INLINE void potential_init_backend(
parameter_file
,
"DiscPatchPotential:surface_density"
);
potential
->
scale_height
=
parser_get_param_double
(
parameter_file
,
"DiscPatchPotential:scale_height"
);
potential
->
z
_disc
=
parser_get_param_double
(
parameter_file
,
"DiscPatchPotential:
z
_disc"
);
potential
->
z
_trunc
=
parser_get_opt_param_double
(
parameter_file
,
"DiscPatchPotential:
z
_trunc"
,
FLT_MAX
);
potential
->
z
_max
=
parser_get_opt_param_double
(
parameter_file
,
"DiscPatchPotential:
z
_max"
,
FLT_MAX
);
potential
->
z
_disc
=
parser_get_param_double
(
parameter_file
,
"DiscPatchPotential:
z
_disc"
);
potential
->
x
_disc
=
parser_get_param_double
(
parameter_file
,
"DiscPatchPotential:
x
_disc"
);
potential
->
x
_trunc
=
parser_get_opt_param_double
(
parameter_file
,
"DiscPatchPotential:
x
_trunc"
,
FLT_MAX
);
potential
->
x
_max
=
parser_get_opt_param_double
(
parameter_file
,
"DiscPatchPotential:
x
_max"
,
FLT_MAX
);
potential
->
x
_disc
=
parser_get_param_double
(
parameter_file
,
"DiscPatchPotential:
x
_disc"
);
potential
->
timestep_mult
=
parser_get_param_double
(
parameter_file
,
"DiscPatchPotential:timestep_mult"
);
potential
->
growth_time
=
parser_get_opt_param_double
(
...
...
@@ -299,22 +299,22 @@ static INLINE void potential_init_backend(
potential
->
growth_time
*=
potential
->
dynamical_time
;
/* Some cross-checks */
if
(
potential
->
z
_trunc
>
potential
->
z
_max
)
error
(
"Potential truncation
z
larger than maximal z"
);
if
(
potential
->
z
_trunc
<
potential
->
scale_height
)
error
(
"Potential truncation
z
smaller than scale height"
);
if
(
potential
->
x
_trunc
>
potential
->
x
_max
)
error
(
"Potential truncation
x
larger than maximal z"
);
if
(
potential
->
x
_trunc
<
potential
->
scale_height
)
error
(
"Potential truncation
x
smaller than scale height"
);
/* Compute derived quantities */
potential
->
scale_height_inv
=
1
.
/
potential
->
scale_height
;
potential
->
norm
=
2
.
*
M_PI
*
phys_const
->
const_newton_G
*
potential
->
surface_density
;
potential
->
norm_over_G
=
2
*
M_PI
*
potential
->
surface_density
;
potential
->
z
_trans
=
potential
->
z
_max
-
potential
->
z
_trunc
;
potential
->
x
_trans
=
potential
->
x
_max
-
potential
->
x
_trunc
;
if
(
potential
->
z
_trans
!=
0
.
f
)
potential
->
z
_trans_inv
=
1
.
/
potential
->
z
_trans
;
if
(
potential
->
x
_trans
!=
0
.
f
)
potential
->
x
_trans_inv
=
1
.
/
potential
->
x
_trans
;
else
potential
->
z
_trans_inv
=
FLT_MAX
;
potential
->
x
_trans_inv
=
FLT_MAX
;