Skip to content
GitLab
Commit 66959f80 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Code formatting

parent 8fb542bd
Hide whitespace changes
Inline Side-by-side
src/runner.c
View file @ 66959f80
......@@ -68,13 +68,18 @@ const float runner_shift[13 * 3] = {
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 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 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
......@@ -88,7 +93,6 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
/* Import the gravity loop functions. */
#include "runner_doiact_grav.h"
/**
* @brief Calculate gravity acceleration from external potential
*
......@@ -101,20 +105,22 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
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];
// struct space *s = r->e->s;
// double CentreOfPotential[3];
TIMER_TIC;
/* /\* location of external gravity point mass - should pass in as paraneter *\/ */
/* /\* 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]; */
message(" (x,y,z) = (%e, %e, %e), M= %e", r->e->potential->point_mass.x, r->e->potential->point_mass.y, r->e->potential->point_mass.z, r->e->potential->point_mass.mass);
message(" (x,y,z) = (%e, %e, %e), M= %e", r->e->potential->point_mass.x,
r->e->potential->point_mass.y, r->e->potential->point_mass.z,
r->e->potential->point_mass.mass);
exit(-1);
/* Recurse? */
if (c->split) {
for (k = 0; k < 8; k++)
......@@ -123,7 +129,7 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
}
#ifdef TASK_VERBOSE
OUT
OUT;
#endif
/* Loop over the parts in this cell. */
for (i = 0; i < gcount; i++) {
......@@ -133,48 +139,57 @@ 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_pointmass(r->e->physical_constants, g);
/* check for energy and angular momentum conservation - begin by synchronizing velocity*/
const float dx = g->x[0]-External_Potential_X;
const float dy = g->x[1]-External_Potential_Y;
const float dz = g->x[2]-External_Potential_Z;
const float dr = sqrtf((dx*dx) + (dy*dy) + (dz*dz));
const float rinv = 1.f / sqrtf(dx*dx + dy*dy + dz*dz);
/* /\* current acceleration *\/ */
/* g->a_grav[0] = - const_G * External_Potential_Mass * dx * rinv * rinv * rinv; */
/* g->a_grav[1] = - const_G * External_Potential_Mass * dy * rinv * rinv * rinv; */
/* g->a_grav[2] = - const_G * External_Potential_Mass * dz * rinv * rinv * rinv; */
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)) - const_G * External_Potential_Mass * rinv;
// external_gravity_pointmass(e->physical_constants,
// potential_constants, g);
external_gravity_pointmass(r->e->physical_constants, g);
/* check for energy and angular momentum conservation - begin by
* synchronizing velocity*/
const float dx = g->x[0] - External_Potential_X;
const float dy = g->x[1] - External_Potential_Y;
const float dz = g->x[2] - External_Potential_Z;
const float dr = sqrtf((dx * dx) + (dy * dy) + (dz * dz));
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
/* /\* current acceleration *\/ */
/* g->a_grav[0] = - const_G * External_Potential_Mass * dx * rinv * rinv
* * rinv; */
/* g->a_grav[1] = - const_G * External_Potential_Mass * dy * rinv * rinv
* * rinv; */
/* g->a_grav[2] = - const_G * External_Potential_Mass * dz * rinv * rinv
* * rinv; */
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)) -
const_G * External_Potential_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 = const_G * External_Potential_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\n", r->e->time, g->tx, g->tv, dt, v2, fg, fga, E, L[2], g->x[0], g->x[1], vx, vy);
if (abs(g->id) == 1) {
float v2 = vx * vx + vy * vy + vz * vz;
float fg = const_G * External_Potential_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\n", r->e->time, g->tx,
g->tv, dt, v2, fg, fga, E, L[2], g->x[0], g->x[1], vx, vy);
// message(" G=%e M=%e\n", const_G, External_Potential_Mass);
}
}
}
TIMER_TOC(timer_dograv_external);
}
/**
* @brief Sort the entries in ascending order using QuickSort.
*
......@@ -780,7 +795,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
TIMER_TIC
#ifdef TASK_VERBOSE
OUT
OUT;
#endif
/* No children? */
......@@ -902,7 +917,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
TIMER_TIC
#ifdef TASK_VERBOSE
OUT
OUT;
#endif
/* No children? */
......@@ -960,7 +975,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
const double dt = (ti_end - ti_start) * timeBase;
const double half_dt = (ti_end - gp->ti_end) * timeBase;
/* Move particle forward in time */
/* Move particle forward in time */
gp->ti_begin = gp->ti_end;
gp->ti_end = gp->ti_begin + new_dti;
......@@ -976,10 +991,9 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
g_updated++;
}
/* Minimal time for next end of time-step */
ti_end_min = min(gp->ti_end, ti_end_min);
ti_end_max = max(gp->ti_end, ti_end_max);
/* Minimal time for next end of time-step */
ti_end_min = min(gp->ti_end, ti_end_min);
ti_end_max = max(gp->ti_end, ti_end_max);
}
/* Now do the hydro ones... */
......@@ -1017,7 +1031,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Compute the next timestep (gravity condition) */
float new_dt_grav = FLT_MAX;
if (p->gpart != NULL)
new_dt_grav = gravity_compute_timestep(r->e->physical_constants, p->gpart);
new_dt_grav =
gravity_compute_timestep(r->e->physical_constants, p->gpart);
float new_dt = fminf(new_dt_hydro, new_dt_grav);
......@@ -1121,7 +1136,6 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
ti_end_min = min(p->ti_end, ti_end_min);
ti_end_max = max(p->ti_end, ti_end_max);
}
}
/* Otherwise, aggregate data from children. */
......@@ -1283,9 +1297,9 @@ void *runner_main(void *data) {
case task_type_grav_down:
runner_dograv_down(r, t->ci);
break;
case task_type_grav_external:
runner_dograv_external(r, t->ci);
break;
case task_type_grav_external:
runner_dograv_external(r, t->ci);
break;
case task_type_psort:
space_do_parts_sort();
break;
......
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