From 2e676ea67852d3ddda4bc7f55537a534e6b89a15 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be>
Date: Thu, 11 Aug 2016 17:03:11 +0100
Subject: [PATCH] Gizmo swift_fixdt works.
---
src/hydro/Gizmo/hydro.h | 26 +++++++++++++++++++-------
src/hydro/Gizmo/hydro_gradients_sph.h | 23 -----------------------
src/hydro/Gizmo/hydro_iact.h | 4 ++--
src/hydro/Gizmo/hydro_io.h | 4 +++-
src/hydro/Gizmo/hydro_part.h | 3 +++
src/hydro/Gizmo/hydro_slope_limiters.h | 4 ++--
6 files changed, 29 insertions(+), 35 deletions(-)
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 1fbec63fca..ca795ab008 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -171,6 +171,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Set the physical time step */
p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
+ /* Set the actual velocity of the particle */
+ p->force.v_full[0] = xp->v_full[0];
+ p->force.v_full[1] = xp->v_full[1];
+ p->force.v_full[2] = xp->v_full[2];
}
/**
@@ -389,18 +393,26 @@ __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->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;
+ if (p->force.dt) {
+ p->a_hydro[0] =
+ (p->conserved.momentum[0] / p->conserved.mass - p->primitives.v[0]) /
+ p->force.dt;
+ p->a_hydro[1] =
+ (p->conserved.momentum[1] / p->conserved.mass - p->primitives.v[1]) /
+ p->force.dt;
+ p->a_hydro[2] =
+ (p->conserved.momentum[2] / p->conserved.mass - p->primitives.v[2]) /
+ p->force.dt;
+ } else {
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+ }
}
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/Gizmo/hydro_gradients_sph.h
index 1671d9e43f..f8dc7ea004 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/Gizmo/hydro_gradients_sph.h
@@ -23,29 +23,6 @@
__attribute__((always_inline)) INLINE static void
hydro_gradients_init_density_loop(struct part *p) {
- /* use the old volumes to estimate new primitive variables to be used for the
- gradient calculation */
- if (p->conserved.mass) {
- p->primitives.rho = p->conserved.mass / p->geometry.volume;
- p->primitives.v[0] = p->conserved.momentum[0] / p->conserved.mass;
- p->primitives.v[1] = p->conserved.momentum[1] / p->conserved.mass;
- p->primitives.v[2] = p->conserved.momentum[2] / p->conserved.mass;
- p->primitives.P =
- hydro_gamma_minus_one *
- (p->conserved.energy -
- 0.5 * (p->conserved.momentum[0] * p->conserved.momentum[0] +
- p->conserved.momentum[1] * p->conserved.momentum[1] +
- p->conserved.momentum[2] * p->conserved.momentum[2]) /
- p->conserved.mass) /
- p->geometry.volume;
- } else {
- p->primitives.rho = 0.0f;
- p->primitives.v[0] = 0.0f;
- p->primitives.v[1] = 0.0f;
- p->primitives.v[2] = 0.0f;
- p->primitives.P = 0.0f;
- }
-
p->primitives.gradients.rho[0] = 0.0f;
p->primitives.gradients.rho[1] = 0.0f;
p->primitives.gradients.rho[2] = 0.0f;
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index d039ccf207..792d5c4ed0 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -135,8 +135,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Bi[k][l] = pi->geometry.matrix_E[k][l];
Bj[k][l] = pj->geometry.matrix_E[k][l];
}
- vi[k] = pi->v[k]; /* particle velocities */
- vj[k] = pj->v[k];
+ vi[k] = pi->force.v_full[k]; /* particle velocities */
+ vj[k] = pj->force.v_full[k];
}
Vi = pi->geometry.volume;
Vj = pj->geometry.volume;
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index dd7f46ceb0..c94451aad9 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -70,7 +70,7 @@ float convert_A(struct engine* e, struct part* p) {
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 12;
+ *num_fields = 13;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -98,6 +98,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
"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);
+ list[12] = io_make_output_field("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
+ parts, conserved.energy);
}
/**
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index 8980bd4a49..ce1d7a5e00 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -159,6 +159,9 @@ struct part {
/* Physical time step of the particle. */
float dt;
+ /* Actual velocity of the particle. */
+ float v_full[3];
+
} force;
/* Particle mass (this field is also part of the conserved quantities...). */
diff --git a/src/hydro/Gizmo/hydro_slope_limiters.h b/src/hydro/Gizmo/hydro_slope_limiters.h
index b341d2f15e..357ba2a84e 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters.h
@@ -20,8 +20,8 @@
#ifndef SWIFT_HYDRO_SLOPE_LIMITERS_H
#define SWIFT_HYDRO_SLOPE_LIMITERS_H
-//#define PER_FACE_LIMITER
-//#define CELL_WIDE_LIMITER
+#define PER_FACE_LIMITER
+#define CELL_WIDE_LIMITER
#ifdef PER_FACE_LIMITER
#include "hydro_slope_limiters_face.h"
--
GitLab