From bfb232eb1c68e274db50ff315d889e9050fb557a Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be>
Date: Thu, 11 Aug 2016 14:07:07 +0100
Subject: [PATCH] 1D and 2D version should now work. Do not produce expected
results however.
---
src/hydro/Gizmo/hydro.h | 13 ++++++----
src/hydro/Gizmo/hydro_gradients.h | 8 +------
src/hydro/Gizmo/hydro_gradients_sph.h | 34 ++++++++++++++-------------
src/hydro/Gizmo/hydro_iact.h | 19 +++++----------
src/hydro/Gizmo/hydro_io.h | 10 +++++++-
src/riemann.h | 2 +-
6 files changed, 43 insertions(+), 43 deletions(-)
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 01ddb6882f..1fbec63fca 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -116,7 +116,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Some smoothing length multiples. */
const float h = p->h;
const float ih = 1.0f / h;
- const float ih2 = ih * ih;
const float ihdim = pow_dimension(ih);
float volume;
@@ -144,7 +143,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
- hydro_gradients_prepare_force_loop(p, ih2, volume);
+ hydro_gradients_prepare_force_loop(p, ih, volume);
/* compute primitive variables */
/* eqns (3)-(5) */
@@ -390,14 +389,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {
+ return;
/* Set the hydro acceleration, based on the new momentum and mass */
/* NOTE: the momentum and mass are only correct for active particles, since
only active particles have received flux contributions from all their
neighbours. Since this method is only called for active particles,
this is indeed the case. */
- p->a_hydro[0] = p->conserved.momentum[0] / p->conserved.mass - p->v[0];
- p->a_hydro[1] = p->conserved.momentum[1] / p->conserved.mass - p->v[1];
- p->a_hydro[2] = p->conserved.momentum[2] / p->conserved.mass - p->v[2];
+ p->a_hydro[0] =
+ (p->conserved.momentum[0] / p->conserved.mass - p->v[0]) / p->force.dt;
+ p->a_hydro[1] =
+ (p->conserved.momentum[1] / p->conserved.mass - p->v[1]) / p->force.dt;
+ p->a_hydro[2] =
+ (p->conserved.momentum[2] / p->conserved.mass - p->v[2]) / p->force.dt;
}
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index f55bde5134..11f79ebff2 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -20,7 +20,7 @@
#ifndef SWIFT_HYDRO_GRADIENTS_H
#define SWIFT_HYDRO_GRADIENTS_H
-#define SPH_GRADIENTS
+//#define SPH_GRADIENTS
#include "hydro_slope_limiters.h"
@@ -168,12 +168,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
- // printf("WL: %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
- // printf("WR: %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
-
- // printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
- // printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
-
Wi[0] += dWi[0];
Wi[1] += dWi[1];
Wi[2] += dWi[2];
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index 29e36ec190..1671d9e43f 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/Gizmo/hydro_gradients_sph.h
@@ -157,28 +157,30 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_density_loop(
* @brief Calculations done before the force loop
*/
__attribute__((always_inline)) INLINE static void
-hydro_gradients_prepare_force_loop(struct part *p, float ih2, float volume) {
+hydro_gradients_prepare_force_loop(struct part *p, float ih, float volume) {
+
+ const float ihdimp1 = pow_dimension_plus_one(ih);
/* finalize gradients by multiplying with volume */
- p->primitives.gradients.rho[0] *= ih2 * ih2 * volume;
- p->primitives.gradients.rho[1] *= ih2 * ih2 * volume;
- p->primitives.gradients.rho[2] *= ih2 * ih2 * volume;
+ p->primitives.gradients.rho[0] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[1] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[2] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][0] *= ih2 * ih2 * volume;
- p->primitives.gradients.v[0][1] *= ih2 * ih2 * volume;
- p->primitives.gradients.v[0][2] *= ih2 * ih2 * volume;
+ p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][0] *= ih2 * ih2 * volume;
- p->primitives.gradients.v[1][1] *= ih2 * ih2 * volume;
- p->primitives.gradients.v[1][2] *= ih2 * ih2 * volume;
+ p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][0] *= ih2 * ih2 * volume;
- p->primitives.gradients.v[2][1] *= ih2 * ih2 * volume;
- p->primitives.gradients.v[2][2] *= ih2 * ih2 * volume;
+ p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
- p->primitives.gradients.P[0] *= ih2 * ih2 * volume;
- p->primitives.gradients.P[1] *= ih2 * ih2 * volume;
- p->primitives.gradients.P[2] *= ih2 * ih2 * volume;
+ p->primitives.gradients.P[0] *= ihdimp1 * volume;
+ p->primitives.gradients.P[1] *= ihdimp1 * volume;
+ p->primitives.gradients.P[2] *= ihdimp1 * volume;
hydro_slope_limit_cell(p);
}
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 82b26f3c78..d039ccf207 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -113,8 +113,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r = sqrtf(r2);
float xi, xj;
- float hi_inv, hi_inv2;
- float hj_inv, hj_inv2;
+ float hi_inv, hi_inv_dim;
+ float hj_inv, hj_inv_dim;
float wi, wj, wi_dx, wj_dx;
int k, l;
float A[3];
@@ -177,13 +177,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Compute kernel of pi. */
hi_inv = 1.0 / hi;
- hi_inv2 = hi_inv * hi_inv;
+ hi_inv_dim = pow_dimension(hi_inv);
xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
/* Compute kernel of pj. */
hj_inv = 1.0 / hj;
- hj_inv2 = hj_inv * hj_inv;
+ hj_inv_dim = pow_dimension(hj_inv);
xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
@@ -193,9 +193,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
for (k = 0; k < 3; k++) {
/* we add a minus sign since dx is pi->x - pj->x */
A[k] = -Vi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) * wi *
- hi_inv * hi_inv2 -
+ hi_inv_dim -
Vj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) * wj *
- hj_inv * hj_inv2;
+ hj_inv_dim;
Anorm += A[k] * A[k];
}
@@ -209,13 +209,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Anorm = sqrtf(Anorm);
for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
-#ifdef PRINT_ID
- if (pi->id == PRINT_ID || pj->id == PRINT_ID) {
- printf("pi: %g %g %g\npj: %g %g %g\nA = %g %g %g\n", pi->x[0], pi->x[1],
- pi->x[2], pj->x[0], pj->x[1], pj->x[2], A[0], A[1], A[2]);
- }
-#endif
-
/* Compute interface position (relative to pi, since we don't need the actual
* position) */
/* eqn. (8) */
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index bba7b5e6b9..dd7f46ceb0 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -56,6 +56,10 @@ float convert_u(struct engine* e, struct part* p) {
return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
}
+float convert_A(struct engine* e, struct part* p) {
+ return p->primitives.P / pow_gamma(p->primitives.rho);
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -66,7 +70,7 @@ float convert_u(struct engine* e, struct part* p) {
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 10;
+ *num_fields = 12;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -90,6 +94,10 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
geometry.volume);
list[9] = io_make_output_field("GradDensity", FLOAT, 3, UNIT_CONV_DENSITY,
parts, primitives.gradients.rho);
+ list[10] = io_make_output_field_convert_part(
+ "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, primitives.P, convert_A);
+ list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
+ parts, primitives.P);
}
/**
diff --git a/src/riemann.h b/src/riemann.h
index d647b02116..06acc2db88 100644
--- a/src/riemann.h
+++ b/src/riemann.h
@@ -27,7 +27,7 @@
#include "stdio.h"
#include "stdlib.h"
-#define HLLC_SOLVER
+#define EXACT_SOLVER
#ifdef EXACT_SOLVER
#include "riemann/riemann_exact.h"
--
GitLab