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

Gadget-2 force loop now correct

parent 6bfdc639
...@@ -1818,7 +1818,6 @@ if ( e->nodeID == 0 ) ...@@ -1818,7 +1818,6 @@ if ( e->nodeID == 0 )
printParticle(e->s->parts, 1000, e->s->nr_parts); printParticle(e->s->parts, 1000, e->s->nr_parts);
printParticle(e->s->parts, 515050, e->s->nr_parts); printParticle(e->s->parts, 515050, e->s->nr_parts);
exit(0);
/* Move forward in time */ /* Move forward in time */
e->timeOld = e->time; e->timeOld = e->time;
......
...@@ -174,7 +174,10 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struc ...@@ -174,7 +174,10 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struc
p->a[2] = 0.0f; p->a[2] = 0.0f;
/* Reset the time derivatives. */ /* Reset the time derivatives. */
p->v_sig = 0.0f; p->entropy_dt = 0.0f;
/* Reset minimal signal velocity */
p->v_sig = FLT_MAX;
} }
...@@ -212,6 +215,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(struct par ...@@ -212,6 +215,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(struct par
*/ */
__attribute__((always_inline)) INLINE static void hydro_end_force(struct part* p) { __attribute__((always_inline)) INLINE static void hydro_end_force(struct part* p) {
p->entropy_dt *= (const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
} }
/** /**
......
...@@ -200,80 +200,82 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( ...@@ -200,80 +200,82 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
__attribute__((always_inline)) INLINE static void runner_iact_force( __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 *pi, struct part *pj) {
/* float r = sqrtf(r2), ri = 1.0f / r; */ float wi, wj, wi_dx, wj_dx;
/* float xi, xj; */
/* float hi_inv, hi2_inv; */ const float fac_mu = 1.f; /* Will change with cosmological integration */
/* float hj_inv, hj2_inv; */
/* float wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, w, dvdr; */ const float r = sqrtf(r2);
/* float mi, mj, POrho2i, POrho2j, rhoi, rhoj; */ const float r_inv = 1.0f / r;
/* float v_sig, omega_ij, Pi_ij; */
/* float f; */ /* Get some values in local variables. */
/* int k; */ //const float mi = pi->mass;
const float mj = pj->mass;
/* /\* Get some values in local variables. *\/ */ const float rhoi = pi->rho;
/* mi = pi->mass; */ const float rhoj = pj->rho;
/* mj = pj->mass; */ const float pressurei = pi->pressure;
/* rhoi = pi->rho; */ const float pressurej = pj->pressure;
/* rhoj = pj->rho; */
/* //POrho2i = pi->force.POrho2; */ /* Get the kernel for hi. */
/* //POrho2j = pj->force.POrho2; */ const float hi_inv = 1.0f / hi;
const float hi2_inv = hi_inv * hi_inv;
/* /\* Get the kernel for hi. *\/ */ const float ui = r * hi_inv;
/* hi_inv = 1.0f / hi; */ kernel_deval(ui, &wi, &wi_dx);
/* hi2_inv = hi_inv * hi_inv; */ const float wi_dr = hi2_inv * hi2_inv * wi_dx;
/* xi = r * hi_inv; */
/* kernel_deval(xi, &wi, &wi_dx); */ /* Get the kernel for hj. */
/* wi_dr = hi2_inv * hi2_inv * wi_dx; */ const float hj_inv = 1.0f / hj;
const float hj2_inv = hj_inv * hj_inv;
/* /\* Get the kernel for hj. *\/ */ const float xj = r * hj_inv;
/* hj_inv = 1.0f / hj; */ kernel_deval(xj, &wj, &wj_dx);
/* hj2_inv = hj_inv * hj_inv; */ const float wj_dr = hj2_inv * hj2_inv * wj_dx;
/* xj = r * hj_inv; */
/* kernel_deval(xj, &wj, &wj_dx); */ /* Compute gradient terms */
/* wj_dr = hj2_inv * hj2_inv * wj_dx; */ const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* /\* Compute dv dot r. *\/ */
/* dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + */ /* Compute sound speeds */
/* (pi->v[2] - pj->v[2]) * dx[2]; */ const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
/* dvdr *= ri; */ const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
float v_sig = ci + cj;
/* /\* Compute the relative velocity. (This is 0 if the particles move away from */
/* * each other and negative otherwise) *\/ */ /* Compute dv dot r. */
/* omega_ij = fminf(dvdr, 0.f); */ 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];
/* /\* Compute signal velocity *\/ */
/* v_sig = pi->force.c + pj->force.c - 3.0f * omega_ij; */ /* Artificial viscosity term */
float visc = 0.f;
/* /\* Compute viscosity tensor *\/ */ if (dvdr < 0.f) {
/* Pi_ij = -const_viscosity_alpha * v_sig * omega_ij / (rhoi + rhoj); */ const float mu_ij = fac_mu * dvdr * r_inv;
v_sig -= 3.f * mu_ij;
/* /\* Apply balsara switch *\/ */ const float rho_ij = 0.5f * (rhoi + rhoj);
/* Pi_ij *= (pi->force.balsara + pj->force.balsara); */ const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v + 0.0001 * ci / fac_mu / hi);
const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v + 0.0001 * cj / fac_mu / hj);
/* /\* Get the common factor out. *\/ */ visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij * (balsara_i + balsara_j);
/* w = ri * */ }
/* ((POrho2i * wi_dr + POrho2j * wj_dr) + 0.25f * Pi_ij * (wi_dr + wj_dr)); */
/* Now, convolve with the kernel */
/* /\* Use the force, Luke! *\/ */ const float visc_term = 0.5f * mj * visc * (wi_dr + wj_dr) * r_inv;
/* for (k = 0; k < 3; k++) { */ const float sph_term = mj * (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
/* f = dx[k] * w; */
/* pi->a[k] -= mj * f; */ /* Eventually got the acceleration */
/* pj->a[k] += mi * f; */ const float acc = visc_term + sph_term;
/* } */
/* Use the force Luke ! */
/* /\* Get the time derivative for u. *\/ */ pi->a[0] -= acc * dx[0];
/* pi->force.u_dt += */ pi->a[1] -= acc * dx[1];
/* mj * dvdr * (POrho2i * wi_dr + 0.125f * Pi_ij * (wi_dr + wj_dr)); */ pi->a[2] -= acc * dx[2];
/* pj->force.u_dt += */
/* mi * dvdr * (POrho2j * wj_dr + 0.125f * Pi_ij * (wi_dr + wj_dr)); */ pj->a[0] += acc * dx[0];
pj->a[1] += acc * dx[1];
/* /\* Get the time derivative for h. *\/ */ pj->a[2] += acc * dx[2];
/* pi->force.h_dt -= mj * dvdr / rhoj * wi_dr; */
/* pj->force.h_dt -= mi * dvdr / rhoi * wj_dr; */ /* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig) ;
/* /\* Update the signal velocity. *\/ */ pj->v_sig = fmaxf(pj->v_sig, v_sig) ;
/* pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig); */
/* pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig); */ /* Change in entropy */
pi->entropy_dt += 0.5f * visc_term * dvdr;
pj->entropy_dt -= 0.5f * visc_term * dvdr;
} }
...@@ -284,76 +286,76 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -284,76 +286,76 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( __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 *pi, struct part *pj) {
/* float r = sqrtf(r2), ri = 1.0f / r; */ float wi, wj, wi_dx, wj_dx;
/* float xi, xj; */
/* float hi_inv, hi2_inv; */ const float fac_mu = 1.f; /* Will change with cosmological integration */
/* float hj_inv, hj2_inv; */
/* float wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, w, dvdr; */ const float r = sqrtf(r2);
/* float /\*mi,*\/ mj, POrho2i, POrho2j, rhoi, rhoj; */ const float r_inv = 1.0f / r;
/* float v_sig, omega_ij, Pi_ij; */
/* float f; */ /* Get some values in local variables. */
/* int k; */ //const float mi = pi->mass;
const float mj = pj->mass;
/* /\* Get some values in local variables. *\/ */ const float rhoi = pi->rho;
/* // mi = pi->mass; */ const float rhoj = pj->rho;
/* mj = pj->mass; */ const float pressurei = pi->pressure;
/* rhoi = pi->rho; */ const float pressurej = pj->pressure;
/* rhoj = pj->rho; */
/* POrho2i = pi->force.POrho2; */ /* Get the kernel for hi. */
/* POrho2j = pj->force.POrho2; */ const float hi_inv = 1.0f / hi;
const float hi2_inv = hi_inv * hi_inv;
/* /\* Get the kernel for hi. *\/ */ const float ui = r * hi_inv;
/* hi_inv = 1.0f / hi; */ kernel_deval(ui, &wi, &wi_dx);
/* hi2_inv = hi_inv * hi_inv; */ const float wi_dr = hi2_inv * hi2_inv * wi_dx;
/* xi = r * hi_inv; */
/* kernel_deval(xi, &wi, &wi_dx); */ /* Get the kernel for hj. */
/* wi_dr = hi2_inv * hi2_inv * wi_dx; */ const float hj_inv = 1.0f / hj;
const float hj2_inv = hj_inv * hj_inv;
/* /\* Get the kernel for hj. *\/ */ const float xj = r * hj_inv;
/* hj_inv = 1.0f / hj; */ kernel_deval(xj, &wj, &wj_dx);
/* hj2_inv = hj_inv * hj_inv; */ const float wj_dr = hj2_inv * hj2_inv * wj_dx;
/* xj = r * hj_inv; */
/* kernel_deval(xj, &wj, &wj_dx); */ /* Compute gradient terms */
/* wj_dr = hj2_inv * hj2_inv * wj_dx; */ const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* /\* Compute dv dot r. *\/ */
/* dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] + */ /* Compute sound speeds */
/* (pi->v[2] - pj->v[2]) * dx[2]; */ const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
/* dvdr *= ri; */ const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
float v_sig = ci + cj;
/* /\* Compute the relative velocity. (This is 0 if the particles move away from */
/* * each other and negative otherwise) *\/ */ /* Compute dv dot r. */
/* omega_ij = fminf(dvdr, 0.f); */ 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];
/* /\* Compute signal velocity *\/ */
/* v_sig = pi->force.c + pj->force.c - 3.0f * omega_ij; */ /* Artificial viscosity term */
float visc = 0.f;
/* /\* Compute viscosity tensor *\/ */ if (dvdr < 0.f) {
/* Pi_ij = -const_viscosity_alpha * v_sig * omega_ij / (rhoi + rhoj); */ const float mu_ij = fac_mu * dvdr * r_inv;
v_sig -= 3.f * mu_ij;
/* /\* Apply balsara switch *\/ */ const float rho_ij = 0.5f * (rhoi + rhoj);
/* Pi_ij *= (pi->force.balsara + pj->force.balsara); */ const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v + 0.0001 * ci / fac_mu / hi);
const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v + 0.0001 * cj / fac_mu / hj);
/* /\* Get the common factor out. *\/ */ visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij * (balsara_i + balsara_j);
/* w = ri * */ }
/* ((POrho2i * wi_dr + POrho2j * wj_dr) + 0.25f * Pi_ij * (wi_dr + wj_dr)); */
/* Now, convolve with the kernel */
/* /\* Use the force, Luke! *\/ */ const float visc_term = 0.5f * mj * visc * (wi_dr + wj_dr) * r_inv;
/* for (k = 0; k < 3; k++) { */ const float sph_term = mj * (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
/* f = dx[k] * w; */
/* pi->a[k] -= mj * f; */ /* Eventually got the acceleration */
/* } */ const float acc = visc_term + sph_term;
/* /\* Get the time derivative for u. *\/ */ /* Use the force Luke ! */
/* pi->force.u_dt += */ pi->a[0] -= acc * dx[0];
/* mj * dvdr * (POrho2i * wi_dr + 0.125f * Pi_ij * (wi_dr + wj_dr)); */ pi->a[1] -= acc * dx[1];
pi->a[2] -= acc * dx[2];
/* /\* Get the time derivative for h. *\/ */
/* pi->force.h_dt -= mj * dvdr / rhoj * wi_dr; */ /* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig) ;
/* /\* Update the signal velocity. *\/ */
/* pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig); */ /* Change in entropy */
/* pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig); */ pi->entropy_dt += 0.5f * visc_term * dvdr;
} }
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment