Commit c46b19fc authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added cosmological terms to the Gadget2 hydro scheme.

parent 51642837
......@@ -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.;
......
......@@ -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;
......
......@@ -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. */
}
......
......@@ -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) {