Skip to content
Snippets Groups Projects
Commit 2d219a25 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

merging mast branch back in.

parents ba05aa17 9f1fd7ec
No related branches found
No related tags found
1 merge request!50Parallel sorting
Showing
with 84 additions and 569 deletions
...@@ -18,13 +18,17 @@ doc/html/ ...@@ -18,13 +18,17 @@ doc/html/
doc/latex/ doc/latex/
doc/man/ doc/man/
doc/Doxyfile doc/Doxyfile
examples/test examples/swift
examples/test_fixdt examples/swift_fixdt
examples/test_fixdt_mpi examples/swift_fixdt_mpi
examples/test_mindt examples/swift_mindt
examples/test_mindt_mpi examples/swift_mindt_mpi
examples/test_mpi examples/swift_mpi
examples/test_single
tests/testGreetings
tests/testReading
tests/input.hdf5
tests/testSingle
m4/libtool.m4 m4/libtool.m4
m4/ltoptions.m4 m4/ltoptions.m4
......
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 pedro.gonnet@durham.ac.uk. # Copyright (c) 2012 pedro.gonnet@durham.ac.uk.
# 2015 matthieu.schaller@durham.ac.uk. # 2015 matthieu.schaller@durham.ac.uk.
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Copyright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # 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 # it under the terms of the GNU Lesser General Public License as published
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
......
#!/bin/bash #!/bin/bash
wget http://www.dur.ac.uk/pedro.gonnet/cosmoVolume.hdf5 wget http://icc.dur.ac.uk/~jlvc76/Files/SWIFT/cosmoVolume.hdf5
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
......
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk). # Matthieu Schaller (matthieu.schaller@durham.ac.uk).
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
...@@ -29,45 +29,40 @@ MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS) ...@@ -29,45 +29,40 @@ MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
MPI_FLAGS = -DWITH_MPI $(METIS_INCS) MPI_FLAGS = -DWITH_MPI $(METIS_INCS)
# Set-up the library # Set-up the library
bin_PROGRAMS = test test_fixdt test_mindt test_single bin_PROGRAMS = swift swift_fixdt swift_mindt
# Build MPI versions as well? # Build MPI versions as well?
if HAVEMPI if HAVEMPI
bin_PROGRAMS += test_mpi test_fixdt_mpi test_mindt_mpi bin_PROGRAMS += swift_mpi swift_fixdt_mpi swift_mindt_mpi
endif endif
# Sources for test # Sources for swift
test_SOURCES = test.c swift_SOURCES = main.c
test_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep | engine_policy_setaffinity" swift_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep | engine_policy_setaffinity"
test_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) swift_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# Sources for test_fixdt # Sources for swift_fixdt
test_fixdt_SOURCES = test.c swift_fixdt_SOURCES = main.c
test_fixdt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep | engine_policy_setaffinity" swift_fixdt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep | engine_policy_setaffinity"
test_fixdt_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) swift_fixdt_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# Sources for test_mindt # Sources for swift_mindt
test_mindt_SOURCES = test.c swift_mindt_SOURCES = main.c
test_mindt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_keep | engine_policy_setaffinity" swift_mindt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_keep | engine_policy_setaffinity"
test_mindt_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) swift_mindt_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# Sources for test_mpi # Sources for swift_mpi
test_mpi_SOURCES = test.c swift_mpi_SOURCES = main.c
test_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep" swift_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_multistep | engine_policy_keep"
test_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS) swift_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
# Sources for test_fixdt_mpi # Sources for swift_fixdt_mpi
test_fixdt_mpi_SOURCES = test.c swift_fixdt_mpi_SOURCES = main.c
test_fixdt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep" swift_fixdt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep"
test_fixdt_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS) swift_fixdt_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
# Sources for test_mindt_mpi # Sources for swift_mindt_mpi
test_mindt_mpi_SOURCES = test.c swift_mindt_mpi_SOURCES = main.c
test_mindt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_keep" swift_mindt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_keep"
test_mindt_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS) swift_mindt_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
# Sources for test_single
test_single_SOURCES = test_single.c
test_single_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
test_single_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
...@@ -63,7 +63,7 @@ for i in range(L): ...@@ -63,7 +63,7 @@ for i in range(L):
v[index,1] = 0. v[index,1] = 0.
v[index,2] = 0. v[index,2] = 0.
m[index] = mass m[index] = mass
h[index] = 2.251 * boxSize / L h[index] = 1.1255 * boxSize / L
u[index] = internalEnergy u[index] = internalEnergy
ids[index] = index ids[index] = index
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
...@@ -67,7 +67,7 @@ for i in range(L): ...@@ -67,7 +67,7 @@ for i in range(L):
v[index,1] = 0. v[index,1] = 0.
v[index,2] = 0. v[index,2] = 0.
m[index] = mass m[index] = mass
h[index] = 2.251 / 2 * boxSize / L h[index] = 1.1255 * boxSize / L
u[index] = internalEnergy u[index] = internalEnergy
ids[index] = index ids[index] = index
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L: if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
...@@ -70,7 +70,7 @@ for i in range(L): ...@@ -70,7 +70,7 @@ for i in range(L):
v[index,1] = 0. v[index,1] = 0.
v[index,2] = 0. v[index,2] = 0.
m[index] = mass m[index] = mass
h[index] = 2.251 / 2 * hbox h[index] = 1.1255 * hbox
u[index] = internalEnergy u[index] = internalEnergy
ids[index] = index ids[index] = index
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox: if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
...@@ -29,7 +29,7 @@ from numpy import * ...@@ -29,7 +29,7 @@ from numpy import *
# Input values? # Input values?
if len(sys.argv) < 3 : if len(sys.argv) < 3 :
print "Useage: " , sys.argv[0] , " <filename> <nr. bins>" print "Usage: " , sys.argv[0] , " <filename> <nr. bins>"
exit() exit()
# Get the input arguments # Get the input arguments
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
...@@ -29,7 +29,7 @@ from numpy import * ...@@ -29,7 +29,7 @@ from numpy import *
# Input values? # Input values?
if len(sys.argv) < 3 : if len(sys.argv) < 3 :
print "Useage: " , sys.argv[0] , " <filename> <nr. bins>" print "Usage: " , sys.argv[0] , " <filename> <nr. bins>"
exit() exit()
# Get the input arguments # Get the input arguments
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
# This program is free software: you can redistribute it and/or modify # This program is free software: you can redistribute it and/or modify
...@@ -74,7 +74,7 @@ ds = grp.create_dataset('Masses', (numPart,1), 'f') ...@@ -74,7 +74,7 @@ ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m ds[()] = m
m = zeros(1) m = zeros(1)
h = full((numPart, 1), 2.251 * boxSize / L) h = full((numPart, 1), 1.1255 * boxSize / L)
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f') ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h ds[()] = h
h = zeros(1) h = zeros(1)
......
...@@ -56,493 +56,8 @@ ...@@ -56,493 +56,8 @@
#define ENGINE_POLICY engine_policy_none #define ENGINE_POLICY engine_policy_none
#endif #endif
/**
* @brief Mapping function to draw a specific cell (gnuplot).
*/
void map_cells_plot(struct cell *c, void *data) {
int depth = *(int *)data;
double *l = c->loc, *h = c->h;
if (c->depth <= depth) {
printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2]);
printf("%.16e %.16e %.16e\n\n\n", l[0], l[1] + h[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n\n\n", l[0], l[1] + h[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n\n\n", l[0], l[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2]);
printf("%.16e %.16e %.16e\n", l[0], l[1] + h[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n", l[0] + h[0], l[1] + h[1], l[2] + h[2]);
printf("%.16e %.16e %.16e\n\n\n", l[0] + h[0], l[1] + h[1], l[2]);
if (!c->split) {
for (int k = 0; k < c->count; k++)
printf("0 0 0 %.16e %.16e %.16e\n", c->parts[k].x[0], c->parts[k].x[1],
c->parts[k].x[2]);
printf("\n\n");
}
/* else
for ( int k = 0 ; k < 8 ; k++ )
if ( c->progeny[k] != NULL )
map_cells_plot( c->progeny[k] , data ); */
}
}
/**
* @brief Mapping function for checking if each part is in its box.
*/
/* void map_check ( struct part *p , struct cell *c , void *data ) {
if ( p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] ||
p->x[0] < c->loc[0] || p->x[0] > c->loc[0]+c->h[0] )
printf( "map_check: particle %i is outside of its box.\n" , p->id );
} */
/**
* @brief Mapping function for neighbour count.
*/
void map_cellcheck(struct cell *c, void *data) {
int k, *count = (int *)data;
struct part *p;
__sync_fetch_and_add(count, c->count);
/* Loop over all parts and check if they are in the cell. */
for (k = 0; k < c->count; k++) {
p = &c->parts[k];
if (p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] ||
p->x[2] > c->loc[2] + c->h[2]) {
printf(
"map_cellcheck: particle at [ %.16e %.16e %.16e ] outside of cell [ "
"%.16e %.16e %.16e ] - [ %.16e %.16e %.16e ].\n",
p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
c->loc[0] + c->h[0], c->loc[1] + c->h[1], c->loc[2] + c->h[2]);
error("particle out of bounds!");
}
}
}
/**
* @brief Mapping function for maxdepth cell count.
*/
void map_maxdepth(struct cell *c, void *data) {
int maxdepth = ((int *)data)[0];
int *count = &((int *)data)[1];
// printf( "%e\n" , p->count );
if (c->depth == maxdepth) *count += 1;
}
/**
* @brief Mapping function for neighbour count.
*/
void map_count(struct part *p, struct cell *c, void *data) {
double *wcount = (double *)data;
// printf( "%i %e %e\n" , p->id , p->count , p->count_dh );
*wcount += p->density.wcount;
}
void map_wcount_min(struct part *p, struct cell *c, void *data) {
struct part **p2 = (struct part **)data;
if (p->density.wcount < (*p2)->density.wcount) *p2 = p;
}
void map_wcount_max(struct part *p, struct cell *c, void *data) {
struct part **p2 = (struct part **)data;
if (p->density.wcount > (*p2)->density.wcount) *p2 = p;
}
void map_h_min(struct part *p, struct cell *c, void *data) {
struct part **p2 = (struct part **)data;
if (p->h < (*p2)->h) *p2 = p;
}
void map_h_max(struct part *p, struct cell *c, void *data) {
struct part **p2 = (struct part **)data;
if (p->h > (*p2)->h) *p2 = p;
}
/**
* @brief Mapping function for neighbour count.
*/
void map_icount(struct part *p, struct cell *c, void *data) {
// int *count = (int *)data;
// printf( "%i\n" , p->icount );
// *count += p->icount;
}
/**
* @brief Mapping function to print the particle position.
*/
void map_dump(struct part *p, struct cell *c, void *data) {
double *shift = (double *)data;
printf("%g\t%g\t%g\n", p->x[0] - shift[0], p->x[1] - shift[1],
p->x[2] - shift[2]);
}
/**
* @brief Compute the average number of pairs per particle using
* a brute-force O(N^2) computation.
*
* @param dim The space dimensions.
* @param parts The #part array.
* @param N The number of parts.
* @param periodic Periodic boundary conditions flag.
*/
void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
int periodic) {
int i, j, k, count = 0;
// int mj, mk;
// double maxratio = 1.0;
double r2, dx[3], rho = 0.0;
double rho_max = 0.0, rho_min = 100;
/* Loop over all particle pairs. */
for (j = 0; j < N; j++) {
if (j % 1000 == 0) {
printf("pairs_n2: j=%i.\n", j);
fflush(stdout);
}
for (k = j + 1; k < N; k++) {
for (i = 0; i < 3; i++) {
dx[i] = parts[j].x[i] - parts[k].x[i];
if (periodic) {
if (dx[i] < -dim[i] / 2)
dx[i] += dim[i];
else if (dx[i] > dim[i] / 2)
dx[i] -= dim[i];
}
}
r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
if (r2 < parts[j].h * parts[j].h || r2 < parts[k].h * parts[k].h) {
runner_iact_density(r2, NULL, parts[j].h, parts[k].h, &parts[j],
&parts[k]);
/* if ( parts[j].h / parts[k].h > maxratio )
{
maxratio = parts[j].h / parts[k].h;
mj = j; mk = k;
}
else if ( parts[k].h / parts[j].h > maxratio )
{
maxratio = parts[k].h / parts[j].h;
mj = j; mk = k;
} */
}
}
}
/* Aggregate the results. */
for (k = 0; k < N; k++) {
// count += parts[k].icount;
rho += parts[k].density.wcount;
rho_min = fmin(parts[k].density.wcount, rho_min);
rho_min = fmax(parts[k].density.wcount, rho_max);
}
/* Dump the result. */
printf("pairs_n2: avg. density per part is %.3f (nr. pairs %.3f).\n",
rho / N + 32.0 / 3, ((double)count) / N);
printf("pairs_n2: densities are in [ %e , %e ].\n", rho_min / N + 32.0 / 3,
rho_max / N + 32.0 / 3);
/* printf( "pairs_n2: maximum ratio between parts %i [%e,%e,%e] and %i
[%e,%e,%e] is %.3f/%.3f\n" ,
mj , parts[mj].x[0] , parts[mj].x[1] , parts[mj].x[2] ,
mk , parts[mk].x[0] , parts[mk].x[1] , parts[mk].x[2] ,
parts[mj].h , parts[mk].h ); fflush(stdout); */
fflush(stdout);
}
void pairs_single_density(double *dim, long long int pid,
struct part *__restrict__ parts, int N,
int periodic) {
int i, k;
// int mj, mk;
// double maxratio = 1.0;
double r2, dx[3];
float fdx[3];
struct part p;
// double ih = 12.0/6.25;
/* Find "our" part. */
for (k = 0; k < N && parts[k].id != pid; k++)
;
if (k == N) error("Part not found.");
p = parts[k];
printf("pairs_single: part[%i].id == %lli.\n", k, pid);
p.rho = 0.0;
p.density.wcount = 0.0;
// p.icount = 0;
p.rho_dh = 0.0;
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
if (parts[k].id == p.id) continue;
for (i = 0; i < 3; i++) {
dx[i] = p.x[i] - parts[k].x[i];
if (periodic) {
if (dx[i] < -dim[i] / 2)
dx[i] += dim[i];
else if (dx[i] > dim[i] / 2)
dx[i] -= dim[i];
}
fdx[i] = dx[i];
}
r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
if (r2 < p.h * p.h) {
runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
/* printf( "pairs_simple: interacting particles %lli [%i,%i,%i] and %lli
[%i,%i,%i], r=%e.\n" ,
pid , (int)(p.x[0]*ih) , (int)(p.x[1]*ih) , (int)(p.x[2]*ih) ,
parts[k].id , (int)(parts[k].x[0]*ih) , (int)(parts[k].x[1]*ih) ,
(int)(parts[k].x[2]*ih) ,
sqrtf(r2) ); */
}
}
/* Dump the result. */
printf("pairs_single: wcount of part %lli (h=%e) is %f.\n", p.id, p.h,
p.density.wcount + 32.0 / 3);
fflush(stdout);
}
void pairs_single_grav(double *dim, long long int pid,
struct gpart *__restrict__ parts, int N, int periodic) {
int i, k;
// int mj, mk;
// double maxratio = 1.0;
double r2, dx[3];
float fdx[3], a[3] = {0.0, 0.0, 0.0}, aabs[3] = {0.0, 0.0, 0.0};
struct gpart pi, pj;
// double ih = 12.0/6.25;
/* Find "our" part. */
for (k = 0; k < N; k++)
if ((parts[k].id > 0 && parts[k].part->id == pid) || parts[k].id == -pid)
break;
if (k == N) error("Part not found.");
pi = parts[k];
pi.a[0] = 0.0f;
pi.a[1] = 0.0f;
pi.a[2] = 0.0f;
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
if (parts[k].id == pi.id) continue;
pj = parts[k];
for (i = 0; i < 3; i++) {
dx[i] = pi.x[i] - pj.x[i];
if (periodic) {
if (dx[i] < -dim[i] / 2)
dx[i] += dim[i];
else if (dx[i] > dim[i] / 2)
dx[i] -= dim[i];
}
fdx[i] = dx[i];
}
r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
runner_iact_grav(r2, fdx, &pi, &pj);
a[0] += pi.a[0];
a[1] += pi.a[1];
a[2] += pi.a[2];
aabs[0] += fabsf(pi.a[0]);
aabs[1] += fabsf(pi.a[1]);
aabs[2] += fabsf(pi.a[2]);
pi.a[0] = 0.0f;
pi.a[1] = 0.0f;
pi.a[2] = 0.0f;
}
/* Dump the result. */
message(
"acceleration on gpart %lli is a=[ %e %e %e ], |a|=[ %.2e %.2e %.2e ].\n",
pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
}
/**
* @brief Test the kernel function by dumping it in the interval [0,1].
*
* @param N number of intervals in [0,1].
*/
void kernel_dump(int N) {
int k;
float x, w, dw_dx;
float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
// float dw_dx4[4] __attribute__ ((aligned (16)));
for (k = 0; k <= N; k++) {
x = ((float)k) / N;
x4[3] = x4[2];
x4[2] = x4[1];
x4[1] = x4[0];
x4[0] = x;
kernel_deval(x, &w, &dw_dx);
// kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
printf(" %e %e %e %e %e %e %e\n", x, w, dw_dx, w4[0], w4[1], w4[2], w4[3]);
}
}
float gadget(float r) {
float fac, h_inv, u, r2 = r * r;
if (r >= const_epsilon)
fac = 1.0f / (r2 * r);
else {
h_inv = 1. / const_epsilon;
u = r * h_inv;
if (u < 0.5)
fac = const_iepsilon3 * (10.666666666667 + u * u * (32.0 * u - 38.4));
else
fac = const_iepsilon3 *
(21.333333333333 - 48.0 * u + 38.4 * u * u -
10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
}
return const_G * fac;
}
void gravity_dump(float r_max, int N) {
int k;
float x, w;
float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
// float dw_dx4[4] __attribute__ ((aligned (16)));
for (k = 1; k <= N; k++) {
x = (r_max * k) / N;
x4[3] = x4[2];
x4[2] = x4[1];
x4[1] = x4[0];
x4[0] = x;
kernel_grav_eval(x, &w);
w *= const_G / (x * x * x);
// blender_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
printf(" %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n", x, w * x, w4[0],
w4[1], w4[2], w4[3], gadget(x) * x);
}
}
/**
* @brief Test the density function by dumping it for two random parts.
*
* @param N number of intervals in [0,1].
*/
void density_dump(int N) {
int k;
float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
struct part *pi[4], *pj[4], Pi[4], Pj[4];
/* Init the interaction parameters. */
for (k = 0; k < 4; k++) {
Pi[k].mass = 1.0f;
Pi[k].rho = 0.0f;
Pi[k].density.wcount = 0.0f;
Pj[k].mass = 1.0f;
Pj[k].rho = 0.0f;
Pj[k].density.wcount = 0.0f;
hi[k] = 1.0;
hj[k] = 1.0;
pi[k] = &Pi[k];
pj[k] = &Pj[k];
}
for (k = 0; k <= N; k++) {
r2[3] = r2[2];
r2[2] = r2[1];
r2[1] = r2[0];
r2[0] = ((float)k) / N;
Pi[0].density.wcount = 0;
Pj[0].density.wcount = 0;
runner_iact_density(r2[0], NULL, hi[0], hj[0], &Pi[0], &Pj[0]);
printf(" %e %e %e", r2[0], Pi[0].density.wcount, Pj[0].density.wcount);
Pi[0].density.wcount = 0;
Pj[0].density.wcount = 0;
Pi[1].density.wcount = 0;
Pj[1].density.wcount = 0;
Pi[2].density.wcount = 0;
Pj[2].density.wcount = 0;
Pi[3].density.wcount = 0;
Pj[3].density.wcount = 0;
runner_iact_vec_density(r2, NULL, hi, hj, pi, pj);
printf(" %e %e %e %e\n", Pi[0].density.wcount, Pi[1].density.wcount,
Pi[2].density.wcount, Pi[3].density.wcount);
}
}
/**
* Factorize a given integer, attempts to keep larger pair of factors.
*/
void factor(int value, int *f1, int *f2) {
int j;
int i;
j = (int)sqrt(value);
for (i = j; i > 0; i--) {
if ((value % i) == 0) {
*f1 = i;
*f2 = value / i;
break;
}
}
}
/** /**
* @brief Main routine that loads a few particles and generates some output. * @brief Main routine that loads a few particles and generates some output.
...@@ -576,14 +91,12 @@ int main(int argc, char *argv[]) { ...@@ -576,14 +91,12 @@ int main(int argc, char *argv[]) {
#ifdef WITH_MPI #ifdef WITH_MPI
/* Start by initializing MPI. */ /* Start by initializing MPI. */
int res, prov; int res, prov;
if ((res = MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &prov)) != if ((res = MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &prov)) != MPI_SUCCESS)
MPI_SUCCESS)
error("Call to MPI_Init failed with error %i.", res); error("Call to MPI_Init failed with error %i.", res);
if (prov != MPI_THREAD_MULTIPLE) if (prov != MPI_THREAD_MULTIPLE)
error( error("MPI does not provide the level of threading required "
"MPI does not provide the level of threading required " "(MPI_THREAD_MULTIPLE).");
"(MPI_THREAD_MULTIPLE)."); if ((res = MPI_Comm_size(MPI_COMM_WORLD, &nr_nodes)) != MPI_SUCCESS)
if ((res = MPI_Comm_size(MPI_COMM_WORLD, &nr_nodes) != MPI_SUCCESS))
error("MPI_Comm_size failed with error %i.", res); error("MPI_Comm_size failed with error %i.", res);
if ((res = MPI_Comm_rank(MPI_COMM_WORLD, &myrank)) != MPI_SUCCESS) if ((res = MPI_Comm_rank(MPI_COMM_WORLD, &myrank)) != MPI_SUCCESS)
error("Call to MPI_Comm_rank failed with error %i.", res); error("Call to MPI_Comm_rank failed with error %i.", res);
...@@ -639,7 +152,7 @@ int main(int argc, char *argv[]) { ...@@ -639,7 +152,7 @@ int main(int argc, char *argv[]) {
break; break;
case 'o': case 'o':
with_outputs = 0; with_outputs = 0;
break; break;
case 'q': case 'q':
if (sscanf(optarg, "%d", &nr_queues) != 1) if (sscanf(optarg, "%d", &nr_queues) != 1)
error("Error parsing number of queues."); error("Error parsing number of queues.");
...@@ -678,17 +191,20 @@ int main(int argc, char *argv[]) { ...@@ -678,17 +191,20 @@ int main(int argc, char *argv[]) {
} }
#if defined(WITH_MPI) #if defined(WITH_MPI)
if (myrank == 0) message("Running with %i thread(s) per node.", nr_threads); if (myrank == 0) {
message("Running with %i thread(s) per node.", nr_threads);
message("grid set to [ %i %i %i ].", grid[0], grid[1], grid[2]);
if (nr_nodes == 1) {
message("WARNING: you are running with one MPI rank.");
message("WARNING: you should use the non-MPI version of this program." );
}
fflush(stdout);
}
#else #else
if (myrank == 0) message("Running with %i thread(s).", nr_threads); if (myrank == 0) message("Running with %i thread(s).", nr_threads);
#endif #endif
#if defined(WITH_MPI)
if (myrank == 0)
message("grid set to [ %i %i %i ].", grid[0], grid[1], grid[2]);
fflush(stdout);
#endif
/* How large are the parts? */ /* How large are the parts? */
if (myrank == 0) { if (myrank == 0) {
message("sizeof(struct part) is %li bytes.", (long int)sizeof(struct part)); message("sizeof(struct part) is %li bytes.", (long int)sizeof(struct part));
...@@ -696,7 +212,7 @@ int main(int argc, char *argv[]) { ...@@ -696,7 +212,7 @@ int main(int argc, char *argv[]) {
(long int)sizeof(struct gpart)); (long int)sizeof(struct gpart));
} }
/* Initilaize unit system */ /* Initialize unit system */
initUnitSystem(&us); initUnitSystem(&us);
if (myrank == 0) { if (myrank == 0) {
message("Unit system: U_M = %e g.", us.UnitMass_in_cgs); message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
...@@ -781,7 +297,7 @@ int main(int argc, char *argv[]) { ...@@ -781,7 +297,7 @@ int main(int argc, char *argv[]) {
// message( "cutoffs in [ %g %g ]." , s.h_min , s.h_max ); fflush(stdout); // message( "cutoffs in [ %g %g ]." , s.h_min , s.h_max ); fflush(stdout);
} }
/* Verify that each particle is in it's propper cell. */ /* Verify that each particle is in it's proper cell. */
if (myrank == 0) { if (myrank == 0) {
icount = 0; icount = 0;
space_map_cells_pre(&s, 0, &map_cellcheck, &icount); space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
...@@ -890,7 +406,7 @@ int main(int argc, char *argv[]) { ...@@ -890,7 +406,7 @@ int main(int argc, char *argv[]) {
#endif #endif
} }
/* Dump a line of agregate output. */ /* Dump a line of aggregate output. */
/* if (myrank == 0) { */ /* if (myrank == 0) { */
/* printf("%i %e %.16e %.16e %.16e %.3e %.3e %i %.3e %.3e", j, e.time, /* printf("%i %e %.16e %.16e %.16e %.3e %.3e %i %.3e %.3e", j, e.time,
*/ */
......
############################################################################### ###############################################################################
# This file is part of SWIFT. # This file is part of SWIFT.
# Coypright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Copyright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Bert Vandenbroucke (bert.vandenbroucke@ugent.be) # Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
# Matthieu Schaller (matthieu.schaller@durham.ac.uk) # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
# #
......
...@@ -34,31 +34,31 @@ do ...@@ -34,31 +34,31 @@ do
# Sod-Shock runs # Sod-Shock runs
if [ ! -e SodShock_mindt_${cpu}.dump ] if [ ! -e SodShock_mindt_${cpu}.dump ]
then then
./test_mindt -c 1.0 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1.0 > SodShock_${cpu}.dump ./swift_mindt -c 1.0 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1.0 > SodShock_${cpu}.dump
fi fi
if [ ! -e SodShock_fixed_${cpu}.dump ] if [ ! -e SodShock_fixed_${cpu}.dump ]
then then
./test_fixdt -r 1000 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1e-4 > SodShock_fixed_${cpu}.dump ./swift_fixdt -r 1000 -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -d 1e-4 > SodShock_fixed_${cpu}.dump
fi fi
# Sedov blast # Sedov blast
if [ ! -e SedovBlast_mindt_${cpu}.dump ] if [ ! -e SedovBlast_mindt_${cpu}.dump ]
then then
./test_mindt -c 0.2 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 1e-10 > SedovBlast_${cpu}.dump ./swift_mindt -c 0.2 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 1e-10 > SedovBlast_${cpu}.dump
fi fi
if [ ! -e SedovBlast_fixed_${cpu}.dump ] if [ ! -e SedovBlast_fixed_${cpu}.dump ]
then then
./test_fixdt -r 4096 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 5e-5 > SedovBlast_fixed_${cpu}.dump ./swift_fixdt -r 4096 -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -d 5e-5 > SedovBlast_fixed_${cpu}.dump
fi fi
# Cosmological volume # Cosmological volume
if [ ! -e CosmoVolume_mindt_${cpu}.dump ] if [ ! -e CosmoVolume_mindt_${cpu}.dump ]
then then
./test_mindt -c 0.01 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 0.01 > CosmoVolume_${cpu}.dump ./swift_mindt -c 0.01 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 0.01 > CosmoVolume_${cpu}.dump
fi fi
if [ ! -e CosmoVolume_fixed_${cpu}.dump ] if [ ! -e CosmoVolume_fixed_${cpu}.dump ]
then then
./test_fixdt -r 256 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 1e-8 > CosmoVolume_fixed_${cpu}.dump ./swift_fixdt -r 256 -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -d 1e-8 > CosmoVolume_fixed_${cpu}.dump
fi fi
done done
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment