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
5f04d794
Commit
5f04d794
authored
Aug 16, 2016
by
Matthieu Schaller
Browse files
Created all the files for the Pressure-Entropy SPH routines.
parent
66d6d525
Changes
10
Hide whitespace changes
Inline
Side-by-side
src/const.h
View file @
5f04d794
...
...
@@ -60,7 +60,8 @@
/* SPH variant to use */
//#define MINIMAL_SPH
#define GADGET2_SPH
//#define GADGET2_SPH
#define HOPKINS_PE_SPH
//#define DEFAULT_SPH
/* Self gravity stuff. */
...
...
src/debug.c
View file @
5f04d794
...
...
@@ -41,6 +41,8 @@
#include
"./hydro/Minimal/hydro_debug.h"
#elif defined(GADGET2_SPH)
#include
"./hydro/Gadget2/hydro_debug.h"
#elif defined(HOPKINS_PE_SPH)
#include
"./hydro/PressureEntropy/hydro_debug.h"
#elif defined(DEFAULT_SPH)
#include
"./hydro/Default/hydro_debug.h"
#else
...
...
src/hydro.h
View file @
5f04d794
...
...
@@ -34,6 +34,10 @@
#include
"./hydro/Gadget2/hydro.h"
#include
"./hydro/Gadget2/hydro_iact.h"
#define SPH_IMPLEMENTATION "Gadget-2 version of SPH (Springel 2005)"
#elif defined(HOPKINS_PE_SPH)
#include
"./hydro/PressureEntropy/hydro.h"
#include
"./hydro/PressureEntropy/hydro_iact.h"
#define SPH_IMPLEMENTATION "Pressure-Entropy formulation of SPH (Hopkins 2013)"
#elif defined(DEFAULT_SPH)
#include
"./hydro/Default/hydro.h"
#include
"./hydro/Default/hydro_iact.h"
...
...
src/hydro/PressureEntropy/hydro.h
0 → 100644
View file @
5f04d794
/*******************************************************************************
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENTROPY_HYDRO_H
#define SWIFT_PRESSURE_ENTROPY_HYDRO_H
/**
* @file PressureEntropy/hydro_part.h
* @brief Pressure-Entropy implementation of SPH (Particle definition)
*
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
* Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
*/
#include
"adiabatic_index.h"
#include
"dimension.h"
#include
"equation_of_state.h"
#include
"hydro_properties.h"
#include
"kernel_hydro.h"
/**
* @brief Returns the internal energy of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_internal_energy
(
const
struct
part
*
restrict
p
,
float
dt
)
{
const
float
entropy
=
p
->
entropy
+
p
->
entropy_dt
*
dt
;
return
gas_internal_energy_from_entropy
(
p
->
rho
,
entropy
);
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_pressure
(
const
struct
part
*
restrict
p
,
float
dt
)
{
const
float
entropy
=
p
->
entropy
+
p
->
entropy_dt
*
dt
;
return
gas_pressure_from_entropy
(
p
->
rho
,
entropy
);
}
/**
* @brief Returns the entropy of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_entropy
(
const
struct
part
*
restrict
p
,
float
dt
)
{
return
p
->
entropy
+
p
->
entropy_dt
*
dt
;
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_get_soundspeed
(
const
struct
part
*
restrict
p
,
float
dt
)
{
return
p
->
force
.
soundspeed
;
}
/**
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
*
* This overrides the current state of the particle but does *not* change its
* time-derivatives
*
* @param p The particle
* @param u The new internal energy
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_set_internal_energy
(
struct
part
*
restrict
p
,
float
u
)
{
p
->
entropy
=
gas_entropy_from_internal_energy
(
p
->
rho
,
u
);
}
/**
* @brief Modifies the thermal state of a particle to the imposed entropy
*
* This overrides the current state of the particle but does *not* change its
* time-derivatives
*
* @param p The particle
* @param S The new entropy
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_set_entropy
(
struct
part
*
restrict
p
,
float
S
)
{
p
->
entropy
=
S
;
}
/**
* @brief Computes the hydro time-step of a given particle
*
* @param p Pointer to the particle data
* @param xp Pointer to the extended particle data
*
*/
__attribute__
((
always_inline
))
INLINE
static
float
hydro_compute_timestep
(
const
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
,
const
struct
hydro_props
*
restrict
hydro_properties
)
{
const
float
CFL_condition
=
hydro_properties
->
CFL_condition
;
/* CFL condition */
const
float
dt_cfl
=
2
.
f
*
kernel_gamma
*
CFL_condition
*
p
->
h
/
p
->
force
.
v_sig
;
return
dt_cfl
;
}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_first_init_part
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
)
{
p
->
ti_begin
=
0
;
p
->
ti_end
=
0
;
xp
->
v_full
[
0
]
=
p
->
v
[
0
];
xp
->
v_full
[
1
]
=
p
->
v
[
1
];
xp
->
v_full
[
2
]
=
p
->
v
[
2
];
}
/**
* @brief Prepares a particle for the density calculation.
*
* Zeroes all the relevant arrays in preparation for the sums taking place in
* the variaous density tasks
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_init_part
(
struct
part
*
restrict
p
)
{
p
->
density
.
wcount
=
0
.
f
;
p
->
density
.
wcount_dh
=
0
.
f
;
p
->
rho
=
0
.
f
;
p
->
rho_dh
=
0
.
f
;
p
->
density
.
div_v
=
0
.
f
;
p
->
density
.
rot_v
[
0
]
=
0
.
f
;
p
->
density
.
rot_v
[
1
]
=
0
.
f
;
p
->
density
.
rot_v
[
2
]
=
0
.
f
;
}
/**
* @brief Finishes the density calculation.
*
* Multiplies the density and number of neighbours by the appropiate constants
* and add the self-contribution term.
*
* @param p The particle to act upon
* @param ti_current The current time (on the integer timeline)
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_end_density
(
struct
part
*
restrict
p
,
int
ti_current
)
{
/* Some smoothing length multiples. */
const
float
h
=
p
->
h
;
const
float
h_inv
=
1
.
0
f
/
h
;
/* 1/h */
const
float
h_inv_dim
=
pow_dimension
(
h_inv
);
/* 1/h^d */
const
float
h_inv_dim_plus_one
=
h_inv_dim
*
h_inv
;
/* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p
->
rho
+=
p
->
mass
*
kernel_root
;
p
->
rho_dh
-=
hydro_dimension
*
p
->
mass
*
kernel_root
;
p
->
density
.
wcount
+=
kernel_root
;
/* Finish the calculation by inserting the missing h-factors */
p
->
rho
*=
h_inv_dim
;
p
->
rho_dh
*=
h_inv_dim_plus_one
;
p
->
density
.
wcount
*=
kernel_norm
;
p
->
density
.
wcount_dh
*=
h_inv
*
kernel_gamma
*
kernel_norm
;
const
float
irho
=
1
.
f
/
p
->
rho
;
/* Compute the derivative term */
p
->
rho_dh
=
1
.
f
/
(
1
.
f
+
hydro_dimension_inv
*
p
->
h
*
p
->
rho_dh
*
irho
);
/* Finish calculation of the velocity curl components */
p
->
density
.
rot_v
[
0
]
*=
h_inv_dim_plus_one
*
irho
;
p
->
density
.
rot_v
[
1
]
*=
h_inv_dim_plus_one
*
irho
;
p
->
density
.
rot_v
[
2
]
*=
h_inv_dim_plus_one
*
irho
;
/* Finish calculation of the velocity divergence */
p
->
density
.
div_v
*=
h_inv_dim_plus_one
*
irho
;
}
/**
* @brief Prepare a particle for the force calculation.
*
* Computes viscosity term, conduction term and smoothing length gradient terms.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param ti_current The current time (on the timeline)
* @param timeBase The minimal time-step size
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_prepare_force
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
,
int
ti_current
,
double
timeBase
)
{
const
float
fac_mu
=
1
.
f
;
/* Will change with cosmological integration */
/* Compute the norm of the curl */
const
float
curl_v
=
sqrtf
(
p
->
density
.
rot_v
[
0
]
*
p
->
density
.
rot_v
[
0
]
+
p
->
density
.
rot_v
[
1
]
*
p
->
density
.
rot_v
[
1
]
+
p
->
density
.
rot_v
[
2
]
*
p
->
density
.
rot_v
[
2
]);
/* Compute the norm of div v */
const
float
abs_div_v
=
fabsf
(
p
->
density
.
div_v
);
/* Compute the pressure */
const
float
half_dt
=
(
ti_current
-
(
p
->
ti_begin
+
p
->
ti_end
)
/
2
)
*
timeBase
;
const
float
pressure
=
hydro_get_pressure
(
p
,
half_dt
);
const
float
irho
=
1
.
f
/
p
->
rho
;
/* Divide the pressure by the density and density gradient */
const
float
P_over_rho2
=
pressure
*
irho
*
irho
*
p
->
rho_dh
;
/* Compute the sound speed */
const
float
soundspeed
=
sqrtf
(
hydro_gamma
*
pressure
*
irho
);
/* Compute the Balsara switch */
const
float
balsara
=
abs_div_v
/
(
abs_div_v
+
curl_v
+
0
.
0001
f
*
soundspeed
/
fac_mu
/
p
->
h
);
/* Update variables. */
p
->
force
.
P_over_rho2
=
P_over_rho2
;
p
->
force
.
soundspeed
=
soundspeed
;
p
->
force
.
balsara
=
balsara
;
}
/**
* @brief Reset acceleration fields of a particle
*
* Resets all hydro acceleration and time derivative fields in preparation
* for the sums taking place in the variaous force tasks
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_reset_acceleration
(
struct
part
*
restrict
p
)
{
/* Reset the acceleration. */
p
->
a_hydro
[
0
]
=
0
.
0
f
;
p
->
a_hydro
[
1
]
=
0
.
0
f
;
p
->
a_hydro
[
2
]
=
0
.
0
f
;
/* Reset the time derivatives. */
p
->
entropy_dt
=
0
.
0
f
;
p
->
force
.
h_dt
=
0
.
0
f
;
/* Reset maximal signal velocity */
p
->
force
.
v_sig
=
0
.
0
f
;
}
/**
* @brief Predict additional particle fields forward in time when drifting
*
* @param p The particle
* @param xp The extended data of the particle
* @param t0 The time at the start of the drift
* @param t1 The time at the end of the drift
* @param timeBase The minimal time-step size
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_predict_extra
(
struct
part
*
restrict
p
,
const
struct
xpart
*
restrict
xp
,
int
t0
,
int
t1
,
double
timeBase
)
{
/* Drift the pressure */
const
float
dt_entr
=
(
t1
-
(
p
->
ti_begin
+
p
->
ti_end
)
/
2
)
*
timeBase
;
const
float
pressure
=
hydro_get_pressure
(
p
,
dt_entr
);
const
float
irho
=
1
.
f
/
p
->
rho
;
/* Divide the pressure by the density and density gradient */
const
float
P_over_rho2
=
pressure
*
irho
*
irho
*
p
->
rho_dh
;
/* Compute the new sound speed */
const
float
soundspeed
=
sqrtf
(
hydro_gamma
*
pressure
*
irho
);
/* Update variables */
p
->
force
.
P_over_rho2
=
P_over_rho2
;
p
->
force
.
soundspeed
=
soundspeed
;
}
/**
* @brief Finishes the force calculation.
*
* Multiplies the forces and accelerationsby the appropiate constants
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_end_force
(
struct
part
*
restrict
p
)
{
p
->
force
.
h_dt
*=
p
->
h
*
hydro_dimension_inv
;
p
->
entropy_dt
*=
0
.
5
f
*
hydro_gamma_minus_one
*
pow_minus_gamma_minus_one
(
p
->
rho
);
}
/**
* @brief Kick the additional variables
*
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt The time-step for this kick
* @param half_dt The half time-step for this kick
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_kick_extra
(
struct
part
*
restrict
p
,
struct
xpart
*
restrict
xp
,
float
dt
,
float
half_dt
)
{
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
const
float
entropy_change
=
p
->
entropy_dt
*
dt
;
if
(
entropy_change
>
-
0
.
5
f
*
p
->
entropy
)
p
->
entropy
+=
entropy_change
;
else
p
->
entropy
*=
0
.
5
f
;
/* Do not 'overcool' when timestep increases */
if
(
p
->
entropy
+
p
->
entropy_dt
*
half_dt
<
0
.
5
f
*
p
->
entropy
)
p
->
entropy_dt
=
-
0
.
5
f
*
p
->
entropy
/
half_dt
;
}
/**
* @brief Converts hydro quantity of a particle at the start of a run
*
* Requires the density to be known
*
* @param p The particle to act upon
*/
__attribute__
((
always_inline
))
INLINE
static
void
hydro_convert_quantities
(
struct
part
*
restrict
p
)
{
/* We read u in the entropy field. We now get S from u */
p
->
entropy
=
gas_entropy_from_internal_energy
(
p
->
rho
,
p
->
entropy
);
}
#endif
/* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
src/hydro/PressureEntropy/hydro_debug.h
0 → 100644
View file @
5f04d794
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (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/>.
*
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H
#define SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H
/**
* @file PressureEntropy/hydro_debug.h
* @brief Pressure-Entropy implementation of SPH (Particle definition)
*
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
* Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
*/
__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],
\n
"
"h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
"P=%.3e, P_over_rho2=%.3e, S=%.3e, dS/dt=%.3e, c=%.3e
\n
"
"divV=%.3e, rotV=[%.3e,%.3e,%.3e], balsara=%.3e
\n
"
"v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d
\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
->
h
,
p
->
density
.
wcount
,
p
->
density
.
wcount_dh
,
p
->
mass
,
p
->
rho_dh
,
p
->
rho
,
hydro_get_pressure
(
p
,
0
.),
p
->
force
.
P_over_rho2
,
p
->
entropy
,
p
->
entropy_dt
,
p
->
force
.
soundspeed
,
p
->
density
.
div_v
,
p
->
density
.
rot_v
[
0
],
p
->
density
.
rot_v
[
1
],
p
->
density
.
rot_v
[
2
],
p
->
force
.
balsara
,
p
->
force
.
v_sig
,
p
->
force
.
h_dt
,
p
->
ti_begin
,
p
->
ti_end
);
}
#endif
/* SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H */
src/hydro/PressureEntropy/hydro_iact.h
0 → 100644
View file @
5f04d794
/*******************************************************************************
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENTROPY_HYDRO_IACT_H
#define SWIFT_PRESSURE_ENTROPY_HYDRO_IACT_H
/**
* @file PressureEntropy/hydro_part.h
* @brief Pressure-Entropy implementation of SPH (Particle definition)
*
* The thermal variable is the entropy (S) and the entropy is smoothed over
* contact discontinuities to prevent spurious surface tension.
*
* Follows Hopkins, P., MNRAS, 2013, Volume 428, Issue 4, pp. 2840-2856
*/
/**
* @brief Density loop
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_density
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
float
wi
,
wi_dx
;
float
wj
,
wj_dx
;
float
dv
[
3
],
curlvr
[
3
];
/* Get the masses. */
const
float
mi
=
pi
->
mass
;
const
float
mj
=
pj
->
mass
;
/* Get r and r inverse. */
const
float
r
=
sqrtf
(
r2
);
const
float
r_inv
=
1
.
0
f
/
r
;
/* Compute the kernel function for pi */
const
float
hi_inv
=
1
.
f
/
hi
;
const
float
ui
=
r
*
hi_inv
;
kernel_deval
(
ui
,
&
wi
,
&
wi_dx
);
/* Compute contribution to the density */
pi
->
rho
+=
mj
*
wi
;
pi
->
rho_dh
-=
mj
*
(
hydro_dimension
*
wi
+
ui
*
wi_dx
);
/* Compute contribution to the number of neighbours */
pi
->
density
.
wcount
+=
wi
;
pi
->
density
.
wcount_dh
-=
ui
*
wi_dx
;
/* Compute the kernel function for pj */
const
float
hj_inv
=
1
.
f
/
hj
;
const
float
uj
=
r
*
hj_inv
;
kernel_deval
(
uj
,
&
wj
,
&
wj_dx
);
/* Compute contribution to the density */
pj
->
rho
+=
mi
*
wj
;
pj
->
rho_dh
-=
mi
*
(
hydro_dimension
*
wj
+
uj
*
wj_dx
);
/* Compute contribution to the number of neighbours */
pj
->
density
.
wcount
+=
wj
;
pj
->
density
.
wcount_dh
-=
uj
*
wj_dx
;
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
];
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 loop (Vectorized version)
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_vec_density
(
float
*
R2
,
float
*
Dx
,
float
*
Hi
,
float
*
Hj
,
struct
part
**
pi
,
struct
part
**
pj
)
{
error
(
"Vectorized version of Pressure-Entropy SPH routine not existatn yet."
);
}
/**
* @brief Density loop (non-symmetric version)
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_density
(
float
r2
,
float
*
dx
,
float
hi
,
float
hj
,
struct
part
*
pi
,
struct
part
*
pj
)
{
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
ri
=
1
.
0
f
/
r
;
/* Compute the kernel function */
const
float
h_inv
=
1
.
0
f
/
hi
;
const
float
u
=
r
*
h_inv
;
kernel_deval
(
u
,
&
wi
,
&
wi_dx
);
/* Compute contribution to the density */
pi
->
rho
+=
mj
*
wi
;
pi
->
rho_dh
-=
mj
*
(
hydro_dimension
*
wi
+
u
*
wi_dx
);
/* Compute contribution to the number of neighbours */
pi
->
density
.
wcount
+=
wi
;
pi
->
density
.
wcount_dh
-=
u
*
wi_dx
;
const
float
fac
=
mj
*
wi_dx
*
ri
;
/* 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
];