Commit fb624a9a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Time integration is now correct. Code still hangs somewhere.

parent 0f955986
......@@ -29,6 +29,10 @@ tests/testGreetings
tests/testReading
tests/input.hdf5
tests/testSingle
tests/testTimeIntegration
tests/testSPHStep
theory/latex/swift.pdf
m4/libtool.m4
m4/ltoptions.m4
......
......@@ -362,7 +362,7 @@ void engine_repartition(struct engine *e) {
if (t->type != task_type_self && t->type != task_type_pair &&
t->type != task_type_sub && t->type != task_type_ghost &&
t->type != task_type_drift && t->type != task_type_kick &&
t->type != task_type_init)
t->type != task_type_init)
continue;
/* Get the task weight. */
......@@ -1370,16 +1370,45 @@ int engine_marktasks(struct engine *e) {
}
/**
* @brief Rebuild the space and tasks.
* @brief Prints the number of tasks in the engine
*
* @param e The #engine.
*/
void engine_rebuild(struct engine *e) {
void engine_print(struct engine *e) {
int k;
struct scheduler *sched = &e->sched;
/* Count and print the number of each task type. */
int counts[task_type_count + 1];
for (k = 0; k <= task_type_count; k++) counts[k] = 0;
for (k = 0; k < sched->nr_tasks; k++)
if (!sched->tasks[k].skip)
counts[(int)sched->tasks[k].type] += 1;
else
counts[task_type_count] += 1;
#ifdef WITH_MPI
printf("[%03i] engine_print: task counts are [ %s=%i", e->nodeID,
taskID_names[0], counts[0]);
#else
printf("engine_print: task counts are [ %s=%i", taskID_names[0], counts[0]);
#endif
for (k = 1; k < task_type_count; k++)
printf(" %s=%i", taskID_names[k], counts[k]);
printf(" skipped=%i ]\n", counts[task_type_count]);
fflush(stdout);
message("nr_parts = %i.", e->s->nr_parts);
}
/**
* @brief Rebuild the space and tasks.
*
* @param e The #engine.
*/
void engine_rebuild(struct engine *e) {
/* Clear the forcerebuild flag, whatever it was. */
e->forcerebuild = 0;
......@@ -1410,25 +1439,8 @@ void engine_rebuild(struct engine *e) {
// message( "engine_marktasks took %.3f ms." , (double)(getticks() -
// tic)/CPU_TPS*1000 );
/* Count and print the number of each task type. */
int counts[task_type_count + 1];
for (k = 0; k <= task_type_count; k++) counts[k] = 0;
for (k = 0; k < sched->nr_tasks; k++)
if (!sched->tasks[k].skip)
counts[(int)sched->tasks[k].type] += 1;
else
counts[task_type_count] += 1;
#ifdef WITH_MPI
printf("[%03i] engine_rebuild: task counts are [ %s=%i", e->nodeID,
taskID_names[0], counts[0]);
#else
printf("engine_rebuild: task counts are [ %s=%i", taskID_names[0], counts[0]);
#endif
for (k = 1; k < task_type_count; k++)
printf(" %s=%i", taskID_names[k], counts[k]);
printf(" skipped=%i ]\n", counts[task_type_count]);
fflush(stdout);
message("nr_parts = %i.", e->s->nr_parts);
/* Print the status of the system */
engine_print(e);
}
/**
......@@ -1746,9 +1758,9 @@ void engine_step(struct engine *e) {
struct cell *c;
struct space *s = e->s;
TIMER_TIC2
TIMER_TIC2;
message("Begin step: %d", e->step);
/* Re-distribute the particles amongst the nodes? */
if (e->forcerepart) engine_repartition(e);
......@@ -1756,33 +1768,20 @@ void engine_step(struct engine *e) {
/* Prepare the space. */
engine_prepare(e);
// engine_single_density( e->s->dim , 3392063069037 , e->s->parts ,
// e->s->nr_parts , e->s->periodic );
engine_print(e);
/* Send off the runners. */
TIMER_TIC
engine_launch(
e, e->nr_threads,
(1 << task_type_sort) | (1 << task_type_self) | (1 << task_type_pair) |
(1 << task_type_sub) | (1 << task_type_ghost) |
(1 << task_type_kick) | (1 << task_type_send) |
(1 << task_type_recv) |
/* (1 << task_type_grav_pp) | (1 << task_type_grav_mm) | */
/* (1 << task_type_grav_up) | (1 << task_type_grav_down) | */
(1 << task_type_link));
TIMER_TIC;
engine_launch(e, e->nr_threads,
(1 << task_type_sort) | (1 << task_type_self) |
(1 << task_type_pair) | (1 << task_type_sub) |
(1 << task_type_init) | (1 << task_type_ghost) |
(1 << task_type_kick) | (1 << task_type_send) |
(1 << task_type_recv) | (1 << task_type_link));
TIMER_TOC(timer_runners);
// engine_single_force( e->s->dim , 8328423931905 , e->s->parts ,
// e->s->nr_parts , e->s->periodic );
// for(k=0; k<10; ++k)
// printParticle(parts, k);
// printParticle( parts , 432626 );
// printParticle( e->s->parts , 3392063069037 , e->s->nr_parts );
// printParticle( e->s->parts , 8328423931905 , e->s->nr_parts );
/* Collect the cell data from the second kick. */
/* Collect the cell data from the kick. */
for (k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) {
c = &s->cells[k];
......@@ -1806,12 +1805,12 @@ void engine_step(struct engine *e) {
out[0] = t_end_min;
if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate dt_min.");
error("Failed to aggregate t_end_min.");
t_end_min = in[0];
out[0] = t_end_max;
if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate dt_max.");
error("Failed to aggregate t_end_max.");
t_end_max = in[0];
out[0] = count;
out[1] = ekin;
......@@ -1830,71 +1829,16 @@ if ( e->nodeID == 0 )
message( "nr_parts=%i." , nr_parts ); */
#endif
e->count_step = count;
e->ekin = ekin;
e->epot = epot;
// printParticle( e->s->parts , 382557 , e->s->nr_parts );
// message( "dt_min/dt_max is %e/%e." , dt_min , dt_max ); fflush(stdout);
// message( "etot is %e (ekin=%e, epot=%e)." , ekin+epot , ekin , epot );
// fflush(stdout);
// message( "total momentum is [ %e , %e , %e ]." , mom[0] , mom[1] , mom[2]
// ); fflush(stdout);
// message( "total angular momentum is [ %e , %e , %e ]." , ang[0] , ang[1] ,
// ang[2] ); fflush(stdout);
// message( "updated %i parts (dt_step=%.3e)." , count , dt_step );
// fflush(stdout);
/* Increase the step. */
/* Move forward in time */
e->timeOld = e->time;
e->time = t_end_min;
e->step += 1;
/* /\* Does the time step need adjusting? *\/ */
/* if (e->policy & engine_policy_fixdt) { */
/* dt = e->dt_orig; */
/* } else { */
/* if (dt == 0) { */
/* e->nullstep += 1; */
/* if (e->dt_orig > 0.0) { */
/* dt = e->dt_orig; */
/* while (dt_min < dt) dt *= 0.5; */
/* while (dt_min > 2 * dt) dt *= 2.0; */
/* } else */
/* dt = dt_min; */
/* for (k = 0; k < s->nr_parts; k++) { */
/* /\* struct part *p = &s->parts[k]; */
/* struct xpart *xp = &s->xparts[k]; */
/* float dt_curr = dt; */
/* for ( int j = (int)( p->dt / dt ) ; j > 1 ; j >>= 1 ) */
/* dt_curr *= 2.0f; */
/* xp->dt_curr = dt_curr; *\/ */
/* s->parts[k].dt = dt; */
/* s->xparts[k].dt_curr = dt; */
/* } */
/* // message( "dt_min=%.3e, adjusting time step to dt=%e." , dt_min ,
* e->dt */
/* // ); */
/* } else { */
/* while (dt_min < dt) { */
/* dt *= 0.5; */
/* e->step *= 2; */
/* e->nullstep *= 2; */
/* // message( "dt_min dropped below time step, adjusting to dt=%e." ,
*/
/* // e->dt ); */
/* } */
/* while (dt_min > 2 * dt && (e->step & 1) == 0) { */
/* dt *= 2.0; */
/* e->step /= 2; */
/* e->nullstep /= 2; */
/* // message( "dt_min is larger than twice the time step, adjusting to
*/
/* // dt=%e." , e->dt ); */
/* } */
/* } */
/* } */
/* e->dt = dt; */
message("e->time=%f", e->time);
message("Ready to drift");
/* /\* Set the system time. *\/ */
/* e->time = dt * (e->step - e->nullstep); */
/* Drift all particles */
engine_launch(e, e->nr_threads, (1 << task_type_drift));
TIMER_TOC2(timer_step);
}
......@@ -2067,7 +2011,6 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
float timeBegin, float timeEnd) {
int k;
float dt_min = dt;
#if defined(HAVE_SETAFFINITY)
int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
int i, j, cpuid[nr_cores];
......
......@@ -137,6 +137,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
int nr_queues, int nr_nodes, int nodeID, int policy,
float timeBegin, float timeEnd);
void engine_prepare(struct engine *e);
void engine_print(struct engine *e);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
void engine_split(struct engine *e, int *grid);
......
......@@ -18,7 +18,6 @@
*
******************************************************************************/
/* Includes. */
#ifndef SWIFT_MAP_H
#define SWIFT_MAP_H
......
......@@ -441,7 +441,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
* @param N_total Total number of particles across all cores
* @param mpi_rank The MPI task rank calling the function
* @param offset Offset in the array where this mpi task starts writing
* @param part A (char*) pointer on the first occurrence of the field of interest
* @param part A (char*) pointer on the first occurrence of the field of
*interest
*in the parts array
* @param field The name (code name) of the field to read from.
* @param us The UnitSystem currently in use
......@@ -582,7 +583,8 @@ void write_output_parallel(struct engine* e, struct UnitSystem* us,
N_total, mpi_rank, offset, u, us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
/* writeArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N, 1, parts, N_total, */
/* writeArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N, 1, parts,
* N_total, */
/* mpi_rank, offset, dt, us, UNIT_CONV_TIME); */
writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts,
N_total, mpi_rank, offset, a, us, UNIT_CONV_ACCELERATION);
......
......@@ -40,9 +40,6 @@ struct xpart {
/* Old density. */
float omega;
/* particle's current time-step. */
float dt_curr;
} __attribute__((aligned(32)));
/* Gravity particle. */
......@@ -115,46 +112,46 @@ struct part {
#endif
/* Store density/force specific stuff. */
union {
// union {
struct {
struct {
/* Particle velocity divergence. */
float div_v;
/* Particle velocity divergence. */
float div_v;
/* Derivative of particle number density. */
float wcount_dh;
/* Derivative of particle number density. */
float wcount_dh;
/* Particle velocity curl. */
float curl_v[3];
/* Particle velocity curl. */
float curl_v[3];
/* Particle number density. */
float wcount;
/* Particle number density. */
float wcount;
} density;
} density;
struct {
struct {
/* Balsara switch */
float balsara;
/* Balsara switch */
float balsara;
/* Aggregate quantities. */
float POrho2;
/* Aggregate quantities. */
float POrho2;
/* Change in particle energy over time. */
float u_dt;
/* Change in particle energy over time. */
float u_dt;
/* Change in smoothing length over time. */
float h_dt;
/* Change in smoothing length over time. */
float h_dt;
/* Signal velocity */
float v_sig;
/* Signal velocity */
float v_sig;
/* Sound speed */
float c;
/* Sound speed */
float c;
} force;
};
} force;
//};
/* Particle pressure. */
// float P;
......
......@@ -518,8 +518,6 @@ void runner_doinit(struct runner *r, struct cell *c) {
/* Get a direct pointer on the part. */
p = &parts[i];
if (p->id == 0) message("init!");
if (p->t_end <= t_end) {
/* Get ready for a density calculation */
......@@ -593,8 +591,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* Final operation on the density. */
p->rho = rho = ih * ih2 * (p->rho + p->mass * kernel_root);
p->rho_dh = rho_dh = (p->rho_dh - 3.0f * p->mass * kernel_root) * ih4;
wcount = (p->density.wcount + kernel_root) *
(4.0f / 3.0 * M_PI * kernel_gamma3);
p->density.wcount = wcount = (p->density.wcount + kernel_root) *
(4.0f / 3.0 * M_PI * kernel_gamma3);
wcount_dh =
p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
......@@ -610,20 +608,12 @@ void runner_doghost(struct runner *r, struct cell *c) {
h_corr = fmaxf(h_corr, -h / 2.f);
}
if (p->id == 0)
message("ghost! h=%f dh=%f wcount=%f rho=%f", p->h, h_corr, wcount,
p->rho);
/* Apply the correction to p->h and to the compact part. */
h = p->h += h_corr;
/* Did we get the right number density? */
if (wcount > kernel_nwneigh + const_delta_nwneigh ||
wcount < kernel_nwneigh - const_delta_nwneigh) {
// message( "particle %lli (h=%e,depth=%i) has bad wcount=%.3f." ,
// p->id , p->h , c->depth , wcount ); fflush(stdout);
// p->h += ( p->density.wcount + kernel_root - kernel_nwneigh ) /
// p->density.wcount_dh;
h = p->h += h_corr;
p->h += (p->density.wcount + kernel_root - kernel_nwneigh) /
p->density.wcount_dh;
pid[redo] = pid[i];
redo += 1;
p->density.wcount = 0.0;
......@@ -672,7 +662,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
(const_viscosity_alpha_max - p->alpha) * S;
/* Update particle's viscosity paramter */
p->alpha += alpha_dot * 1; // p->dt;
p->alpha += alpha_dot * (p->t_end - p->t_begin);
#endif
/* Reset the acceleration. */
......@@ -768,8 +758,6 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
p = &parts[k];
xp = &xparts[k];
if (p->id == 0) message("drift !");
/* Get local copies of particle data. */
h = p->h;
ih = 1.0f / h;
......@@ -816,11 +804,9 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
/* Predict gradient term */
p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (rho * xp->omega);
}
}
if (timer) {
#ifdef TIMER_VERBOSE
message("runner %02i: %i parts at depth %i took %.3f ms.", r->id, c->count,
......@@ -830,11 +816,8 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
TIMER_TOC(timer_drift);
#endif
}
}
/**
* @brief Combined second and first kick for fixed dt.
*
......@@ -879,15 +862,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
m = p->mass;
x[0] = p->x[0], x[1] = p->x[1], x[2] = p->x[2];
if (p->id == 0)
message("Kick ! t_beg=%f t_end=%f t_cur=%f", p->t_begin, p->t_end,
t_current);
/* if (p->id == 0) */
/* message("Kick ! t_beg=%f t_end=%f t_cur=%f", p->t_begin, p->t_end, */
/* t_current); */
/* If particle needs to be kicked */
if (p->t_end <= t_current) {
if (p->id == 0) message("Kicking like a boss");
/* First, finish the force loop */
p->force.h_dt *= p->h * 0.333333333f;
......@@ -919,16 +900,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
t_end = p->t_end + 0.5f * new_dt;
dt = t_end - t_start;
if (p->id == 0)
message("Kick ! dt=%f t_beg=%f t_end=%f", dt, p->t_begin, p->t_end);
/* Move particle forward in time */
p->t_begin = p->t_end;
p->t_end = p->t_begin + dt;
if (p->id == 0)
message("Kick ! dt=%f t_beg=%f t_end=%f", dt, p->t_begin, p->t_end);
/* Kick particles in momentum space */
xp->v_full[0] += p->a[0] * dt;
xp->v_full[1] += p->a[1] * dt;
......
......@@ -314,7 +314,8 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
u, COMPULSORY);
readArray(h_grp, "ParticleIDs", ULONGLONG, *N, 1, *parts, N_total, offset,
id, COMPULSORY);
/* readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, N_total, offset, dt, */
/* readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, N_total, offset, dt,
*/
/* OPTIONAL); */
readArray(h_grp, "Acceleration", FLOAT, *N, 3, *parts, N_total, offset, a,
OPTIONAL);
......@@ -487,7 +488,8 @@ void writeArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
* @param type The #DATA_TYPE of the array.
* @param N The number of particles to write.
* @param dim The dimension of the data (1 for scalar, 3 for vector)
* @param part A (char*) pointer on the first occurrence of the field of interest
* @param part A (char*) pointer on the first occurrence of the field of
*interest
*in the parts array
* @param field The name (code name) of the field to read from.
* @param us The UnitSystem currently in use
......@@ -627,7 +629,8 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
prepareArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N_total, 1,
us, UNIT_CONV_NO_UNITS);
/* prepareArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N_total, 1, us, */
/* prepareArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N_total, 1, us,
*/
/* UNIT_CONV_TIME); */
prepareArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N_total, 3,
us, UNIT_CONV_ACCELERATION);
......
......@@ -353,7 +353,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
* @param type The #DATA_TYPE of the array.
* @param N The number of particles to write.
* @param dim The dimension of the data (1 for scalar, 3 for vector)
* @param part A (char*) pointer on the first occurrence of the field of interest
* @param part A (char*) pointer on the first occurrence of the field of
*interest
*in the parts array
* @param field The name (code name) of the field to read from.
* @param us The UnitSystem currently in use
......
......@@ -43,9 +43,9 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
"none", "sort", "self", "pair", "sub",
"ghost", "kick1", "kick2", "send", "recv",
"link", "grav_pp", "grav_mm", "grav_up", "grav_down"};
"none", "sort", "self", "pair", "sub", "init",
"ghost", "drift", "kick", "send", "recv", "link",
"grav_pp", "grav_mm", "grav_up", "grav_down"};
/**
* @brief Unlock the cell held by this task.
......
......@@ -19,7 +19,6 @@
*
******************************************************************************/
#include <math.h>
#include <stdlib.h>
#include <stddef.h>
......@@ -31,7 +30,6 @@
#include "tools.h"
#include "swift.h"
/**
* Factorize a given integer, attempts to keep larger pair of factors.
*/
......@@ -49,8 +47,6 @@ void factor(int value, int *f1, int *f2) {
}
}
/**
* @brief Compute the average number of pairs per particle using
* a brute-force O(N^2) computation.
......@@ -234,7 +230,6 @@ void pairs_single_grav(double *dim, long long int pid,
pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
}
/**
* @brief Test the density function by dumping it for two random parts.
*
......
......@@ -24,9 +24,7 @@ void density_dump(int N);
void pairs_single_grav(double *dim, long long int pid,
struct gpart *__restrict__ parts, int N, int periodic);
void pairs_single_density(double *dim, long long int pid,
struct part *__restrict__ parts, int N,
int periodic);
struct part *__restrict__ parts, int N, int periodic);
void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
int periodic);
......@@ -20,10 +20,10 @@ AM_CFLAGS = -I../src -DCPU_TPS=2.67e9 $(HDF5_CPPFLAGS)
AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# List of programs and scripts to run in the test suite
TESTS = testGreetings testReading.sh testSingle testTimeIntegration
TESTS = testGreetings testReading.sh testSingle testTimeIntegration testSPHStep
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration testSPHStep
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
......@@ -32,7 +32,8 @@ testReading_SOURCES = testReading.c
testTimeIntegration_SOURCES = testTimeIntegration.c
# Sources for test_single
testSPHStep_SOURCES = testSPHStep.c
testSingle_SOURCES = testSingle.c
testSingle_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
testSingle_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (C) 2015 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/>.
*
******************************************************************************/