From 71ea8dbb0e2d3c6c8977eeca914f31006991b83d Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Wed, 19 Aug 2020 09:00:35 +0200 Subject: [PATCH] Refactor. Fix bug in atom initialization. --- include_CLANG.mk | 2 +- src/allocate.c | 3 +- src/atom.c | 22 ++++++++--- src/includes/atom.h | 1 + src/includes/neighbor.h | 2 +- src/main.c | 86 ++++++++++++++++++++++++++++++----------- src/neighbor.c | 2 +- src/pbc.c | 32 +++++++++------ src/thermo.c | 1 + util/mdBench.c | 2 +- 10 files changed, 108 insertions(+), 45 deletions(-) diff --git a/include_CLANG.mk b/include_CLANG.mk index d2992b7..0ff0922 100644 --- a/include_CLANG.mk +++ b/include_CLANG.mk @@ -7,7 +7,7 @@ ANSI_CFLAGS += -std=c99 ANSI_CFLAGS += -pedantic ANSI_CFLAGS += -Wextra -CFLAGS = -O3 #-g $(ANSI_CFLAGS) +CFLAGS = -O3 $(ANSI_CFLAGS) #-g ASFLAGS = -masm=intel CXXFLAGS = $(CFLAGS) FCFLAGS = diff --git a/src/allocate.c b/src/allocate.c index dca0c97..a3727c9 100644 --- a/src/allocate.c +++ b/src/allocate.c @@ -64,8 +64,7 @@ void* reallocate ( size_t newBytesize, size_t oldBytesize) { - void* newarray; - newarray = allocate(alignment, newBytesize); + void* newarray = allocate(alignment, newBytesize); if(ptr != NULL) { memcpy(newarray, ptr, oldBytesize); diff --git a/src/atom.c b/src/atom.c index b61060c..7476a17 100644 --- a/src/atom.c +++ b/src/atom.c @@ -24,15 +24,27 @@ * * ======================================================================================= */ +#include +#include #include -#include -#include #include +#include #include #define DELTA 20000 +void initAtom(Atom *atom) +{ + atom->x = NULL; atom->y = NULL; atom->z = NULL; + atom->vx = NULL; atom->vy = NULL; atom->vz = NULL; + atom->fx = NULL; atom->fy = NULL; atom->fz = NULL; + atom->Natoms = 0; + atom->Nlocal = 0; + atom->Nghost = 0; + atom->Nmax = 0; +} + void createAtom(Atom *atom, Parameter *param) { double xlo = 0.0; double xhi = param->xprd; @@ -128,9 +140,9 @@ void growAtom(Atom *atom) int nold = atom->Nmax; atom->Nmax += DELTA; - atom->x = (double*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); - atom->y = (double*) reallocate(atom->y, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); - atom->z = (double*) reallocate(atom->z, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->x = (double*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->y = (double*) reallocate(atom->y, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->z = (double*) reallocate(atom->z, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); atom->vx = (double*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); atom->vy = (double*) reallocate(atom->vy, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); atom->vz = (double*) reallocate(atom->vz, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); diff --git a/src/includes/atom.h b/src/includes/atom.h index 168e556..63f3512 100644 --- a/src/includes/atom.h +++ b/src/includes/atom.h @@ -36,6 +36,7 @@ typedef struct { double *fx, *fy, *fz; } Atom; +extern void initAtom(Atom*); extern void createAtom(Atom*, Parameter*); extern void growAtom(Atom*); #endif diff --git a/src/includes/neighbor.h b/src/includes/neighbor.h index 8dc1016..dfad5d1 100644 --- a/src/includes/neighbor.h +++ b/src/includes/neighbor.h @@ -39,7 +39,7 @@ typedef struct { } Neighbor; extern void initNeighbor(Neighbor*, Parameter*); -extern void setupNeighbor(Parameter*); +extern void setupNeighbor(); extern void binatoms(Atom*); extern void buildNeighbor(Atom*, Neighbor*); #endif diff --git a/src/main.c b/src/main.c index 841bf4f..73b28de 100644 --- a/src/main.c +++ b/src/main.c @@ -33,6 +33,7 @@ #include #include +#include #include #include #include @@ -42,16 +43,16 @@ #define HLINE "----------------------------------------------------------------------------\n" -#define FACTOR 0.999 -#define SMALL 1.0e-6 -#define DELTA 20000 +typedef enum { + TOTAL = 0, + NEIGH, + FORCE, + NUMTIMER +} timertype; -void init(Atom *atom, Parameter *param) + +void init(Parameter *param) { - atom->x = NULL; atom->y = NULL; atom->z = NULL; - atom->vx = NULL; atom->vy = NULL; atom->vz = NULL; - atom->fx = NULL; atom->fy = NULL; atom->fz = NULL; - param->epsilon = 1.0; param->sigma6 = 1.0; param->rho = 0.8442; @@ -73,6 +74,19 @@ void init(Atom *atom, Parameter *param) param->zprd = param->nz * lattice; } +void setup(Parameter *param, Atom *atom, Neighbor *neighbor){ + initAtom(atom); + initNeighbor(neighbor, param); + initPbc(); + setupNeighbor(); + createAtom(atom, param); + setupThermo(param, atom->Natoms); + adjustThermo(param, atom); + setupPbc(atom, param); + updatePbc(atom, param); + buildNeighbor(atom, neighbor); +} + void initialIntegrate(Parameter *param, Atom *atom) { double* x = atom->x; double* y = atom->y; double* z = atom->z; @@ -101,7 +115,7 @@ void finalIntegrate(Parameter *param, Atom *atom) } } -void computeForce(Neighbor *neighbor, Atom *atom, Parameter *param) +void computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) { int Nlocal = atom->Natoms; int* neighs; @@ -164,26 +178,52 @@ void printAtomState(Atom *atom) } } - int main (int argc, char** argv) { + double timer[NUMTIMER]; + double S; Atom atom; Neighbor neighbor; Parameter param; - init(&atom, ¶m); - initNeighbor(&neighbor, ¶m); - initPbc(); - setupNeighbor(¶m); - createAtom(&atom, ¶m); - setupThermo(¶m, atom.Natoms); - adjustThermo(¶m, &atom); - /* printAtomState(&atom); */ - setupPbc(&atom, ¶m); - updatePbc(&atom, ¶m); - buildNeighbor(&atom, &neighbor); + init(¶m); + + for(int i = 0; i < argc; i++) + { + if((strcmp(argv[i], "-n") == 0) || (strcmp(argv[i], "--nsteps") == 0)) + { + param.ntimes = atoi(argv[++i]); + continue; + } + if((strcmp(argv[i], "-nx") == 0)) + { + param.nx = atoi(argv[++i]); + continue; + } + if((strcmp(argv[i], "-ny") == 0)) + { + param.ny = atoi(argv[++i]); + continue; + } + if((strcmp(argv[i], "-nz") == 0)) + { + param.nz = atoi(argv[++i]); + continue; + } + if((strcmp(argv[i], "-h") == 0) || (strcmp(argv[i], "--help") == 0)) + { + printf("MD Bench: A minimalistic re-implementation of miniMD\n"); + printf(HLINE); + printf("-n / --nsteps : set number of timesteps for simulation\n"); + printf("-nx/-ny/-nz : set linear dimension of systembox in x/y/z direction\n"); + printf(HLINE); + exit(EXIT_SUCCESS); + } + } + + setup(¶m, &atom, &neighbor); computeThermo(0, ¶m, &atom); - computeForce(&neighbor, &atom, ¶m); + computeForce(¶m, &atom, &neighbor); for(int n = 0; n < param.ntimes; n++) { @@ -199,7 +239,7 @@ int main (int argc, char** argv) buildNeighbor(&atom, &neighbor); } - computeForce(&neighbor, &atom, ¶m); + computeForce(¶m, &atom, &neighbor); finalIntegrate(¶m, &atom); if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { diff --git a/src/neighbor.c b/src/neighbor.c index a39cb48..aae4f11 100644 --- a/src/neighbor.c +++ b/src/neighbor.c @@ -76,7 +76,7 @@ void initNeighbor(Neighbor *neighbor, Parameter *param) neighbor->neighbors = NULL; } -void setupNeighbor(Parameter *param) +void setupNeighbor() { double coord; int mbinxhi, mbinyhi, mbinzhi; diff --git a/src/pbc.c b/src/pbc.c index ba4a0ae..a2f8011 100644 --- a/src/pbc.c +++ b/src/pbc.c @@ -39,17 +39,7 @@ static int *PBCx, *PBCy, *PBCz; static void growPbc(); -void growPbc() -{ - int nold = NmaxGhost; - NmaxGhost += DELTA; - - BorderMap = (int*) reallocate(BorderMap, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); - PBCx = (int*) reallocate(PBCx, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); - PBCy = (int*) reallocate(PBCy, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); - PBCz = (int*) reallocate(PBCz, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); -} - +/* exported subroutines */ void initPbc() { NmaxGhost = 0; @@ -57,6 +47,8 @@ void initPbc() PBCx = NULL; PBCy = NULL; PBCz = NULL; } +/* update coordinates of ghost atoms */ +/* uses mapping created in setupPbc */ void updatePbc(Atom *atom, Parameter *param) { int nlocal = atom->Nlocal; @@ -74,6 +66,8 @@ void updatePbc(Atom *atom, Parameter *param) } } +/* relocate atoms that have left domain according + * to periodic boundary conditions */ void updateAtomsPbc(Atom *atom, Parameter *param) { double* x = atom->x; @@ -105,6 +99,10 @@ void updateAtomsPbc(Atom *atom, Parameter *param) } } +/* setup periodic boundary conditions by + * defining ghost atoms around domain + * only creates mapping and coordinate corrections + * that are then enforced in updatePbc */ #define ADDGHOST(dx,dy,dz) Nghost++; BorderMap[Nghost] = i; PBCx[Nghost] = dx; PBCy[Nghost] = dy; PBCz[Nghost] = dz; void setupPbc(Atom *atom, Parameter *param) { @@ -159,3 +157,15 @@ void setupPbc(Atom *atom, Parameter *param) // increase by one to make it the ghost atom count atom->Nghost = Nghost + 1; } + +/* internal subroutines */ +void growPbc() +{ + int nold = NmaxGhost; + NmaxGhost += DELTA; + + BorderMap = (int*) reallocate(BorderMap, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); + PBCx = (int*) reallocate(PBCx, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); + PBCy = (int*) reallocate(PBCy, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); + PBCz = (int*) reallocate(PBCz, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); +} diff --git a/src/thermo.c b/src/thermo.c index b3bd651..b35884d 100644 --- a/src/thermo.c +++ b/src/thermo.c @@ -44,6 +44,7 @@ static double p_act; static double e_act; static int mstat; +/* exported subroutines */ void setupThermo(Parameter *param, int natoms) { int maxstat = param->ntimes / param->nstat + 2; diff --git a/util/mdBench.c b/util/mdBench.c index 8138ef6..91d040a 100644 --- a/util/mdBench.c +++ b/util/mdBench.c @@ -325,7 +325,7 @@ void init(Neighbor *neighbor, Parameter *param) param->dt = 0.005; param->nx = 32; param->ny = 32; - param->nz = 32; + param->nz = 64; param->cutforce = 2.5; param->temp = 1.44; param->nstat = 100;