From f5acadad7e70d14eca045c6d0b62086d460932af Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Wed, 25 Sep 2019 18:35:15 +0200
Subject: [PATCH] Re-introduced rho_dh in the Pressure-Energy SPH hydro case to
allow for the alternative neighbour number definition.
---
src/hydro/PressureEnergy/hydro.h | 4 ++++
src/hydro/PressureEnergy/hydro_iact.h | 4 ++++
src/hydro/PressureEnergy/hydro_part.h | 3 +++
3 files changed, 11 insertions(+)
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index 6ce6c5c526..0dcad24277 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -517,6 +517,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
+ p->density.rho_dh = 0.f;
p->pressure_bar = 0.f;
p->density.pressure_bar_dh = 0.f;
@@ -550,6 +551,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
p->pressure_bar += p->mass * p->u * kernel_root;
p->density.pressure_bar_dh -= hydro_dimension * p->mass * p->u * kernel_root;
p->density.wcount += kernel_root;
@@ -557,6 +559,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
+ p->density.rho_dh *= h_inv_dim_plus_one;
p->pressure_bar *= (h_inv_dim * hydro_gamma_minus_one);
p->density.pressure_bar_dh *= (h_inv_dim_plus_one * hydro_gamma_minus_one);
p->density.wcount *= h_inv_dim;
@@ -601,6 +604,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
p->pressure_bar =
p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
p->density.wcount = kernel_root * h_inv_dim;
+ p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f;
p->density.pressure_bar_dh = 0.f;
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index 62106bb35c..f0bd991522 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -69,6 +69,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
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 -=
@@ -82,6 +83,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
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);
@@ -147,6 +150,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
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;
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
index 6231ad4e13..6458d37e82 100644
--- a/src/hydro/PressureEnergy/hydro_part.h
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -133,6 +133,9 @@ struct part {
/*! Derivative of the neighbour number with respect to h. */
float wcount_dh;
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
/*! Derivative of the weighted pressure with respect to h */
float pressure_bar_dh;
--
GitLab