diff --git a/src/cosmology.c b/src/cosmology.c
index c94fc3c63b2afc8005117648a19f5d5d8dba0bbb..5d5b45c2d00cd614e92da3de3335e7b558a727ea 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -457,7 +457,7 @@ void cosmology_init_no_cosmo(const struct swift_params *params,
c->log_a_begin = 0.;
c->log_a_end = 0.;
- c->H = 1.;
+ c->H = 0.;
c->a = 1.;
c->a_inv = 1.;
c->a2_inv = 1.;
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index e27d63dd6bcf3cfead622613537b778e97f841f7..e5396988333d77801122e2508f9c2ae6c1ab69fc 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -36,10 +36,20 @@
#include "minmax.h"
/**
- * @brief Density loop
+ * @brief Density interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_density(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ float r2, float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
float wi, wi_dx;
float wj, wj_dx;
@@ -117,10 +127,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
}
/**
- * @brief Density loop (non-symmetric version)
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__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 r2, float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
float wi, wi_dx;
float dv[3], curlvr[3];
@@ -399,10 +419,20 @@ runner_iact_nonsym_2_vec_density(float *R2, float *Dx, float *Dy, float *Dz,
#endif
/**
- * @brief Force loop
+ * @brief Force interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_force(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ float r2, float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
float wi, wj, wi_dx, wj_dx;
@@ -508,10 +538,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
}
/**
- * @brief Force loop (non-symmetric version)
+ * @brief Force interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
- float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ float r2, float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
float wi, wj, wi_dx, wj_dx;
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 07d24c8df984634e19cb1aaf8bc072158e1b4c95..4303fbc235c55ac5a9ab977011c1dcc198d0d59d 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -138,6 +138,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
@@ -149,6 +150,10 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
struct part *restrict parts_i = ci->parts;
struct part *restrict parts_j = cj->parts;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
for (int k = 0; k < 3; k++) {
@@ -196,7 +201,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2 && pi_active) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
}
if (r2 < hjg2 && pj_active) {
@@ -204,7 +209,7 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
dx[1] = -dx[1];
dx[2] = -dx[2];
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over the parts in cj. */
@@ -226,6 +231,7 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
@@ -237,6 +243,10 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
struct part *restrict parts_i = ci->parts;
struct part *restrict parts_j = cj->parts;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
for (int k = 0; k < 3; k++) {
@@ -286,17 +296,17 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
if (pi_active && pj_active) {
- IACT(r2, dx, hi, hj, pi, pj);
+ IACT(r2, dx, hi, hj, pi, pj, a, H);
} else if (pi_active) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
} else if (pj_active) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
}
} /* loop over the parts in cj. */
@@ -316,12 +326,17 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_hydro(c, e)) return;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
const int count = c->count;
struct part *restrict parts = c->parts;
@@ -365,17 +380,17 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
/* Hit or miss? */
if (doi && doj) {
- IACT(r2, dx, hi, hj, pi, pj);
+ IACT(r2, dx, hi, hj, pi, pj, a, H);
} else if (doi) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
} else if (doj) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -394,12 +409,17 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_hydro(c, e)) return;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
const int count = c->count;
struct part *restrict parts = c->parts;
@@ -443,17 +463,17 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
/* Hit or miss? */
if (doi && doj) {
- IACT(r2, dx, hi, hj, pi, pj);
+ IACT(r2, dx, hi, hj, pi, pj, a, H);
} else if (doi) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
} else if (doj) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -480,11 +500,18 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
int count, struct cell *restrict cj,
const double *shift) {
+ const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
+
TIMER_TIC;
const int count_j = cj->count;
struct part *restrict parts_j = cj->parts;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
/* Loop over the parts_i. */
for (int pid = 0; pid < count; pid++) {
@@ -496,7 +523,6 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
const float hig2 = hi * hi * kernel_gamma2;
#ifdef SWIFT_DEBUG_CHECKS
- const struct engine *e = r->e;
if (!part_is_active(pi, e))
error("Trying to correct smoothing length of inactive particle !");
#endif
@@ -526,7 +552,7 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2) {
- IACT_NONSYM(r2, dx, hi, pj->h, pi, pj);
+ IACT_NONSYM(r2, dx, hi, pj->h, pi, pj, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -553,11 +579,18 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
struct cell *restrict cj, const int sid, const int flipped,
const double *shift) {
+ const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
+
TIMER_TIC;
const int count_j = cj->count;
struct part *restrict parts_j = cj->parts;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
/* Pick-out the sorted lists. */
const struct entry *restrict sort_j = cj->sort[sid];
const float dxj = cj->dx_max_sort;
@@ -593,7 +626,6 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
- const struct engine *e = r->e;
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
@@ -604,7 +636,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -641,7 +673,6 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
- const struct engine *e = r->e;
/* Check that particles have been drifted to the current time */
if (pi->ti_drift != e->ti_current)
error("Particle pi not drifted to current time");
@@ -652,7 +683,7 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 < hig2) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -732,12 +763,15 @@ void DOPAIR_SUBSET_BRANCH(struct runner *r, struct cell *restrict ci,
void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
struct part *restrict parts, int *restrict ind, int count) {
-#ifdef SWIFT_DEBUG_CHECKS
const struct engine *e = r->e;
-#endif
+ const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
const int count_i = ci->count;
struct part *restrict parts_j = ci->parts;
@@ -778,7 +812,7 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
/* Hit or miss? */
if (r2 > 0.f && r2 < hig2) {
- IACT_NONSYM(r2, dx, hi, pj->h, pi, pj);
+ IACT_NONSYM(r2, dx, hi, pj->h, pi, pj, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -820,6 +854,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const double *shift) {
const struct engine *restrict e = r->e;
+ const struct cosmology *restrict cosmo = e->cosmology;
TIMER_TIC;
@@ -852,6 +887,10 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const double dj_min = sort_j[0].d;
const float dx_max = (ci->dx_max_sort + cj->dx_max_sort);
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
if (cell_is_active_hydro(ci, e)) {
/* Loop over the parts in ci. */
@@ -926,7 +965,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
/* Hit or miss? */
if (r2 < hig2) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
}
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -1006,7 +1045,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
/* Hit or miss? */
if (r2 < hjg2) {
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over the parts in ci. */
} /* loop over the parts in cj. */
@@ -1109,7 +1148,8 @@ void DOPAIR1_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const double *shift) {
- struct engine *restrict e = r->e;
+ const struct engine *restrict e = r->e;
+ const struct cosmology *restrict cosmo = e->cosmology;
TIMER_TIC;
@@ -1139,6 +1179,10 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
struct part *restrict parts_i = ci->parts;
struct part *restrict parts_j = cj->parts;
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
/* Maximal displacement since last rebuild */
const double dx_max = (ci->dx_max_sort + cj->dx_max_sort);
@@ -1269,7 +1313,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
/* Hit or miss?
(note that we will do the other condition in the reverse loop) */
if (r2 < hig2) {
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over the active parts in cj. */
}
@@ -1331,9 +1375,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
/* Does pj need to be updated too? */
if (part_is_active(pj, e))
- IACT(r2, dx, hi, hj, pi, pj);
+ IACT(r2, dx, hi, hj, pi, pj, a, H);
else
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
}
} /* loop over the parts in cj. */
} /* Is pi active? */
@@ -1419,7 +1463,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
/* Hit or miss?
(note that we must avoid the r2 < hig2 cases we already processed) */
if (r2 < hjg2 && r2 >= hig2) {
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
}
} /* loop over the active parts in ci. */
}
@@ -1484,9 +1528,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
/* Does pi need to be updated too? */
if (part_is_active(pi, e))
- IACT(r2, dx, hj, hi, pj, pi);
+ IACT(r2, dx, hj, hi, pj, pi, a, H);
else
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over the parts in ci. */
} /* Is pj active? */
@@ -1592,6 +1636,7 @@ void DOPAIR2_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
void DOSELF1(struct runner *r, struct cell *restrict c) {
const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
@@ -1610,6 +1655,10 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
countdt += 1;
}
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
/* Loop over the particles in the cell. */
for (int pid = 0; pid < count; pid++) {
@@ -1651,7 +1700,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
/* Hit or miss? */
if (r2 < hj * hj * kernel_gamma2) {
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over all other particles. */
}
@@ -1692,14 +1741,14 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
/* Which parts need to be updated? */
if (r2 < hig2 && doj)
- IACT(r2, dx, hi, hj, pi, pj);
+ IACT(r2, dx, hi, hj, pi, pj, a, H);
else if (!doj)
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
else {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
}
} /* loop over all other particles. */
@@ -1752,6 +1801,7 @@ void DOSELF1_BRANCH(struct runner *r, struct cell *c) {
void DOSELF2(struct runner *r, struct cell *restrict c) {
const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
TIMER_TIC;
@@ -1770,6 +1820,10 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
countdt += 1;
}
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
/* Loop over the particles in the cell. */
for (int pid = 0; pid < count; pid++) {
@@ -1811,7 +1865,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
/* Hit or miss? */
if (r2 < hig2 || r2 < hj * hj * kernel_gamma2) {
- IACT_NONSYM(r2, dx, hj, hi, pj, pi);
+ IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
}
} /* loop over all other particles. */
}
@@ -1850,9 +1904,9 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
/* Does pj need to be updated too? */
if (part_is_active(pj, e))
- IACT(r2, dx, hi, hj, pi, pj);
+ IACT(r2, dx, hi, hj, pi, pj, a, H);
else
- IACT_NONSYM(r2, dx, hi, hj, pi, pj);
+ IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
}
} /* loop over all other particles. */
}
diff --git a/src/tools.c b/src/tools.c
index 8162f22774615e6b481f4a113b986be50cdbf7aa..682ebd0ad7cc0514cb5904dbd112b2b1bc1da05d 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -74,6 +74,7 @@ void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic) {
// double maxratio = 1.0;
double r2, dx[3], rho = 0.0;
double rho_max = 0.0, rho_min = 100;
+ float a = 1.f, H = 0.f;
/* Loop over all particle pairs. */
for (j = 0; j < N; j++) {
@@ -94,7 +95,7 @@ void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic) {
r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
if (r2 < parts[j].h * parts[j].h || r2 < parts[k].h * parts[k].h) {
runner_iact_density(r2, NULL, parts[j].h, parts[k].h, &parts[j],
- &parts[k]);
+ &parts[k], a, H);
/* if ( parts[j].h / parts[k].h > maxratio )
{
maxratio = parts[j].h / parts[k].h;
@@ -133,12 +134,10 @@ void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic) {
void pairs_single_density(double *dim, long long int pid,
struct part *restrict parts, int N, int periodic) {
int i, k;
- // int mj, mk;
- // double maxratio = 1.0;
double r2, dx[3];
float fdx[3];
struct part p;
- // double ih = 12.0/6.25;
+ float a = 1.f, H = 0.f;
/* Find "our" part. */
for (k = 0; k < N && parts[k].id != pid; k++)
@@ -164,7 +163,7 @@ void pairs_single_density(double *dim, long long int pid,
}
r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
if (r2 < p.h * p.h) {
- runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
+ runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k], a, H);
/* printf( "pairs_simple: interacting particles %lli [%i,%i,%i] and %lli
[%i,%i,%i], r=%e.\n" ,
pid , (int)(p.x[0]*ih) , (int)(p.x[1]*ih) , (int)(p.x[2]*ih) ,
@@ -186,6 +185,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *pi, *pj;
const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
const struct engine *e = r->e;
+ float a = 1.f, H = 0.f;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
@@ -213,7 +213,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
if (r2 < hig2) {
/* Interact */
- runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj);
+ runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj, a, H);
}
}
}
@@ -244,7 +244,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
if (r2 < hjg2) {
/* Interact */
- runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi);
+ runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi, a, H);
}
}
}
@@ -256,6 +256,7 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *pi, *pj;
const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
const struct engine *e = r->e;
+ float a = 1.f, H = 0.f;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
@@ -285,7 +286,7 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
if (r2 < hig2 || r2 < hjg2) {
/* Interact */
- runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj);
+ runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj, a, H);
}
}
}
@@ -318,7 +319,7 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
if (r2 < hjg2 || r2 < hig2) {
/* Interact */
- runner_iact_nonsym_force(r2, dx, hj, pi->h, pj, pi);
+ runner_iact_nonsym_force(r2, dx, hj, pi->h, pj, pi, a, H);
}
}
}
@@ -328,6 +329,7 @@ void self_all_density(struct runner *r, struct cell *ci) {
float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
struct part *pi, *pj;
const struct engine *e = r->e;
+ float a = 1.f, H = 0.f;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
@@ -355,7 +357,7 @@ void self_all_density(struct runner *r, struct cell *ci) {
if (r2 < hig2 && part_is_active(pi, e)) {
/* Interact */
- runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj);
+ runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj, a, H);
}
/* Hit or miss? */
@@ -366,7 +368,7 @@ void self_all_density(struct runner *r, struct cell *ci) {
dxi[2] = -dxi[2];
/* Interact */
- runner_iact_nonsym_density(r2, dxi, hj, hi, pj, pi);
+ runner_iact_nonsym_density(r2, dxi, hj, hi, pj, pi, a, H);
}
}
}
@@ -375,6 +377,7 @@ void self_all_density(struct runner *r, struct cell *ci) {
void self_all_force(struct runner *r, struct cell *ci) {
float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
struct part *pi, *pj;
+ float a = 1.f, H = 0.f;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
@@ -402,7 +405,7 @@ void self_all_force(struct runner *r, struct cell *ci) {
if (r2 < hig2 || r2 < hjg2) {
/* Interact */
- runner_iact_force(r2, dxi, hi, hj, pi, pj);
+ runner_iact_force(r2, dxi, hi, hj, pi, pj, a, H);
}
}
}
@@ -417,6 +420,7 @@ void engine_single_density(double *dim, long long int pid,
double r2, dx[3];
float fdx[3];
struct part p;
+ float a = 1.f, H = 0.f;
/* Find "our" part. */
int k;
@@ -443,7 +447,7 @@ void engine_single_density(double *dim, long long int pid,
}
r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
if (r2 < p.h * p.h * kernel_gamma2) {
- runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
+ runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k], a, H);
}
}
@@ -460,6 +464,7 @@ void engine_single_force(double *dim, long long int pid,
double r2, dx[3];
float fdx[3];
struct part p;
+ float a = 1.f, H = 0.f;
/* Find "our" part. */
for (k = 0; k < N && parts[k].id != pid; k++)
@@ -488,7 +493,7 @@ void engine_single_force(double *dim, long long int pid,
if (r2 < p.h * p.h * kernel_gamma2 ||
r2 < parts[k].h * parts[k].h * kernel_gamma2) {
hydro_reset_acceleration(&p);
- runner_iact_nonsym_force(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
+ runner_iact_nonsym_force(r2, fdx, p.h, parts[k].h, &p, &parts[k], a, H);
}
}