Commit 3add41b7 authored by Johannes Blaschke's avatar Johannes Blaschke
Browse files

implement simulation resuming with tryps

parent f456dd6d
......@@ -17,11 +17,11 @@ if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Cray")
message("-- cray compiler detected")
if("${CMAKE_BUILD_TYPE}" STREQUAL Debug)
message("-- building Debug Version")
message("----> building Debug Version")
add_definitions(-DDEBUG)
set(CMAKE_C_FLAGS "-g")
else()
message("-- building Release Version")
message("----> building Release Version")
add_definitions(-DNDEBUG)
# aggress, autoprefetch, fusion2, intrinsics, list=a, noomp_trance, display_opt
set(CMAKE_C_FLAGS "-O3 -hfp3 -h c99,pl=./compiler_information,wp")
......
......@@ -413,7 +413,7 @@ int read_in_squirmer_data(squirmer_state* squirmers)
}
char buffer[BUFFERLENGTH];
while (fgets (buffer, BUFFERLENGTH, fp) != NULL)
while(fgets(buffer, BUFFERLENGTH, fp) != NULL)
{
squirmer s;
int sscanf_res = sscanf(buffer, "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
......@@ -604,7 +604,7 @@ initialize_world(world_state * world, int numtasks)
if(config.trypanosome_on)
{
if(init_trypanosome(world, tconfig.trypanosome_head_initpos))//Remember: this returns 1 if CRASH
if(init_trypanosome(world, tconfig.trypanosome_head_initpos)) //Remember: this returns 1 if CRASH
return 1;
}
......@@ -615,10 +615,13 @@ initialize_world(world_state * world, int numtasks)
printf("Continuing simulation \n");
}
if(read_in_particle_data(&world->part)) return EXIT_FAILURE;
if(config.squirmeron)
{
if(config.squirmeron == 1){
if(read_in_squirmer_data(&world->squirmers)) return EXIT_FAILURE;
}
if(config.trypanosome_on == 1){
tryp_FILE tryp_fp = tryp_support__access_fp("end load");
if(tryp_support__init_load(world, tryp_fp)) return EXIT_FAILURE;
}
if(read_in_conserved_quantities_and_check_for_conservation(get_particles(world), get_nparticles(world),get_squirmers(world), get_nsquirmers(world), &world->timestep))
return EXIT_FAILURE;
......@@ -669,8 +672,7 @@ initialize_world(world_state * world, int numtasks)
return 0;
}
int initialize_everything(int numtasks, world_state * world, char * simnumber)
{
int initialize_everything(int numtasks, world_state * world, char * simnumber){
START(initialization);
int numberprocs[2] = { 0, 0 };
......@@ -678,7 +680,7 @@ int initialize_everything(int numtasks, world_state * world, char * simnumber)
return 1;
optimize_parallelization(numberprocs, numtasks);
set_nprocs(&config.nprocs, numberprocs);
set_nprocs(& config.nprocs, numberprocs);
MPI_Cart_create(MPI_COMM_WORLD, 2, numberprocs, config.periods, 1, &cartcomm);
initialize_mpi_structs();
......@@ -700,9 +702,10 @@ int initialize_everything(int numtasks, world_state * world, char * simnumber)
if(config.ppercell > 0)
assert(cpu.nparticles > 0);
#ifdef DEBUG
if(config.debug) initialize_debug();
#endif
#ifdef DEBUG
if(config.debug) initialize_debug();
#endif
initialize_localtemperature();
if(config.pflow == 1)
initialize_slices();
......
......@@ -349,8 +349,8 @@ void run_n_timesteps(int n, world_state* world, int ext_control)
if(cpu.rank==0) printf("Beginning time loop, begin timestep = %d end timestep = %d \n", world->timestep, n);
//all_cpus_print_active_and_inactive_squirmers(get_squirmers(world), get_nsquirmers(world), "Squirmers[beginning]");
if(config.trypanosome_on) //TODO: this should be abstracted into object
init_tryp_fp();
// if(config.trypanosome_on) //TODO: this should be abstracted into object
// init_tryp_fp();
for (; world->timestep<=n; world->timestep++)
{
......@@ -385,30 +385,26 @@ void run_n_timesteps(int n, world_state* world, int ext_control)
fflush(0);
if((cpu.rank == 0) && (config.trypanosome_on) && (world->timestep%10 == 0))
{
FILE * tryp_fp_r = access_tryp_pos_fp();
FILE * tryp_fp_v = access_tryp_vel_fp();
FILE * tryp_fp_f = access_tryp_f_fp();
FILE * tryp_fp_flagellum = access_tryp_flagellum_fp();
print_trypanosome_to_file(world->trypanosome, tryp_fp_r, tryp_fp_v, tryp_fp_f, tryp_fp_flagellum , world->timestep);//TODO: implement averaging
if((cpu.rank == 0) && (config.trypanosome_on == 1) && (world->timestep % config.naverage == 0)){
tryp_FILE tryp_fp = tryp_support__access_fp("sim data");
print_trypanosome_to_file(world->trypanosome, tryp_fp, world->timestep); //TODO: implement averaging
}
if(world->timestep>0)
{
if(world->timestep % (500000) == 0 || last)
{
if(cpu.rank==0)
{
printf("\n Printing world state to file \n");
}
if(world->timestep>0){
if((world->timestep % 500000 == 0 ) || last){
if(cpu.rank == 0) printf("\n Printing world state to file \n");
print_particle_enddata_to_file(get_particles(world), get_nparticles(world));
if(config.squirmeron)
{
if(config.squirmeron == 1)
print_squirmer_enddata_to_file(get_squirmers(world), get_nsquirmers(world));
}
print_conserved_quantities_to_file(get_particles(world), get_nparticles(world),get_squirmers(world), get_nsquirmers(world), world->timestep);
if((cpu.rank == 0) && (config.trypanosome_on == 1)){
tryp_FILE tryp_fp = tryp_support__access_fp("end data");
print_trypanosome_to_file(world->trypanosome, tryp_fp, -1);
}
}
}
}
......@@ -469,22 +465,21 @@ void main_support__print_core_domains(FILE * file_dest){
free(rec_max_z);
}
int main(int argc, char * argv[])
{
signal(SIGSEGV, &printBacktrace);
signal(SIGABRT, &printBacktrace);
int main(int argc, char * argv[]){
signal(SIGSEGV, & printBacktrace);
signal(SIGABRT, & printBacktrace);
int numtasks;
world_state world;
mpcc_rngzig_init(); // initialize random number generator
MPI_Init(&argc, &argv);
MPI_Init(& argc, & argv);
initialize_time_measure();
START(entire_simulation);
MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
MPI_Comm_size(MPI_COMM_WORLD, & numtasks);
int rank = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_rank(MPI_COMM_WORLD, & rank);
if(rank == 0){
fprintf(stdout, "INIT: MPI-Comms initialized:\n");
......
......@@ -295,6 +295,8 @@ obstacle_time_driven_md(world_state * world, int n_md_step, const double t, cons
// fluid particle streaming ************************************************************************
stream_particles(world, dt, i, last);
// this has to be done again as tryp particle v might have changed....
assign_trypanosome_nodes_to_particles(tryp_cloud, world->trypanosome, cpu.gshift, world);
t_remaining -= dt;
}
......
......@@ -1125,93 +1125,328 @@ trypanosome_is_in_core_domain(trypanosome * tryp, vector3 gshift)
return tryp_in_domain;
}
void
init_tryp_fp()
{
char path[100] = "tryp_properties.dat";
FILE * fp = fopen( path, "w" );
fprintf(fp, "n_step: (r_x,r_y,r_z,v_x,v_y,v_z,f_x,f_y,f_z) ... Naf times \n");
fflush(fp);
fclose(fp);
FILE * access_tryp_pos_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_r_properties.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_r_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_r_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
FILE *
access_tryp_pos_fp()
{
char path[100] = "tryp_r_properties.dat";
return fopen(path, "a" );
FILE * access_tryp_vel_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_v_properties.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_v_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_v_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
FILE *
access_tryp_vel_fp()
{
char path[100] = "tryp_v_properties.dat";
return fopen(path, "a" );
FILE * access_tryp_f_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_f_properties.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_f_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_f_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
FILE *
access_tryp_f_fp()
{
char path[100] = "tryp_f_properties.dat";
return fopen(path, "a" );
FILE * access_tryp_f_ext_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_f_ext_properties.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_f_ext_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_f_ext_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
FILE *
access_tryp_flagellum_fp()
{
char path[100] = "tryp_r_flagellum.dat";
return fopen(path, "a" );
FILE * access_tryp_flagellum_r_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_r_flagellum.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_r_flagellum_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_r_flagellum_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
FILE * access_tryp_flagellum_v_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_v_flagellum.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_v_flagellum_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_v_flagellum_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
FILE * access_tryp_flagellum_f_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_f_flagellum.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_f_flagellum_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_f_flagellum_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
FILE * access_tryp_flagellum_f_ext_fp(const char * mode){
if(strcmp(mode, "sim data") == 0){
char * path = "tryp_f_ext_flagellum.dat";
return fopen(path, "a" );
} else if(strcmp(mode, "end data") == 0){
char * path = "tryp_f_ext_flagellum_enddata.dat";
return fopen(path, "w" );
} else if(strcmp(mode, "end load") == 0){
char * path = "tryp_f_ext_flagellum_enddata.dat";
return fopen(path, "r");
}
return NULL;
}
void
print_trypanosome_to_file(trypanosome * tryp, FILE * fp_r, FILE * fp_v, FILE* fp_f, FILE * fp_flag, int timestep)
{
fprintf(fp_r,"%d: ", timestep);
fprintf(fp_v,"%d: ", timestep);
fprintf(fp_f,"%d: ", timestep);
fprintf(fp_flag,"%d: ", timestep);
void tryp_support__close_fp(tryp_FILE tryp_fp){
fprintf(tryp_fp.fp_r, "\n");
fflush(tryp_fp.fp_r);
fclose(tryp_fp.fp_r);
fprintf(tryp_fp.fp_v, "\n");
fflush(tryp_fp.fp_v);
fclose(tryp_fp.fp_v);
fprintf(tryp_fp.fp_f, "\n");
fflush(tryp_fp.fp_f);
fclose(tryp_fp.fp_f);
fprintf(tryp_fp.fp_f_ext, "\n");
fflush(tryp_fp.fp_f_ext);
fclose(tryp_fp.fp_f_ext);
char node_format[] = "(%f,%f,%f)";
fprintf(tryp_fp.fp_flag_r, "\n");
fflush(tryp_fp.fp_flag_r);
fclose(tryp_fp.fp_flag_r);
fprintf(tryp_fp.fp_flag_v, "\n");
fflush(tryp_fp.fp_flag_v);
fclose(tryp_fp.fp_flag_v);
fprintf(tryp_fp.fp_flag_f, "\n");
fflush(tryp_fp.fp_flag_f);
fclose(tryp_fp.fp_flag_f);
fprintf(tryp_fp.fp_flag_f_ext, "\n");
fflush(tryp_fp.fp_flag_f_ext);
fclose(tryp_fp.fp_flag_f_ext);
}
tryp_FILE tryp_support__access_fp(const char * mode){
tryp_FILE tryp_fp;
tryp_fp.fp_r = access_tryp_pos_fp(mode);
tryp_fp.fp_v = access_tryp_vel_fp(mode);
tryp_fp.fp_f = access_tryp_f_fp(mode);
tryp_fp.fp_f_ext = access_tryp_f_ext_fp(mode);
tryp_fp.fp_flag_r = access_tryp_flagellum_r_fp(mode);
tryp_fp.fp_flag_v = access_tryp_flagellum_v_fp(mode);
tryp_fp.fp_flag_f = access_tryp_flagellum_f_fp(mode);
tryp_fp.fp_flag_f_ext = access_tryp_flagellum_f_ext_fp(mode);
return tryp_fp;
}
void print_trypanosome_to_file(trypanosome * tryp, tryp_FILE tryp_fp, int timestep){
if(timestep > -1){
fprintf(tryp_fp.fp_r, "%d: ", timestep);
fprintf(tryp_fp.fp_v, "%d: ", timestep);
fprintf(tryp_fp.fp_f, "%d: ", timestep);
fprintf(tryp_fp.fp_f_ext, "%d: ", timestep);
fprintf(tryp_fp.fp_flag_r, "%d: ", timestep);
fprintf(tryp_fp.fp_flag_v, "%d: ", timestep);
fprintf(tryp_fp.fp_flag_f, "%d: ", timestep);
fprintf(tryp_fp.fp_flag_f_ext, "%d: ", timestep);
}
char node_format[] = "(%lf,%lf,%lf)";
vertex * vertices = tryp->tryp_net.skeleton.vertices;
int n_vertices = tryp->tryp_net.skeleton.n_vertex;
int n_vertices = tryp->tryp_net.skeleton.n_vertex;
curve * flagellum = &(tryp->tryp_net.flagellum);
int n_flagellum = tryp->tryp_net.flagellum.length;
int n_flagellum = tryp->tryp_net.flagellum.length;
for(int i=0; i<n_vertices; i++)
{
vector3 r = vertices[i].r;
vector3 v = vertices[i].v;
vector3 f = vertices[i].f;
for(int i = 0; i < n_vertices; i++){
vector3 r = vertices[i].r;
vector3 v = vertices[i].v;
vector3 f = vertices[i].f;
vector3 f_ext = vertices[i].f_ext;
vector3 f_total;
vadd(&f_total, &f_ext, &f);
fprintf(fp_r, node_format, r.x, r.y, r.z);
fprintf(fp_v, node_format, v.x, v.y, v.z);
fprintf(fp_f, node_format, f_total.x, f_total.y, f_total.z);
fprintf(tryp_fp.fp_r, node_format, r.x, r.y, r.z);
fprintf(tryp_fp.fp_v, node_format, v.x, v.y, v.z);
fprintf(tryp_fp.fp_f, node_format, f.x, f.y, f.z);
fprintf(tryp_fp.fp_f_ext, node_format, f_ext.x, f_ext.y, f_ext.z);
}
for(int i=0; i<n_flagellum; i++)
{
for(int i = 0; i < n_flagellum; i++){
vector3 r = flagellum->edges[i].start->r;
fprintf(fp_flag, "(%f,%f,%f)", r.x, r.y, r.z);
fprintf(tryp_fp.fp_flag_r, node_format, r.x, r.y, r.z);
vector3 v = flagellum->edges[i].start->v;
fprintf(tryp_fp.fp_flag_v, node_format, v.x, v.y, v.z);
vector3 f = flagellum->edges[i].start->f;
fprintf(tryp_fp.fp_flag_f, node_format, f.x, f.y, f.z);
vector3 f_ext = flagellum->edges[i].start->f_ext;
fprintf(tryp_fp.fp_flag_f_ext, node_format, f_ext.x, f_ext.y, f_ext.z);
}
vector3 r = flagellum->edges[flagellum->length-1].end->r;
fprintf(fp_flag, "(%f,%f,%f)", r.x, r.y, r.z);
int i_end = flagellum->length-1;
vector3 r = flagellum->edges[i_end].end->r;
fprintf(tryp_fp.fp_flag_r, node_format, r.x, r.y, r.z);
vector3 v = flagellum->edges[i_end].start->v;
fprintf(tryp_fp.fp_flag_v, node_format, v.x, v.y, v.z);
vector3 f = flagellum->edges[i_end].start->f;
fprintf(tryp_fp.fp_flag_f, node_format, f.x, f.y, f.z);
fprintf(fp_r, "\n");
fflush(fp_r);
fclose(fp_r);
vector3 f_ext = flagellum->edges[i_end].start->f_ext;
fprintf(tryp_fp.fp_flag_f_ext, node_format, f_ext.x, f_ext.y, f_ext.z);
fprintf(fp_v, "\n");
fflush(fp_v);
fclose(fp_v);
tryp_support__close_fp(tryp_fp);
}
int tryp_support__init_load(world_state * world, tryp_FILE tryp_fp){
char node_format[] = "(%lf,%lf,%lf)";
trypanosome * tryp = world->trypanosome;
vertex * vertices = tryp->tryp_net.skeleton.vertices;
int n_vertices = tryp->tryp_net.skeleton.n_vertex;
curve * flagellum = & tryp->tryp_net.flagellum;
int n_flagellum = tryp->tryp_net.flagellum.length;
if(cpu.rank == 0){
for(int i = 0; i < n_vertices; i++){
vector3 * r = & vertices[i].r;
vector3 * v = & vertices[i].v;
vector3 * f = & vertices[i].f;
vector3 * f_ext = & vertices[i].f_ext;
if(fscanf(tryp_fp.fp_r, node_format, & r->x, & r->y, & r->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
if(fscanf(tryp_fp.fp_v, node_format, & v->x, & v->y, & v->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
if(fscanf(tryp_fp.fp_f, node_format, & f->x, & f->y, & f->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
if(fscanf(tryp_fp.fp_f_ext, node_format, & f_ext->x, & f_ext->y, & f_ext->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
}
for(int i = 0; i < n_flagellum; i++){
vector3 * r = & flagellum->edges[i].start->r;
fprintf(fp_f, "\n");
fflush(fp_f);
fclose(fp_f);
if(fscanf(tryp_fp.fp_flag_r, node_format, & r->x, & r->y, & r->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
vector3 * v = & flagellum->edges[i].start->v;
if(fscanf(tryp_fp.fp_flag_v, node_format, & v->x, & v->y, & v->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
vector3 * f = & flagellum->edges[i].start->f;
if(fscanf(tryp_fp.fp_flag_f, node_format, & f->x, & f->y, & f->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
vector3 * f_ext = & flagellum->edges[i].start->f_ext;
if(fscanf(tryp_fp.fp_flag_f_ext, node_format, & f_ext->x, & f_ext->y, & f_ext->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
}
int i_end = flagellum->length-1;
fprintf(fp_flag, "\n");
fflush(fp_flag);
fclose(fp_flag);
vector3 * r = & flagellum->edges[i_end].end->r;
if(fscanf(tryp_fp.fp_flag_r, node_format, & r->x, & r->y, & r->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
vector3 * v = & flagellum->edges[i_end].start->v;
if(fscanf(tryp_fp.fp_flag_v, node_format, & v->x, & v->y, & v->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
vector3 * f = & flagellum->edges[i_end].start->f;
if(fscanf(tryp_fp.fp_flag_f, node_format, & f->x, & f->y, & f->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
vector3 * f_ext = & flagellum->edges[i_end].start->f_ext;
if(fscanf(tryp_fp.fp_flag_f_ext, node_format, & f_ext->x, & f_ext->y, & f_ext->z) != 3){
tryp_support__close_fp(tryp_fp);
return 1;
}
tryp_support__close_fp(tryp_fp);
trypanosome_network__update_metadata(& tryp->tryp_net);
}
MPI_Bcast_trypanosome_r(world->trypanosome);
assign_trypanosome_nodes_to_particles(tryp->tryp_cloud, tryp, cpu.gshift, world);
return 0;
}
double wca_force(double r, double e, double d)
{
return 4*e*( 0.25*pow( (d/r) ,12) - 0.5*pow( (d/r) ,6))+e;
......
......@@ -35,12 +35,29 @@ void tryp_cloud_to_del_p(particle * tryp_cloud, trypanosome * tryp);
int trypanosome_is_in_core_domain(trypanosome * tryp, vector3 gshift);
void init_tryp_fp(void);
FILE * access_tryp_pos_fp(void);
FILE * access_tryp_vel_fp(void);
FILE * access_tryp_f_fp(void);
FILE * access_tryp_flagellum_fp(void);
void print_trypanosome_to_file(trypanosome * tryp, FILE * fp_r, FILE * fp_v, FILE* fp_f, FILE * fp_flag, int timestep);
typedef struct {
FILE * fp_r;
FILE * fp_v;
FILE * fp_f;
FILE * fp_f_ext;
FILE * fp_flag_r;
FILE * fp_flag_v;
FILE * fp_flag_f;
FILE * fp_flag_f_ext;
} tryp_FILE;
FILE * access_tryp_pos_fp(const char * mode);
FILE * access_tryp_vel_fp(const char * mode);
FILE * access_tryp_f_fp(const char * mode);
FILE * access_tryp_f_ext_fp(const char * mode);
FILE * access_tryp_flagellum_r_fp(const char * mode);
FILE * access_tryp_flagellum_v_fp(const char * mode);
FILE * access_tryp_flagellum_f_fp(const char * mode);
FILE * access_tryp_flagellum_f_ext_fp(const char * mode);
tryp_FILE tryp_support__access_fp(const char * mode);
void print_trypanosome_to_file(trypanosome * tryp, tryp_FILE tryp_fp , int timestep);
int tryp_support__init_load(world_state * world, tryp_FILE tryp_fp);