From e38490cf71025d05023aaa43309535e2470fc1ae Mon Sep 17 00:00:00 2001 From: Rob van der Wijngaart Date: Tue, 29 Dec 2015 15:13:16 -0800 Subject: [PATCH 01/14] Prepping for LCG in PIC. --- common/make.common | 3 +++ 1 file changed, 3 insertions(+) diff --git a/common/make.common b/common/make.common index c49837f24..3cf56936e 100644 --- a/common/make.common +++ b/common/make.common @@ -23,6 +23,9 @@ endif wtime.o:$(COMMON)/wtime.c $(CCOMPILER) $(CFLAGS) $(TUNEFLAGS) $(INCLUDEPATHSPLUS) -c $< +lcg.o:$(COMMON)/lcg.c + $(CCOMPILER) $(CFLAGS) $(TUNEFLAGS) $(INCLUDEPATHSPLUS) -c $< + MPI_bail_out.o:$(COMMON)/MPI_bail_out.c $(CCOMPILER) $(CFLAGS) $(TUNEFLAGS) $(INCLUDEPATHSPLUS) -c $< From 2a7d5c88d17aeaa9356b15a5d23cbc4a3c9f3cde Mon Sep 17 00:00:00 2001 From: Rob van der Wijngaart Date: Tue, 29 Dec 2015 16:50:37 -0800 Subject: [PATCH 02/14] Introducing serial PIC PRK, which requires the lcg (Linear Congruential Generator) stuff. --- Makefile | 2 + SERIAL/PIC/Makefile | 33 ++ SERIAL/PIC/pic.c | 924 +++++++++++++++++++++++++++++++++++++++ SHMEM/Stencil/Makefile | 4 +- SHMEM/Synch_p2p/Makefile | 4 +- SHMEM/Transpose/Makefile | 4 +- common/lcg.c | 152 +++++++ include/lcg.h | 45 ++ scripts/small/runserial | 1 + 9 files changed, 1163 insertions(+), 6 deletions(-) create mode 100644 SERIAL/PIC/Makefile create mode 100644 SERIAL/PIC/pic.c create mode 100644 common/lcg.c create mode 100644 include/lcg.h diff --git a/Makefile b/Makefile index 1151e29ea..fdc781f22 100644 --- a/Makefile +++ b/Makefile @@ -185,6 +185,7 @@ allserial: cd SERIAL/Branch; $(MAKE) branch "DEFAULT_OPT_FLAGS = $(default_opt_flags)" \ "MATRIX_RANK = $(matrix_rank)" \ "NUMBER_OF_FUNCTIONS = $(number_of_functions)" + cd SERIAL/PIC; $(MAKE) pic "DEFAULT_OPT_FLAGS = $(default_opt_flags)" clean: cd MPI1/DGEMM; $(MAKE) clean @@ -260,6 +261,7 @@ clean: cd SERIAL/Sparse; $(MAKE) clean cd SERIAL/Synch_p2p; $(MAKE) clean cd SERIAL/Branch; $(MAKE) clean + cd SERIAL/PIC; $(MAKE) clean rm -f stats.json veryclean: clean diff --git a/SERIAL/PIC/Makefile b/SERIAL/PIC/Makefile new file mode 100644 index 000000000..a3e9e9ec6 --- /dev/null +++ b/SERIAL/PIC/Makefile @@ -0,0 +1,33 @@ +include ../../common/SERIAL.defs +COMOBJS += lcg.o + +##### User configurable options ##### +#uncomment any of the following flags (or change values) to change defaults + +OPTFLAGS = $(DEFAULT_OPT_FLAGS) +#OPTFLAGS=-axCORE-AVX2 -O3 -restrict +#description: change above into something that is a decent optimization on you system + +#RESTRICTFLAG = -DRESTRICT_KEYWORD +#description: the "restrict" keyword can be used on IA platforms to disambiguate +# data accessed through pointers + +DEBUGFLAG = -DVERBOSE +#description: default diagnostic style is silent + +USERFLAGS = +#description: parameter to specify optional flags + +#set the following variables for custom libraries and/or other objects +EXTOBJS = +LIBS = +LIBPATHS = +INCLUDEPATHS = + +### End User configurable options ### +TUNEFLAGS = $(RESTRICTFLAG) $(DEBUGFLAG) $(USERFLAGS) +PROGRAM = pic +OBJS = $(PROGRAM).o $(COMOBJS) + +include ../../common/make.common + diff --git a/SERIAL/PIC/pic.c b/SERIAL/PIC/pic.c new file mode 100644 index 000000000..41ccf4d61 --- /dev/null +++ b/SERIAL/PIC/pic.c @@ -0,0 +1,924 @@ +/* +Copyright (c) 2015, Intel Corporation + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +* Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +* Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. +* Neither the name of Intel Corporation nor the names of its + contributors may be used to endorse or promote products + derived from this software without specific prior written + permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. +*/ + +/******************************************************************* + +NAME: PIC + +PURPOSE: This program tests the efficiency with which a cloud of + charged particles can be moved through a spatially fixed + collection of charges located at the vertices of a square + equi-spaced grid. It is a proxy for a component of a + particle-in-cell methjod + +USAGE: <#simulation steps> <#particles> \ + \ + \ + [ ] + + + The output consists of diagnostics to make sure the + algorithm worked, and of timing statistics. + +FUNCTIONS CALLED: + + Other than standard C functions, the following functions are used in + this program: + wtime() + bad_path() + LCG_next() + +HISTORY: - Written by Evangelos Georganas, August 2015. + - RvdW: Refactored to make the code PRK conforming, December 2015 + +**********************************************************************************/ + +#include +#include + +#define VERIF_AUX + +#include +#include +#include + +#define QG(i,j) Qgrid[(j)*(g)+i] +#define mass_inv 1.0 +#define Q 1.0 +#define epsilon 0.000001 +#define dt 1.0 + +#define SUCCESS 1 +#define FAILURE 0 + +#define GEOMETRIC 0 +#define SINUSOIDAL 1 +#define LINEAR 2 +#define PATCH 3 +#define UNDEFINED 4 + +typedef struct { + int64_t xleft; + int64_t xright; + int64_t ybottom; + int64_t ytop; +} bbox_t; + +/* Particle data structure */ +typedef struct particle_t { + double x; + double y; + double v_x; + double v_y; + double q; + /* The following variables are used only for verification/debug purposes */ +#ifdef VERIF_AUX + double x0; + double y0; + int64_t k; // determines how many cells particles move per time step in the x direction + int64_t m; // determines how many cells particles move per time step in the y direction + int64_t initTimestamp; +#endif +} particle_t; + +/* Initializes the grid of charges + We follow a column major format for the grid. Note that this may affect cache performance, depending on access pattern of particles. */ + +/* The grid is indexed in this way: + + y ^ + | + | + | + | + | + | + (0,0)--------------> x */ + +double *initializeGrid(int64_t g) +{ + double *Qgrid; + int64_t y, x; + + Qgrid = (double*) malloc(g*g*sizeof(double)); + if (Qgrid == NULL) { + printf("ERROR: Could not allocate space for grid\n"); + exit(EXIT_FAILURE); + } + + /* initialization with dipoles */ + for (x=0; x 0) ? p.x0 + total_steps * (2*p.k+1) : p.x0 - total_steps * (2*p.k+1) ; + y_T = p.y0 + p.m * total_steps; + + x_periodic = fmod(x_T+total_steps *(2*p.k+1)*L, L); + y_periodic = fmod(y_T+total_steps *fabs(p.m)*L, L); + + if ( fabs(p.x - x_periodic) > epsilon || fabs(p.y - y_periodic) > epsilon) { + return FAILURE; + } + + return SUCCESS; +} + +/* Removes particles from a specified area of the simulation domain */ +particle_t *remove_particles(int64_t removal_timestep, bbox_t patch, int64_t *n, particle_t *particles, int *partial_correctness, double *Qgrid, int64_t g) +{ + int64_t pos = 0, i; + /* The boundaries of the simulation domain where we have to remove the particles */ + double left_boundary = patch.xleft; + double right_boundary = patch.xright; + double top_boundary = patch.ytop; + double bottom_boundary = patch.ybottom; + particle_t *new_particles_array; + + for (i = 0; i < (*n); i++) { + if ( (particles[i].x > left_boundary) && (particles[i].x < right_boundary) && + (particles[i].y > bottom_boundary) && (particles[i].y < top_boundary)) { + /* We should remove and verify this particle */ + (*partial_correctness) *= verifyParticle(particles[i], removal_timestep, Qgrid, g); + } else { + /* We should keep the particle */ + particles[pos] = particles[i]; + pos++; + } + } + + (*n) = pos; + + + return particles; +} + + +/* Computes the Coulomb force among two charges q1 and q2 */ +int computeCoulomb(double x_dist, double y_dist, double q1, double q2, double *fx, double *fy) +{ + double r2 = x_dist * x_dist + y_dist * y_dist; + double r = sqrt(r2); + double cos_theta = x_dist / r; + double sin_theta = y_dist / r; + double f_coulomb = q1 * q2 / r2; + + (*fx) = f_coulomb * cos_theta; + (*fy) = f_coulomb * sin_theta; + + return 0; +} + +/* Computes the total Coulomb force on a particle exerted from the charges of the corresponding cell */ +int computeTotalForce(particle_t p, int64_t g, double *Qgrid, double *fx, double *fy) +{ + int64_t y, x, k; + double tmp_fx, tmp_fy, rel_y, rel_x; + double tmp_res_x = 0.0; + double tmp_res_y = 0.0; + + /* Coordinates of the cell containing the particle */ + y = (int64_t) floor(p.y); + x = (int64_t) floor(p.x); + rel_x = p.x - x; + rel_y = p.y - y; + + /* Coulomb force from top-left charge */ + computeCoulomb(rel_x, rel_y, p.q, QG(y,x), &tmp_fx, &tmp_fy); + tmp_res_x += tmp_fx; + tmp_res_y += tmp_fy; + + /* Coulomb force from bottom-left charge */ + computeCoulomb(rel_x, 1.0-rel_y, p.q, QG(y+1,x), &tmp_fx, &tmp_fy); + tmp_res_x += tmp_fx; + tmp_res_y -= tmp_fy; + + /* Coulomb force from top-right charge */ + computeCoulomb(1.0-rel_x, rel_y, p.q, QG(y,x+1), &tmp_fx, &tmp_fy); + tmp_res_x -= tmp_fx; + tmp_res_y += tmp_fy; + + /* Coulomb force from bottom-right charge */ + computeCoulomb(1.0-rel_x, 1.0-rel_y, p.q, QG(y+1,x+1), &tmp_fx, &tmp_fy); + tmp_res_x -= tmp_fx; + tmp_res_y -= tmp_fy; + + (*fx) = tmp_res_x; + (*fy) = tmp_res_y; + + return 0; +} + +/* Moves a particle given the total acceleration */ +int moveParticle(particle_t *particle, double ax, double ay, double L) +{ + + /* Update particle positions, taking into account periodic boundaries */ + particle->x = fmod(particle->x + particle->v_x*dt + 0.5*ax*dt*dt + L, L); + particle->y = fmod(particle->y + particle->v_y*dt + 0.5*ay*dt*dt + L, L); + + /* Update velocities */ + particle->v_x += ax * dt; + particle->v_y += ay * dt; + + return 0; +} + + +int bad_patch(bbox_t *patch, bbox_t *patch_contain) { + if (patch->xleft>=patch->xright || patch->ybottom>=patch->ytop) return(1); + if (patch_contain) { + if (patch->xleft xleft || patch->xright>patch_contain->xright) return(2); + if (patch->ybottomybottom || patch->ytop >patch_contain->ytop) return(3); + } + return(0); +} + +int main(int argc, char ** argv) { + + int args_used = 1; // keeps track of # consumed arguments + int64_t g; // dimension of grid in points + int64_t T; // total number of simulation steps + int64_t n; // total number of particles in the simulation + int64_t n_old; // number of particled before removal + char *init_mode; // Initialization mode for particles + double rho; // rho parameter for the initial geometric particle distribution + int64_t k, m; // determine initial horizontal and vertical velocityof particles -- + // (2*k)+1 cells per time step + int64_t particle_mode, alpha, beta; + bbox_t grid_patch, init_patch, injection_patch, removal_patch; + int removal_mode = 0, injection_mode = 0, injection_timestep, removal_timestep, particles_per_cell; + int partial_correctness = 1; + int64_t L; // dimension of grid in cells + double *Qgrid;// the grid is represented as an array of charges + particle_t *particles; // the particles array + int64_t t, i, j; + double fx, fy, ax, ay, simulation_time; + int correct_simulation = 1; + int error; + int64_t particle_steps, particles_added; + double avg_time; + + printf("Parallel Research Kernels Version %s\n", PRKVERSION); + printf("Serial Particle-in-Cell execution on 2D grid\n"); + + /******************************************************************************* + ** process and test input parameters + ********************************************************************************/ + + if (argc<6) { + printf("Usage: %s <#simulation steps> <#particles> ", argv[0]); + printf("\n"); + printf(" [ ]\n"); + printf(" init mode \"GEOMETRIC\" parameters: \n"); + printf(" \"SINUSOIDAL\" parameters: none\n"); + printf(" \"LINEAR\" parameters: \n"); + printf(" \"PATCH\" parameters: \n"); + printf(" population change mode \"INJECTION\" parameters: <# particles>