Skip to content
Snippets Groups Projects
Commit 9813ef6e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Don't be too verbose...

parent 1a991b27
No related branches found
No related tags found
1 merge request!143Gravity particles
Hide whitespace changes
Inline Side-by-side
src/gravity/Default/gravity.h
+ 14
− 3
View file @ 9813ef6e
......@@ -20,19 +20,25 @@
#include <float.h>
#include "potentials.h"
/**
* @brief Computes the gravity time-step of a given particle.
*
* This function only branches towards the potential chosen by the user.
*
* @param potential The properties of the external potential.
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float gravity_compute_timestep(
const struct external_potential* potential, const struct phys_const* const phys_const, const struct gpart* const gp) {
const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* const gp) {
float dt = FLT_MAX;
#ifdef EXTERNAL_POTENTIAL_POINTMASS
dt = fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
dt =
fminf(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
#endif
return dt;
......@@ -79,10 +85,15 @@ __attribute__((always_inline))
/**
* @brief Computes the gravitational acceleration induced by external potentials
*
* This function only branches towards the potential chosen by the user.
*
* @param potential The properties of the external potential.
* @param phys_const The physical constants in internal units.
* @param gp The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void external_gravity(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* gp) {
__attribute__((always_inline)) INLINE static void external_gravity(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* gp) {
#ifdef EXTERNAL_POTENTIAL_POINTMASS
external_gravity_pointmass(potential, phys_const, gp);
......
src/gravity/Default/gravity_debug.h
+ 2
− 1
View file @ 9813ef6e
......@@ -24,5 +24,6 @@ __attribute__((always_inline))
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"mass=%.3e 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->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin, p->ti_end);
p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin,
p->ti_end);
}
src/gravity/Default/gravity_part.h
+ 3
− 4
View file @ 9813ef6e
......@@ -40,10 +40,9 @@ struct gpart {
/* Particle time of end of time-step. */
int ti_end;
/* current time of x, and of v_full */
float tx;
float tv;
/* /\* current time of x, and of v_full *\/ */
/* float tx; */
/* float tv; */
/* Anonymous union for id/part. */
union {
......
src/runner.c
+ 57
− 63
View file @ 9813ef6e
......@@ -71,19 +71,20 @@ const float runner_shift[13 * 3] = {
/* Does the axis need flipping ? */
const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
#define MY_CELL 428428428
#define DX 0.1
#define NX 1000
#define CELL_ID \
((int)(c->loc[0] / DX) * NX *NX + (int)(c->loc[1] / DX) * NX + \
(int)(c->loc[2] / DX))
#define OUT \
if (CELL_ID == MY_CELL) { \
message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount, r->e->time); \
}
//#define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time);
//#define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d
//%f\n",c->loc[0],c->loc[1],c->loc[2], CELL_ID, c->count, r->e->time);
/* #define MY_CELL 428428428 */
/* #define DX 0.1 */
/* #define NX 1000 */
/* #define CELL_ID \ */
/* ((int)(c->loc[0] / DX) * NX *NX + (int)(c->loc[1] / DX) * NX + \ */
/* (int)(c->loc[2] / DX)) */
/* #define OUT \ */
/* if (CELL_ID == MY_CELL) { \ */
/* message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount, r->e->time); \ */
/* } */
/* #define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time); */
/* #define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d */
/* %f\n",c->loc[0],c->loc[1],c->loc[2], CELL_ID, c->count, r->e->time); */
/* Import the density loop functions. */
#define FUNCTION density
......@@ -106,20 +107,11 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
void runner_dograv_external(struct runner *r, struct cell *c) {
struct gpart *g, *gparts = c->gparts;
float L[3], E;
int i, k, gcount = c->gcount;
const int ti_current = r->e->ti_current;
// struct space *s = r->e->s;
// double CentreOfPotential[3];
TIMER_TIC;
/* /\* location of external gravity point mass - should pass in as
* paraneter *\/ */
/* CentreOfPotential[0] = 0.5 * s->dim[0]; */
/* CentreOfPotential[1] = 0.5 * s->dim[1]; */
/* CentreOfPotential[2] = 0.5 * s->dim[2]; */
/* Recurse? */
if (c->split) {
for (k = 0; k < 8; k++)
......@@ -130,6 +122,7 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
#ifdef TASK_VERBOSE
OUT;
#endif
/* Loop over the parts in this cell. */
for (i = 0; i < gcount; i++) {
......@@ -138,49 +131,50 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
/* Is this part within the time step? */
if (g->ti_end <= ti_current) {
// external_gravity_pointmass(e->physical_constants,
// potential_constants, g);
external_gravity(r->e->potential, r->e->physical_constants, g);
/* check for energy and angular momentum conservation - begin by
* synchronizing velocity*/
const float dx = g->x[0] - r->e->potential->point_mass.x;
const float dy = g->x[1] - r->e->potential->point_mass.y;
const float dz = g->x[2] - r->e->potential->point_mass.z;
const float dr = sqrtf((dx * dx) + (dy * dy) + (dz * dz));
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const int current_dti = g->ti_end - g->ti_begin;
const float dt = 0.5f * current_dti * r->e->timeBase;
const float vx = g->v_full[0] + dt * g->a_grav[0];
const float vy = g->v_full[1] + dt * g->a_grav[1];
const float vz = g->v_full[2] + dt * g->a_grav[2];
/* E/L */
E = 0.5 * ((vx * vx) + (vy * vy) + (vz * vz)) -
r->e->physical_constants->newton_gravity *
r->e->potential->point_mass.mass * rinv;
L[0] = dy * vz - dz * vy;
L[1] = dz * vx - dx * vz;
L[2] = dx * vy - dy * vx;
if (abs(g->id) == 1) {
float v2 = vx * vx + vy * vy + vz * vz;
float fg = r->e->physical_constants->newton_gravity *
r->e->potential->point_mass.mass * rinv;
float fga = sqrtf((g->a_grav[0] * g->a_grav[0]) +
(g->a_grav[1] * g->a_grav[1]) +
(g->a_grav[2] * g->a_grav[2])) *
dr;
// message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]= %f
// x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2], g->x[0],
// g->x[1], vx, vy);
message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time,
g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1],
vx, vy);
/* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity,
* r->e->potential->point_mass.mass); */
/* exit(-1); */
}
/* /\* check for energy and angular momentum conservation - begin by */
/* * synchronizing velocity*\/ */
/* const float dx = g->x[0] - r->e->potential->point_mass.x; */
/* const float dy = g->x[1] - r->e->potential->point_mass.y; */
/* const float dz = g->x[2] - r->e->potential->point_mass.z; */
/* const float dr = sqrtf((dx * dx) + (dy * dy) + (dz * dz)); */
/* const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz); */
/* const int current_dti = g->ti_end - g->ti_begin; */
/* const float dt = 0.5f * current_dti * r->e->timeBase; */
/* const float vx = g->v_full[0] + dt * g->a_grav[0]; */
/* const float vy = g->v_full[1] + dt * g->a_grav[1]; */
/* const float vz = g->v_full[2] + dt * g->a_grav[2]; */
/* /\* E/L *\/ */
/* float L[3], E; */
/* E = 0.5 * ((vx * vx) + (vy * vy) + (vz * vz)) - */
/* r->e->physical_constants->newton_gravity * */
/* r->e->potential->point_mass.mass * rinv; */
/* L[0] = dy * vz - dz * vy; */
/* L[1] = dz * vx - dx * vz; */
/* L[2] = dx * vy - dy * vx; */
/* if (abs(g->id) == 1) { */
/* float v2 = vx * vx + vy * vy + vz * vz; */
/* float fg = r->e->physical_constants->newton_gravity * */
/* r->e->potential->point_mass.mass * rinv; */
/* float fga = sqrtf((g->a_grav[0] * g->a_grav[0]) + */
/* (g->a_grav[1] * g->a_grav[1]) + */
/* (g->a_grav[2] * g->a_grav[2])) * */
/* dr; */
/* // message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]= %f */
/* // x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2], g->x[0], */
/* // g->x[1], vx, vy); */
/* message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time, */
/* g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1], */
/* vx, vy); */
/* /\* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity, */
/* * r->e->potential->point_mass.mass); *\/ */
/* /\* exit(-1); *\/ */
/* } */
}
}
TIMER_TOC(timer_dograv_external);
......
src/task.h
+ 0
− 1
View file @ 9813ef6e
......@@ -28,7 +28,6 @@
#include "cycle.h"
/* Some constants. */
#define TASK_VERBOSE
#define task_maxwait 3
#define task_maxunlock 15
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment