Commit 13db1cd5 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Star particles are now correctly initialised at the start of a run.

parent b9304880
......@@ -77,6 +77,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \
riemann/riemann_hllc.h riemann/riemann_trrs.h \
riemann/riemann_exact.h riemann/riemann_vacuum.h \
stars/Default/star.h stars/Default/star_iact.h stars/Default/star_io.h \
stars/Default/star_debug.h stars/Default/star_part.h \
potential/none/potential.h potential/point_mass/potential.h \
potential/isothermal/potential.h potential/disc_patch/potential.h \
potential/softened_isothermal/potential.h \
......
......@@ -52,6 +52,7 @@
#include "memswap.h"
#include "minmax.h"
#include "runner.h"
#include "stars.h"
#include "threadpool.h"
#include "tools.h"
......@@ -1748,6 +1749,32 @@ void space_init_gparts(struct space *s) {
}
}
/**
* @brief Initialises all the s-particles by setting them into a valid state
*
* Calls star_first_init_spart() on all the particles
*/
void space_init_sparts(struct space *s) {
const size_t nr_sparts = s->nr_sparts;
struct spart *restrict sp = s->sparts;
for (size_t i = 0; i < nr_sparts; ++i) {
#ifdef HYDRO_DIMENSION_2D
sp[i].x[2] = 0.f;
sp[i].v[2] = 0.f;
#endif
#ifdef HYDRO_DIMENSION_1D
sp[i].x[1] = sp[i].x[2] = 0.f;
sp[i].v[1] = sp[i].v[2] = 0.f;
#endif
star_first_init_spart(&sp[i]);
}
}
/**
* @brief Split the space into cells given the array of particles.
*
......@@ -1916,7 +1943,7 @@ void space_init(struct space *s, const struct swift_params *params,
space_init_parts(s);
space_init_xparts(s);
space_init_gparts(s);
// space_init_sparts(s); // MATTHIEU
space_init_sparts(s);
/* Init the space lock. */
if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
......
......@@ -189,6 +189,7 @@ void space_do_parts_sort();
void space_do_gparts_sort();
void space_init_parts(struct space *s);
void space_init_gparts(struct space *s);
void space_init_sparts(struct space *s);
void space_link_cleanup(struct space *s);
void space_check_drift_point(struct space *s, int ti_current);
void space_clean(struct space *s);
......
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_STAR_H
#define SWIFT_STAR_H
/* Config parameters. */
#include "../config.h"
/* So far only one model here */
/* Straight-forward import */
#include "./stars/Default/star.h"
#include "./stars/Default/star_iact.h"
#endif
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 2016 Tom Theuns (tom.theuns@durham.ac.uk)
* Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
......@@ -17,89 +16,64 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_GRAVITY_H
#define SWIFT_DEFAULT_GRAVITY_H
#ifndef SWIFT_DEFAULT_STAR_H
#define SWIFT_DEFAULT_STAR_H
#include <float.h>
#include "minmax.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
* @brief Computes the gravity time-step of a given star particle.
*
* @param gp Pointer to the g-particle data.
* @param sp Pointer to the s-particle data.
*/
__attribute__((always_inline)) INLINE static float
gravity_compute_timestep_self(const struct gpart* const gp) {
__attribute__((always_inline)) INLINE static float star_compute_timestep_self(
const struct spart* const sp) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
gp->a_grav[2] * gp->a_grav[2];
const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
const float dt = sqrtf(2.f * const_gravity_eta * gp->epsilon / ac);
return dt;
return FLT_MAX;
}
/**
* @brief Initialises the g-particles for the first time
* @brief Initialises the s-particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param gp The particle to act upon
* @param sp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp) {
__attribute__((always_inline)) INLINE static void star_first_init_spart(
struct spart* sp) {
gp->ti_begin = 0;
gp->ti_end = 0;
gp->epsilon = 0.; // MATTHIEU
sp->ti_begin = 0;
sp->ti_end = 0;
}
/**
* @brief Prepares a g-particle for the gravity calculation
*
* Zeroes all the relevant arrays in preparation for the sums taking place in
* the variaous tasks
* @brief Prepares a s-particle for its interactions
*
* @param gp The particle to act upon
* @param sp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void gravity_init_gpart(
struct gpart* gp) {
/* Zero the acceleration */
gp->a_grav[0] = 0.f;
gp->a_grav[1] = 0.f;
gp->a_grav[2] = 0.f;
}
__attribute__((always_inline)) INLINE static void star_init_spart(
struct spart* sp) {}
/**
* @brief Finishes the gravity calculation.
* @brief Finishes the calculation of (non-gravity) forces acting on stars
*
* Multiplies the forces and accelerations by the appropiate constants
*
* @param gp The particle to act upon
* @param const_G Newton's constant in internal units
* @param sp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void gravity_end_force(
struct gpart* gp, float const_G) {
/* Let's get physical... */
gp->a_grav[0] *= const_G;
gp->a_grav[1] *= const_G;
gp->a_grav[2] *= const_G;
}
__attribute__((always_inline)) INLINE static void star_end_force(
struct spart* sp) {}
/**
* @brief Kick the additional variables
*
* @param gp The particle to act upon
* @param sp The particle to act upon
* @param dt The time-step for this kick
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
struct gpart* gp, float dt, float half_dt) {}
__attribute__((always_inline)) INLINE static void star_kick_extra(
struct spart* sp, float dt, float half_dt) {}
#endif /* SWIFT_DEFAULT_GRAVITY_H */
#endif /* SWIFT_DEFAULT_STAR_H */
......@@ -23,9 +23,9 @@ __attribute__((always_inline)) INLINE static void star_debug_particle(
const struct spart* p) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v_full=[%.3e,%.3e,%.3e] \n t_begin=%d, t_end=%d\n",
"v_full=[%.3e,%.3e,%.3e] p->mass=%.3e \n t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
p->ti_begin, p->ti_end);
p->mass, p->ti_begin, p->ti_end);
}
#endif /* SWIFT_DEFAULT_STAR_DEBUG_H */
Markdown is supported
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