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
4d69cea1
Commit
4d69cea1
authored
Nov 07, 2018
by
Josh Borrow
Committed by
Matthieu Schaller
Nov 07, 2018
Browse files
Add P-U with M&M Switch
parent
188f5537
Changes
23
Expand all
Hide whitespace changes
Inline
Side-by-side
configure.ac
View file @
4d69cea1
...
...
@@ -1150,7 +1150,7 @@ esac
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
[Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
[Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy,
pressure-energy-monaghan,
default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
...
...
@@ -1177,6 +1177,9 @@ case "$with_hydro" in
pressure-energy)
AC_DEFINE([HOPKINS_PU_SPH], [1], [Pressure-Energy SPH])
;;
pressure-energy-monaghan)
AC_DEFINE([HOPKINS_PU_SPH_MONAGHAN], [1], [Pressure-Energy SPH with M&M Variable A.V.])
;;
default)
AC_DEFINE([DEFAULT_SPH], [1], [Default SPH])
;;
...
...
examples/KelvinHelmholtz_2D/makeMovie.py
View file @
4d69cea1
...
...
@@ -91,6 +91,7 @@ if __name__ == "__main__":
# Creation of first frame
fig
,
ax
=
plt
.
subplots
(
1
,
1
,
figsize
=
(
1
,
1
),
frameon
=
False
)
ax
.
axis
(
"off"
)
# Remove annoying black frame.
data_x
,
data_y
,
density
=
load_and_extract
(
"kelvinhelmholtz_0000.hdf5"
)
...
...
examples/SodShock_1D/plotSolution.py
View file @
4d69cea1
...
...
@@ -70,11 +70,11 @@ snap = int(sys.argv[1])
sim
=
h5py
.
File
(
"sodShock_%04d.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"
]
scheme
=
str
(
sim
[
"/HydroScheme"
].
attrs
[
"Scheme"
]
)
kernel
=
str
(
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"
]
git
=
str
(
sim
[
"Code"
].
attrs
[
"Git Revision"
]
)
x
=
sim
[
"/PartType0/Coordinates"
][:,
0
]
v
=
sim
[
"/PartType0/Velocities"
][:,
0
]
...
...
@@ -82,6 +82,11 @@ u = sim["/PartType0/InternalEnergy"][:]
S
=
sim
[
"/PartType0/Entropy"
][:]
P
=
sim
[
"/PartType0/Pressure"
][:]
rho
=
sim
[
"/PartType0/Density"
][:]
try
:
alpha
=
sim
[
"/PartType0/Viscosity"
][:]
plot_alpha
=
True
except
:
plot_alpha
=
False
N
=
1000
# Number of points
x_min
=
-
1.
...
...
@@ -259,14 +264,23 @@ ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
xlim
(
-
0.5
,
0.5
)
ylim
(
0.8
,
2.2
)
# Entropy profile ---------------------------------
# Entropy
/alpha
profile ---------------------------------
subplot
(
235
)
plot
(
x
,
S
,
'.'
,
color
=
'r'
,
ms
=
4.0
)
plot
(
x_s
,
s_s
,
'--'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
if
plot_alpha
:
plot
(
x
,
alpha
,
'.'
,
color
=
'r'
,
ms
=
4.0
)
ylabel
(
r
"${\rm{Viscosity}}~\alpha$"
,
labelpad
=
0
)
# Show location of shock
plot
([
x_56
,
x_56
],
[
-
100
,
100
],
color
=
"k"
,
alpha
=
0.5
,
ls
=
"dashed"
,
lw
=
1.2
)
ylim
(
0
,
1
)
else
:
plot
(
x
,
S
,
'.'
,
color
=
'r'
,
ms
=
4.0
)
plot
(
x_s
,
s_s
,
'--'
,
color
=
'k'
,
alpha
=
0.8
,
lw
=
1.2
)
ylabel
(
"${
\\
rm{Entropy}}~S$"
,
labelpad
=
0
)
ylim
(
0.8
,
3.8
)
xlabel
(
"${
\\
rm{Position}}~x$"
,
labelpad
=
0
)
ylabel
(
"${
\\
rm{Entropy}}~S$"
,
labelpad
=
0
)
xlim
(
-
0.5
,
0.5
)
ylim
(
0.8
,
3.8
)
# Information -------------------------------------
subplot
(
236
,
frameon
=
False
)
...
...
@@ -284,5 +298,6 @@ ylim(0, 1)
xticks
([])
yticks
([])
tight_layout
()
savefig
(
"SodShock.png"
,
dpi
=
200
)
src/debug.c
View file @
4d69cea1
...
...
@@ -48,6 +48,8 @@
#include
"./hydro/PressureEntropy/hydro_debug.h"
#elif defined(HOPKINS_PU_SPH)
#include
"./hydro/PressureEnergy/hydro_debug.h"
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#include
"./hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h"
#elif defined(DEFAULT_SPH)
#include
"./hydro/Default/hydro_debug.h"
#elif defined(GIZMO_MFV_SPH)
...
...
src/hydro.h
View file @
4d69cea1
...
...
@@ -45,6 +45,12 @@
#include
"./hydro/PressureEnergy/hydro.h"
#include
"./hydro/PressureEnergy/hydro_iact.h"
#define SPH_IMPLEMENTATION "Pressure-Energy SPH (Hopkins 2013)"
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#include
"./hydro/PressureEnergyMorrisMonaghanAV/hydro.h"
#include
"./hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h"
#define SPH_IMPLEMENTATION \
"Pressure-Energy SPH (Hopkins 2013) with a Morris and Monaghan (1997) " \
"variable artificial viscosity."
#elif defined(DEFAULT_SPH)
#include
"./hydro/Default/hydro.h"
#include
"./hydro/Default/hydro_iact.h"
...
...
src/hydro/PressureEnergy/hydro.h
View file @
4d69cea1
...
...
@@ -472,8 +472,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p
->
density
.
rot_v
[
1
]
*=
h_inv_dim_plus_one
*
a_inv2
*
rho_inv
;
p
->
density
.
rot_v
[
2
]
*=
h_inv_dim_plus_one
*
a_inv2
*
rho_inv
;
/* Finish calculation of the velocity divergence */
p
->
density
.
div_v
*=
h_inv_dim_plus_one
*
rho_inv
*
a_inv2
;
/* Finish calculation of the velocity divergence, including hubble flow term
*/
p
->
density
.
div_v
*=
h_inv_dim_plus_one
*
rho_inv
*
a_inv2
+
cosmo
->
H
*
hydro_dimension
;
}
/**
...
...
src/hydro/PressureEnergy/hydro_iact.h
View file @
4d69cea1
...
...
@@ -284,7 +284,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
wj_dr
*
dvdr
*
r_inv
;
/* Viscosity term */
const
float
visc_du_term
=
0
.
5
f
*
visc_acc_term
*
dvdr
;
const
float
visc_du_term
=
0
.
5
f
*
visc_acc_term
*
dvdr
_Hubble
;
/* Assemble the energy equation term */
const
float
du_dt_i
=
sph_du_term_i
+
visc_du_term
;
...
...
@@ -404,7 +404,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
wi_dr
*
dvdr
*
r_inv
;
/* Viscosity term */
const
float
visc_du_term
=
0
.
5
f
*
visc_acc_term
*
dvdr
;
const
float
visc_du_term
=
0
.
5
f
*
visc_acc_term
*
dvdr
_Hubble
;
/* Assemble the energy equation term */
const
float
du_dt_i
=
sph_du_term_i
+
visc_du_term
;
...
...
src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h
0 → 100644
View file @
4d69cea1
This diff is collapsed.
Click to expand it.
src/hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h
0 → 100644
View file @
4d69cea1
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
* Josh Borrow (joshua.borrow@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/>.
*
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H
#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H
/**
* @file PressureEnergy/hydro_debug.h
* @brief P-U conservative implementation of SPH (Debugging routines)
*
* The thermal variable is the internal energy (u). A simple variable
* viscosity term (Morris & Monaghan 1997) with a Balsara switch is
* implemented.
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_debug_particle
(
const
struct
part
*
p
,
const
struct
xpart
*
xp
)
{
printf
(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e]
\n
a=[%.3e,%.3e,%.3e], "
"u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e
\n
"
"h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e,
\n
"
"p_dh=%.3e, p_bar=%.3e
\n
"
"time_bin=%d, alpha=%.3e
\n
"
,
p
->
x
[
0
],
p
->
x
[
1
],
p
->
x
[
2
],
p
->
v
[
0
],
p
->
v
[
1
],
p
->
v
[
2
],
xp
->
v_full
[
0
],
xp
->
v_full
[
1
],
xp
->
v_full
[
2
],
p
->
a_hydro
[
0
],
p
->
a_hydro
[
1
],
p
->
a_hydro
[
2
],
p
->
u
,
p
->
u_dt
,
p
->
force
.
v_sig
,
hydro_get_comoving_pressure
(
p
),
p
->
h
,
p
->
force
.
h_dt
,
(
int
)
p
->
density
.
wcount
,
p
->
mass
,
p
->
density
.
rho_dh
,
p
->
rho
,
p
->
density
.
pressure_bar_dh
,
p
->
pressure_bar
,
p
->
time_bin
,
p
->
alpha
);
}
#endif
/* SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_DEBUG_H */
src/hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h
0 → 100644
View file @
4d69cea1
/*******************************************************************************
* This file is part* of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) &
* Josh Borrow (joshua.borrow@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/>.
*
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IACT_H
#define SWIFT_PRESSURE_ENERGY_MORRIS_HYDRO_IACT_H
/**
* @file PressureEnergy/hydro_iact.h
* @brief P-U implementation of SPH (Neighbour loop equations)
*
* The thermal variable is the internal energy (u). A simple variable
* viscosity term (Morris & Monaghan 1997) with a Balsara switch is
* implemented.
*
* No thermal conduction term is implemented.
*
* See PressureEnergy/hydro.h for references.
*/
#include
"adiabatic_index.h"
#include
"minmax.h"
/**
* @brief Density interaction between two part*icles.
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_density
(
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
,
float
a
,
float
H
)
{
float
wi
,
wj
,
wi_dx
,
wj_dx
;
float
dv
[
3
],
curlvr
[
3
];
const
float
r
=
sqrtf
(
r2
);
/* Get the masses. */
const
float
mi
=
pi
->
mass
;
const
float
mj
=
pj
->
mass
;
/* Compute density of pi. */
const
float
hi_inv
=
1
.
f
/
hi
;
const
float
ui
=
r
*
hi_inv
;
kernel_deval
(
ui
,
&
wi
,
&
wi_dx
);
pi
->
rho
+=
mj
*
wi
;
pi
->
density
.
rho_dh
-=
mj
*
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
pi
->
pressure_bar
+=
mj
*
wi
*
pj
->
u
;
pi
->
density
.
pressure_bar_dh
-=
mj
*
pj
->
u
*
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
pi
->
density
.
wcount
+=
wi
;
pi
->
density
.
wcount_dh
-=
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
/* Compute density of pj. */
const
float
hj_inv
=
1
.
f
/
hj
;
const
float
uj
=
r
*
hj_inv
;
kernel_deval
(
uj
,
&
wj
,
&
wj_dx
);
pj
->
rho
+=
mi
*
wj
;
pj
->
density
.
rho_dh
-=
mi
*
(
hydro_dimension
*
wj
+
uj
*
wj_dx
);
pj
->
pressure_bar
+=
mi
*
wj
*
pi
->
u
;
pj
->
density
.
pressure_bar_dh
-=
mi
*
pi
->
u
*
(
hydro_dimension
*
wj
+
uj
*
wj_dx
);
pj
->
density
.
wcount
+=
wj
;
pj
->
density
.
wcount_dh
-=
(
hydro_dimension
*
wj
+
uj
*
wj_dx
);
/* Now we need to compute the div terms */
const
float
r_inv
=
1
.
f
/
r
;
const
float
faci
=
mj
*
wi_dx
*
r_inv
;
const
float
facj
=
mi
*
wj_dx
*
r_inv
;
/* Compute dv dot r */
dv
[
0
]
=
pi
->
v
[
0
]
-
pj
->
v
[
0
];
dv
[
1
]
=
pi
->
v
[
1
]
-
pj
->
v
[
1
];
dv
[
2
]
=
pi
->
v
[
2
]
-
pj
->
v
[
2
];
const
float
dvdr
=
dv
[
0
]
*
dx
[
0
]
+
dv
[
1
]
*
dx
[
1
]
+
dv
[
2
]
*
dx
[
2
];
pi
->
density
.
div_v
-=
faci
*
dvdr
;
pj
->
density
.
div_v
-=
facj
*
dvdr
;
/* Compute dv cross r */
curlvr
[
0
]
=
dv
[
1
]
*
dx
[
2
]
-
dv
[
2
]
*
dx
[
1
];
curlvr
[
1
]
=
dv
[
2
]
*
dx
[
0
]
-
dv
[
0
]
*
dx
[
2
];
curlvr
[
2
]
=
dv
[
0
]
*
dx
[
1
]
-
dv
[
1
]
*
dx
[
0
];
pi
->
density
.
rot_v
[
0
]
+=
faci
*
curlvr
[
0
];
pi
->
density
.
rot_v
[
1
]
+=
faci
*
curlvr
[
1
];
pi
->
density
.
rot_v
[
2
]
+=
faci
*
curlvr
[
2
];
/* Negative because of the change in sign of dx & dv. */
pj
->
density
.
rot_v
[
0
]
+=
facj
*
curlvr
[
0
];
pj
->
density
.
rot_v
[
1
]
+=
facj
*
curlvr
[
1
];
pj
->
density
.
rot_v
[
2
]
+=
facj
*
curlvr
[
2
];
}
/**
* @brief Density interaction between two part*icles (non-symmetric).
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_density
(
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
const
struct
part
*
pj
,
float
a
,
float
H
)
{
float
wi
,
wi_dx
;
float
dv
[
3
],
curlvr
[
3
];
/* Get the masses. */
const
float
mj
=
pj
->
mass
;
/* Get r and r inverse. */
const
float
r
=
sqrtf
(
r2
);
const
float
h_inv
=
1
.
f
/
hi
;
const
float
ui
=
r
*
h_inv
;
kernel_deval
(
ui
,
&
wi
,
&
wi_dx
);
pi
->
rho
+=
mj
*
wi
;
pi
->
density
.
rho_dh
-=
mj
*
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
pi
->
pressure_bar
+=
mj
*
wi
*
pj
->
u
;
pi
->
density
.
pressure_bar_dh
-=
mj
*
pj
->
u
*
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
pi
->
density
.
wcount
+=
wi
;
pi
->
density
.
wcount_dh
-=
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
const
float
r_inv
=
1
.
f
/
r
;
const
float
faci
=
mj
*
wi_dx
*
r_inv
;
/* Compute dv dot r */
dv
[
0
]
=
pi
->
v
[
0
]
-
pj
->
v
[
0
];
dv
[
1
]
=
pi
->
v
[
1
]
-
pj
->
v
[
1
];
dv
[
2
]
=
pi
->
v
[
2
]
-
pj
->
v
[
2
];
const
float
dvdr
=
dv
[
0
]
*
dx
[
0
]
+
dv
[
1
]
*
dx
[
1
]
+
dv
[
2
]
*
dx
[
2
];
pi
->
density
.
div_v
-=
faci
*
dvdr
;
/* Compute dv cross r */
curlvr
[
0
]
=
dv
[
1
]
*
dx
[
2
]
-
dv
[
2
]
*
dx
[
1
];
curlvr
[
1
]
=
dv
[
2
]
*
dx
[
0
]
-
dv
[
0
]
*
dx
[
2
];
curlvr
[
2
]
=
dv
[
0
]
*
dx
[
1
]
-
dv
[
1
]
*
dx
[
0
];
pi
->
density
.
rot_v
[
0
]
+=
faci
*
curlvr
[
0
];
pi
->
density
.
rot_v
[
1
]
+=
faci
*
curlvr
[
1
];
pi
->
density
.
rot_v
[
2
]
+=
faci
*
curlvr
[
2
];
}
/**
* @brief Force interaction between two part*icles.
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_force
(
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
,
float
a
,
float
H
)
{
/* Cosmological factors entering the EoMs */
const
float
fac_mu
=
pow_three_gamma_minus_five_over_two
(
a
);
const
float
a2_Hubble
=
a
*
a
*
H
;
const
float
r
=
sqrtf
(
r2
);
const
float
r_inv
=
1
.
0
f
/
r
;
/* Recover some data */
const
float
mj
=
pj
->
mass
;
const
float
mi
=
pi
->
mass
;
const
float
miui
=
mi
*
pi
->
u
;
const
float
mjuj
=
mj
*
pj
->
u
;
const
float
rhoi
=
pi
->
rho
;
const
float
rhoj
=
pj
->
rho
;
/* Compute gradient terms */
const
float
f_ij
=
1
.
f
-
(
pi
->
force
.
f
/
mjuj
);
const
float
f_ji
=
1
.
f
-
(
pj
->
force
.
f
/
miui
);
/* Get the kernel for hi. */
const
float
hi_inv
=
1
.
0
f
/
hi
;
const
float
hid_inv
=
pow_dimension_plus_one
(
hi_inv
);
/* 1/h^(d+1) */
const
float
xi
=
r
*
hi_inv
;
float
wi
,
wi_dx
;
kernel_deval
(
xi
,
&
wi
,
&
wi_dx
);
const
float
wi_dr
=
hid_inv
*
wi_dx
;
/* Get the kernel for hj. */
const
float
hj_inv
=
1
.
0
f
/
hj
;
const
float
hjd_inv
=
pow_dimension_plus_one
(
hj_inv
);
/* 1/h^(d+1) */
const
float
xj
=
r
*
hj_inv
;
float
wj
,
wj_dx
;
kernel_deval
(
xj
,
&
wj
,
&
wj_dx
);
const
float
wj_dr
=
hjd_inv
*
wj_dx
;
/* Compute dv dot r. */
const
float
dvdr
=
(
pi
->
v
[
0
]
-
pj
->
v
[
0
])
*
dx
[
0
]
+
(
pi
->
v
[
1
]
-
pj
->
v
[
1
])
*
dx
[
1
]
+
(
pi
->
v
[
2
]
-
pj
->
v
[
2
])
*
dx
[
2
];
/* Includes the hubble flow term; not used for du/dt */
const
float
dvdr_Hubble
=
dvdr
+
a2_Hubble
*
r2
;
/* Are the part*icles moving towards each others ? */
const
float
omega_ij
=
min
(
dvdr_Hubble
,
0
.
f
);
const
float
mu_ij
=
fac_mu
*
r_inv
*
omega_ij
;
/* This is 0 or negative */
/* Compute sound speeds and signal velocity */
const
float
ci
=
pi
->
force
.
soundspeed
;
const
float
cj
=
pj
->
force
.
soundspeed
;
const
float
v_sig
=
ci
+
cj
-
3
.
f
*
mu_ij
;
/* Balsara term */
const
float
balsara_i
=
pi
->
force
.
balsara
;
const
float
balsara_j
=
pj
->
force
.
balsara
;
/* Construct the full viscosity term */
const
float
rho_ij
=
0
.
5
f
*
(
rhoi
+
rhoj
);
const
float
alpha
=
0
.
5
f
*
(
pi
->
alpha
+
pj
->
alpha
);
const
float
visc
=
-
0
.
25
f
*
alpha
*
v_sig
*
mu_ij
*
(
balsara_i
+
balsara_j
)
/
rho_ij
;
/* Convolve with the kernel */
const
float
visc_acc_term
=
0
.
5
f
*
visc
*
(
wi_dr
+
wj_dr
)
*
r_inv
;
/* SPH acceleration term */
const
float
sph_acc_term
=
pj
->
u
*
pi
->
u
*
hydro_gamma_minus_one
*
hydro_gamma_minus_one
*
((
f_ij
/
pi
->
pressure_bar
)
*
wi_dr
+
(
f_ji
/
pj
->
pressure_bar
)
*
wj_dr
)
*
r_inv
;
/* Assemble the acceleration */
const
float
acc
=
sph_acc_term
+
visc_acc_term
;
/* Use the force Luke ! */
pi
->
a_hydro
[
0
]
-=
mj
*
acc
*
dx
[
0
];
pi
->
a_hydro
[
1
]
-=
mj
*
acc
*
dx
[
1
];
pi
->
a_hydro
[
2
]
-=
mj
*
acc
*
dx
[
2
];
pj
->
a_hydro
[
0
]
+=
mi
*
acc
*
dx
[
0
];
pj
->
a_hydro
[
1
]
+=
mi
*
acc
*
dx
[
1
];
pj
->
a_hydro
[
2
]
+=
mi
*
acc
*
dx
[
2
];
/* Get the time derivative for u. */
const
float
sph_du_term_i
=
hydro_gamma_minus_one
*
hydro_gamma_minus_one
*
pj
->
u
*
pi
->
u
*
(
f_ij
/
pi
->
pressure_bar
)
*
wi_dr
*
dvdr
*
r_inv
;
const
float
sph_du_term_j
=
hydro_gamma_minus_one
*
hydro_gamma_minus_one
*
pi
->
u
*
pj
->
u
*
(
f_ji
/
pj
->
pressure_bar
)
*
wj_dr
*
dvdr
*
r_inv
;
/* Viscosity term */
const
float
visc_du_term
=
0
.
5
f
*
visc_acc_term
*
dvdr_Hubble
;
/* Assemble the energy equation term */
const
float
du_dt_i
=
sph_du_term_i
+
visc_du_term
;
const
float
du_dt_j
=
sph_du_term_j
+
visc_du_term
;
/* Internal energy time derivative */
pi
->
u_dt
+=
du_dt_i
*
mj
;
pj
->
u_dt
+=
du_dt_j
*
mi
;
/* Get the time derivative for h. */
pi
->
force
.
h_dt
-=
mj
*
dvdr
*
r_inv
/
rhoj
*
wi_dr
;
pj
->
force
.
h_dt
-=
mi
*
dvdr
*
r_inv
/
rhoi
*
wj_dr
;
/* Update the signal velocity. */
pi
->
force
.
v_sig
=
max
(
pi
->
force
.
v_sig
,
v_sig
);
pj
->
force
.
v_sig
=
max
(
pj
->
force
.
v_sig
,
v_sig
);
}
/**
* @brief Force interaction between two part*icles (non-symmetric).
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_force
(
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
const
struct
part
*
pj
,
float
a
,
float
H
)
{
/* Cosmological factors entering the EoMs */
const
float
fac_mu
=
pow_three_gamma_minus_five_over_two
(
a
);
const
float
a2_Hubble
=
a
*
a
*
H
;
const
float
r
=
sqrtf
(
r2
);
const
float
r_inv
=
1
.
0
f
/
r
;
/* Recover some data */
// const float mi = pi->mass;
const
float
mj
=
pj
->
mass
;
const
float
mi
=
pi
->
mass
;
const
float
miui
=
mi
*
pi
->
u
;
const
float
mjuj
=
mj
*
pj
->
u
;
const
float
rhoi
=
pi
->
rho
;
const
float
rhoj
=
pj
->
rho
;
/* Compute gradient terms */
const
float
f_ij
=
1
.
f
-
(
pi
->
force
.
f
/
mjuj
);
const
float
f_ji
=
1
.
f
-
(
pj
->
force
.
f
/
miui
);
/* Get the kernel for hi. */
const
float
hi_inv
=
1
.
0
f
/
hi
;
const
float
hid_inv
=
pow_dimension_plus_one
(
hi_inv
);
/* 1/h^(d+1) */
const
float
xi
=
r
*
hi_inv
;
float
wi
,
wi_dx
;
kernel_deval
(
xi
,
&
wi
,
&
wi_dx
);
const
float
wi_dr
=
hid_inv
*
wi_dx
;
/* Get the kernel for hj. */
const
float
hj_inv
=
1
.
0
f
/
hj
;
const
float
hjd_inv
=
pow_dimension_plus_one
(
hj_inv
);
/* 1/h^(d+1) */
const
float
xj
=
r
*
hj_inv
;
float
wj
,
wj_dx
;
kernel_deval
(
xj
,
&
wj
,
&
wj_dx
);
const
float
wj_dr
=
hjd_inv
*
wj_dx
;
/* Compute dv dot r. */
const
float
dvdr
=
(
pi
->
v
[
0
]
-
pj
->
v
[
0
])
*
dx
[
0
]
+
(
pi
->
v
[
1
]
-
pj
->
v
[
1
])
*
dx
[
1
]
+
(
pi
->
v
[
2
]
-
pj
->
v
[
2
])
*
dx
[
2
];
/* Includes the hubble flow term; not used for du/dt */
const
float
dvdr_Hubble
=
dvdr
+
a2_Hubble
*
r2
;
/* Are the part*icles moving towards each others ? */
const
float
omega_ij
=
min
(
dvdr_Hubble
,
0
.
f
);
const
float
mu_ij
=
fac_mu
*
r_inv
*
omega_ij
;
/* This is 0 or negative */
/* Compute sound speeds and signal velocity */
const
float
ci
=
pi
->
force
.
soundspeed
;
const
float
cj
=
pj
->
force
.
soundspeed
;
const
float
v_sig
=
ci
+
cj
-
3
.
f
*
mu_ij
;
/* Balsara term */
const
float
balsara_i
=
pi
->
force
.
balsara
;
const
float
balsara_j
=
pj
->
force
.
balsara
;
/* Construct the full viscosity term */
const
float
rho_ij
=
0
.
5
f
*
(
rhoi
+
rhoj
);
const
float
alpha
=
0
.
5
f
*
(
pi
->
alpha
+
pj
->
alpha
);
const
float
visc
=
-
0
.
25
f
*
alpha
*
v_sig
*
mu_ij
*
(
balsara_i
+
balsara_j
)
/
rho_ij
;