Skip to content
GitLab
Commit 46ea8d13 authored by Tom Theuns's avatar Tom Theuns
Browse files

worked in progress

parent 48e3f35c
Hide whitespace changes
Inline Side-by-side
src/const.h
View file @ 46ea8d13
......@@ -20,6 +20,8 @@
#ifndef SWIFT_CONST_H
#define SWIFT_CONST_H
#include "physical_constants.h"
/* Hydrodynamical constants. */
#define const_hydro_gamma (5.0f / 3.0f)
......@@ -56,7 +58,6 @@
0.57735f /* Opening criteria, which is the ratio of the \
cell distance over the cell width. */
#define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */
#define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */
#define const_iepsilon2 (const_iepsilon* const_iepsilon)
......@@ -79,8 +80,18 @@
/* System of units */
#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
#define const_unit_mass_in_cgs 1 /* 1.9891e33 /\* 1 M_sun *\/ */
#define const_unit_velocity_in_cgs 1 /* 1e5 /\* km s^-1 *\/ */
#define const_unit_length_in_cgs (1000 * PARSEC_IN_CGS ) /* kpc */
#define const_unit_mass_in_cgs (SOLAR_MASS_IN_CGS) /* solar mass */
#define const_unit_velocity_in_cgs (1e5) /* km/s */
/* Derived constants */
#define const_unit_time_in_cgs (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
#define const_G (NEWTON_GRAVITY_CGS * const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/ (const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs))
/* External Potential Constants */
#define External_Potential_X (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs)
#define External_Potential_Y (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs)
#define External_Potential_Z (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs)
#define External_Potential_Mass (1e10 * SOLAR_MASS_IN_CGS / const_unit_mass_in_cgs)
#endif /* SWIFT_CONST_H */
src/engine.c
View file @ 46ea8d13
......@@ -119,6 +119,14 @@ void engine_mkghosts(struct engine *e, struct cell *c, struct cell *super) {
/* Add the kick task. */
c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0, c,
NULL, 0);
/* /\* Add the gravity tasks *\/ */
/* c->grav_external = scheduler_addtask(s, task_type_grav_external, task_subtype_none, 0, 0, */
/* c, NULL, 0); */
/* /\* Enforce gravity calculated before kick *\/ */
/* scheduler_addunlock(s, c->grav_external, c->kick); */
}
}
......@@ -775,6 +783,7 @@ void engine_maketasks(struct engine *e) {
}
}
#ifdef GRAVITY
#ifdef DEFAULT_GRAVITY
/* Add the gravity mm tasks. */
for (i = 0; i < nr_cells; i++)
......@@ -786,6 +795,20 @@ void engine_maketasks(struct engine *e) {
scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
&cells[i], &cells[j], 0);
}
#endif
/* #ifdef EXTERNAL_POTENTIAL */
/* /\* Add the external potential gravity *\/ */
/* for (i = 0; i < nr_cells; i++) */
/* if (cells[i].gcount > 0) { */
/* scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0, */
/* &cells[i], NULL, 0); */
/* for (j = i + 1; j < nr_cells; j++) */
/* if (cells[j].gcount > 0) */
/* scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0, */
/* &cells[i], &cells[j], 0); */
/* } */
/* #endif */
#endif
scheduler_print_tasks(sched, "sched.txt");
......
src/gravity/Default/runner_iact_grav.h
View file @ 46ea8d13
......@@ -21,9 +21,10 @@
#define SWIFT_RUNNER_IACT_GRAV_H
/* Includes. */
#include "const.h"
#include "kernel.h"
#include "vector.h"
#include "../../const.h"
#include "../../kernel.h"
#include "../../vector.h"
#include "../../timers.h"
/**
* @file runner_iact_grav.h
......@@ -31,6 +32,8 @@
*
*/
/**
* @brief Gravity potential
*/
......
src/gravity/ExternalPotential/runner_iact_grav.h
View file @ 46ea8d13
......@@ -26,7 +26,7 @@
#include "vector.h"
/**
* @file runner_iact_grav.h
* @file runner_iact_grav.c
* @brief Gravity interaction functions.
*
*/
......
src/runner.c
View file @ 46ea8d13
......@@ -68,6 +68,11 @@ 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 OUT if(CELL_ID == MY_CELL) {message(" cell %d %d %f \n",CELL_ID, c->count, 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
#include "runner_doiact.h"
......@@ -80,6 +85,69 @@ 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
*
* @param runner task
* @param cell
*/
void runner_dograv_external(struct runner *r, struct cell *c) {
struct part *p, *parts = c->parts;
float rinv;
float L[3], E;
int i, k, count = c->count;
float t_end = r->e->time;
//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++)
if (c->progeny[k] != NULL) runner_dograv_external(r, c->progeny[k]);
return;
}
#ifdef TASK_VERBOSE
OUT
#endif
/* /\* Loop over the parts in this cell. *\/ */
/* for (i = 0; i < count; i++) { */
/* /\* Get a direct pointer on the part. *\/ */
/* p = &parts[i]; */
/* /\* Is this part within the time step? *\/ */
/* if (p->t_end <= t_end) { */
/* rinv = 1 / sqrtf((p->x[0]-External_Potential_X)*(p->x[0]-External_Potential_X) + (p->x[1]-External_Potential_Y)*(p->x[1]-External_Potential_Y) + (p->x[2]-External_Potential_Z)*(p->x[2]-External_Potential_Z)); */
/* /\* check for energy and angular momentum conservation *\/ */
/* E = 0.5 * ((p->v[0]*p->v[0]) + (p->v[1]*p->v[1]) + (p->v[2]*p->v[2])) - const_G * External_Potential_Mass * rinv; */
/* L[0] = (p->x[1] - External_Potential_X)*p->v[2] - (p->x[2] - External_Potential_X)*p->v[1]; */
/* L[1] = (p->x[2] - External_Potential_Y)*p->v[0] - (p->x[0] - External_Potential_Y)*p->v[2]; */
/* L[2] = (p->x[0] - External_Potential_Z)*p->v[1] - (p->x[1] - External_Potential_Z)*p->v[0]; */
/* if(p->id == 0) { */
/* message("update %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\n", r->e->time, 1./rinv, p->x[0], p->x[1], p->x[2], E, L[0], L[1], L[2]); */
/* message(" ... %f %f %f\n", p->v[0], p->v[1], p->v[2]); */
/* message(" ... %e %e\n", const_G, External_Potential_Mass); */
/* } */
/* p->grav_accel[0] = - const_G * External_Potential_Mass * (p->x[0] - External_Potential_X) * rinv * rinv * rinv; */
/* p->grav_accel[1] = - const_G * External_Potential_Mass * (p->x[1] - External_Potential_Y) * rinv * rinv * rinv; */
/* p->grav_accel[2] = - const_G * External_Potential_Mass * (p->x[2] - External_Potential_Z) * rinv * rinv * rinv; */
/* } */
/* } */
/* TIMER_TOC(timer_dograv_external); */
}
/**
* @brief Sort the entries in ascending order using QuickSort.
*
......
src/task.c
View file @ 46ea8d13
......@@ -45,7 +45,7 @@
const char *taskID_names[task_type_count] = {
"none", "sort", "self", "pair", "sub", "init",
"ghost", "drift", "kick", "send", "recv", "grav_pp",
"grav_mm", "grav_up", "grav_down", "psort", "split_cell", "rewait"};
"grav_mm", "grav_up", "grav_down", "grav_external", "psort", "split_cell", "rewait"};
const char *subtaskID_names[task_type_count] = {"none", "density",
"force", "grav"};
......
src/task.h
View file @ 46ea8d13
......@@ -25,6 +25,7 @@
#include "cycle.h"
/* Some constants. */
//#define TASK_VERBOSE
#define task_maxwait 3
#define task_maxunlock 15
......@@ -45,6 +46,7 @@ enum task_types {
task_type_grav_mm,
task_type_grav_up,
task_type_grav_down,
task_type_grav_external,
task_type_psort,
task_type_split_cell,
task_type_rewait,
......
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