Add first version of force calculation with cluster scheme

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2022-01-28 18:07:41 +01:00
parent eedcc97e4a
commit 6691803910
8 changed files with 85 additions and 85 deletions

View File

@ -45,8 +45,9 @@ void initAtom(Atom *atom)
{ {
atom->x = NULL; atom->y = NULL; atom->z = NULL; atom->x = NULL; atom->y = NULL; atom->z = NULL;
atom->vx = NULL; atom->vy = NULL; atom->vz = NULL; atom->vx = NULL; atom->vy = NULL; atom->vz = NULL;
atom->fx = NULL; atom->fy = NULL; atom->fz = NULL;
atom->cl_x = NULL; atom->cl_x = NULL;
atom->cl_v = NULL;
atom->cl_f = NULL;
atom->Natoms = 0; atom->Natoms = 0;
atom->Nlocal = 0; atom->Nlocal = 0;
atom->Nghost = 0; atom->Nghost = 0;
@ -275,9 +276,6 @@ void growAtom(Atom *atom)
atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->vy = (MD_FLOAT*) reallocate(atom->vy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); atom->vy = (MD_FLOAT*) reallocate(atom->vy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->vz = (MD_FLOAT*) reallocate(atom->vz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); atom->vz = (MD_FLOAT*) reallocate(atom->vz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->fy = (MD_FLOAT*) reallocate(atom->fy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int)); atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
} }
@ -287,4 +285,6 @@ void growClusters(Atom *atom)
atom->Nclusters_max += DELTA; atom->Nclusters_max += DELTA;
atom->clusters = (Cluster*) reallocate(atom->clusters, ALIGNMENT, atom->Nclusters_max * sizeof(Cluster), nold * sizeof(Cluster)); atom->clusters = (Cluster*) reallocate(atom->clusters, ALIGNMENT, atom->Nclusters_max * sizeof(Cluster), nold * sizeof(Cluster));
atom->cl_x = (MD_FLOAT*) reallocate(atom->cl_x, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT)); atom->cl_x = (MD_FLOAT*) reallocate(atom->cl_x, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT));
atom->cl_f = (MD_FLOAT*) reallocate(atom->cl_f, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT));
atom->cl_v = (MD_FLOAT*) reallocate(atom->cl_v, ALIGNMENT, atom->Nclusters_max * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT), nold * CLUSTER_DIM_N * 3 * sizeof(MD_FLOAT));
} }

View File

@ -33,6 +33,7 @@
#include <util.h> #include <util.h>
double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbor, Stats *stats) { double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbor, Stats *stats) {
/*
if(eam->nmax < atom->Nmax) { if(eam->nmax < atom->Nmax) {
eam->nmax = atom->Nmax; eam->nmax = atom->Nmax;
if(eam->fp != NULL) { free(eam->fp); } if(eam->fp != NULL) { free(eam->fp); }
@ -45,9 +46,11 @@ double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbo
MD_FLOAT* rhor_spline = eam->rhor_spline; MD_FLOAT* frho_spline = eam->frho_spline; MD_FLOAT* z2r_spline = eam->z2r_spline; MD_FLOAT* rhor_spline = eam->rhor_spline; MD_FLOAT* frho_spline = eam->frho_spline; MD_FLOAT* z2r_spline = eam->z2r_spline;
MD_FLOAT rdr = eam->rdr; int nr = eam->nr; int nr_tot = eam->nr_tot; MD_FLOAT rdrho = eam->rdrho; MD_FLOAT rdr = eam->rdr; int nr = eam->nr; int nr_tot = eam->nr_tot; MD_FLOAT rdrho = eam->rdrho;
int nrho = eam->nrho; int nrho_tot = eam->nrho_tot; int nrho = eam->nrho; int nrho_tot = eam->nrho_tot;
*/
double S = getTimeStamp(); double S = getTimeStamp();
LIKWID_MARKER_START("force_eam_fp"); LIKWID_MARKER_START("force_eam_fp");
/*
#pragma omp parallel for #pragma omp parallel for
for(int i = 0; i < Nlocal; i++) { for(int i = 0; i < Nlocal; i++) {
neighs = &neighbor->neighbors[i * neighbor->maxneighs]; neighs = &neighbor->neighbors[i * neighbor->maxneighs];
@ -207,6 +210,7 @@ double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbo
addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH); addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH);
} }
*/
LIKWID_MARKER_STOP("force_eam"); LIKWID_MARKER_STOP("force_eam");
double E = getTimeStamp(); double E = getTimeStamp();
return E-S; return E-S;

View File

@ -31,55 +31,45 @@
double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) { double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
int Nlocal = atom->Nlocal; int Nlocal = atom->Nlocal;
int* neighs; int* neighs;
MD_FLOAT* fx = atom->fx;
MD_FLOAT* fy = atom->fy;
MD_FLOAT* fz = atom->fz;
#ifndef EXPLICIT_TYPES
MD_FLOAT cutforcesq = param->cutforce * param->cutforce; MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
MD_FLOAT sigma6 = param->sigma6; MD_FLOAT sigma6 = param->sigma6;
MD_FLOAT epsilon = param->epsilon; MD_FLOAT epsilon = param->epsilon;
#endif
for(int i = 0; i < Nlocal; i++) { for(int ci = 0; ci < atom->Nclusters_local; ci++) {
fx[i] = 0.0; MD_FLOAT *fptr = cluster_force_ptr(ci);
fy[i] = 0.0; for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
fz[i] = 0.0; cluster_x(fptr, cii) = 0.0;
cluster_y(fptr, cii) = 0.0;
cluster_z(fptr, cii) = 0.0;
}
} }
double S = getTimeStamp(); double S = getTimeStamp();
LIKWID_MARKER_START("force"); LIKWID_MARKER_START("force");
#pragma omp parallel for #pragma omp parallel for
for(int i = 0; i < Nlocal; i++) { for(int ci = 0; ci < atom->Nclusters_local; ci++) {
neighs = &neighbor->neighbors[i * neighbor->maxneighs]; MD_FLOAT *ciptr = cluster_pos_ptr(ci);
int numneighs = neighbor->numneigh[i]; MD_FLOAT *cifptr = cluster_force_ptr(ci);
MD_FLOAT xtmp = atom_x(i); neighs = &neighbor->neighbors[ci * neighbor->maxneighs];
MD_FLOAT ytmp = atom_y(i); int numneighs = neighbor->numneigh[ci];
MD_FLOAT ztmp = atom_z(i);
for(int k = 0; k < numneighs; k++) {
int cj = neighs[k];
MD_FLOAT *cjptr = cluster_pos_ptr(cj);
for(int cii = 0; cii < CLUSTER_DIM_M; cii++) {
MD_FLOAT xtmp = cluster_x(ciptr, cii);
MD_FLOAT ytmp = cluster_y(ciptr, cii);
MD_FLOAT ztmp = cluster_z(ciptr, cii);
MD_FLOAT fix = 0; MD_FLOAT fix = 0;
MD_FLOAT fiy = 0; MD_FLOAT fiy = 0;
MD_FLOAT fiz = 0; MD_FLOAT fiz = 0;
#ifdef EXPLICIT_TYPES for(int cjj = 0; cjj < CLUSTER_DIM_N; cjj++) {
const int type_i = atom->type[i]; MD_FLOAT delx = xtmp - cluster_x(cjptr, cjj);
#endif MD_FLOAT dely = ytmp - cluster_y(cjptr, cjj);
MD_FLOAT delz = ztmp - cluster_z(cjptr, cjj);
for(int k = 0; k < numneighs; k++) {
int j = neighs[k];
MD_FLOAT delx = xtmp - atom_x(j);
MD_FLOAT dely = ytmp - atom_y(j);
MD_FLOAT delz = ztmp - atom_z(j);
MD_FLOAT rsq = delx * delx + dely * dely + delz * delz; MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
#ifdef EXPLICIT_TYPES
const int type_j = atom->type[j];
const int type_ij = type_i * atom->ntypes + type_j;
const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij];
const MD_FLOAT sigma6 = atom->sigma6[type_ij];
const MD_FLOAT epsilon = atom->epsilon[type_ij];
#endif
if(rsq < cutforcesq) { if(rsq < cutforcesq) {
MD_FLOAT sr2 = 1.0 / rsq; MD_FLOAT sr2 = 1.0 / rsq;
MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6; MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6;
@ -90,9 +80,11 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
} }
} }
fx[i] += fix; cluster_x(cifptr, cii) += fix;
fy[i] += fiy; cluster_y(cifptr, cii) += fiy;
fz[i] += fiz; cluster_z(cifptr, cii) += fiz;
}
}
addStat(stats->total_force_neighs, numneighs); addStat(stats->total_force_neighs, numneighs);
addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH); addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH);

View File

@ -42,8 +42,9 @@ typedef struct {
int Nclusters, Nclusters_local, Nclusters_ghost, Nclusters_max; int Nclusters, Nclusters_local, Nclusters_ghost, Nclusters_max;
MD_FLOAT *x, *y, *z; MD_FLOAT *x, *y, *z;
MD_FLOAT *vx, *vy, *vz; MD_FLOAT *vx, *vy, *vz;
MD_FLOAT *fx, *fy, *fz;
MD_FLOAT *cl_x; MD_FLOAT *cl_x;
MD_FLOAT *cl_v;
MD_FLOAT *cl_f;
int *border_map; int *border_map;
int *type; int *type;
int ntypes; int ntypes;
@ -61,7 +62,9 @@ extern int readAtom(Atom*, Parameter*);
extern void growAtom(Atom*); extern void growAtom(Atom*);
extern void growClusters(Atom*); extern void growClusters(Atom*);
#define cluster_ptr(ci) &(atom->cl_x[(ci) * CLUSTER_DIM_N * 3]) #define cluster_pos_ptr(ci) &(atom->cl_x[(ci) * CLUSTER_DIM_N * 3])
#define cluster_velocity_ptr(ci) &(atom->cl_v[(ci) * CLUSTER_DIM_N * 3])
#define cluster_force_ptr(ci) &(atom->cl_f[(ci) * CLUSTER_DIM_N * 3])
#ifdef AOS #ifdef AOS
#define POS_DATA_LAYOUT "AoS" #define POS_DATA_LAYOUT "AoS"

View File

@ -129,30 +129,33 @@ double reneighbour(
return E-S; return E-S;
} }
void initialIntegrate(Parameter *param, Atom *atom) void initialIntegrate(Parameter *param, Atom *atom) {
{ for(int ci = 0; ci < atom->Nclusters_local; ci++) {
MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; MD_FLOAT *ciptr = cluster_pos_ptr(ci);
MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz; MD_FLOAT *civptr = cluster_velocity_ptr(ci);
MD_FLOAT *cifptr = cluster_force_ptr(ci);
for(int i = 0; i < atom->Nlocal; i++) { for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
vx[i] += param->dtforce * fx[i]; cluster_x(civptr, cii) += param->dtforce * cluster_x(cifptr, cii);
vy[i] += param->dtforce * fy[i]; cluster_y(civptr, cii) += param->dtforce * cluster_y(cifptr, cii);
vz[i] += param->dtforce * fz[i]; cluster_z(civptr, cii) += param->dtforce * cluster_z(cifptr, cii);
atom_x(i) = atom_x(i) + param->dt * vx[i]; cluster_x(ciptr, cii) += param->dt * cluster_x(civptr, cii);
atom_y(i) = atom_y(i) + param->dt * vy[i]; cluster_y(ciptr, cii) += param->dt * cluster_y(civptr, cii);
atom_z(i) = atom_z(i) + param->dt * vz[i]; cluster_z(ciptr, cii) += param->dt * cluster_z(civptr, cii);
}
} }
} }
void finalIntegrate(Parameter *param, Atom *atom) void finalIntegrate(Parameter *param, Atom *atom) {
{ for(int ci = 0; ci < atom->Nclusters_local; ci++) {
MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; MD_FLOAT *civptr = cluster_velocity_ptr(ci);
MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz; MD_FLOAT *cifptr = cluster_force_ptr(ci);
for(int i = 0; i < atom->Nlocal; i++) { for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
vx[i] += param->dtforce * fx[i]; cluster_x(civptr, cii) += param->dtforce * cluster_x(cifptr, cii);
vy[i] += param->dtforce * fy[i]; cluster_y(civptr, cii) += param->dtforce * cluster_y(cifptr, cii);
vz[i] += param->dtforce * fz[i]; cluster_z(civptr, cii) += param->dtforce * cluster_z(cifptr, cii);
}
} }
} }
@ -257,13 +260,11 @@ int main(int argc, char** argv)
#if defined(MEM_TRACER) || defined(INDEX_TRACER) #if defined(MEM_TRACER) || defined(INDEX_TRACER)
traceAddresses(&param, &atom, &neighbor, n + 1); traceAddresses(&param, &atom, &neighbor, n + 1);
#endif #endif
/*
if(param.force_field == FF_EAM) { if(param.force_field == FF_EAM) {
timer[FORCE] = computeForceEam(&eam, &param, &atom, &neighbor, &stats); timer[FORCE] = computeForceEam(&eam, &param, &atom, &neighbor, &stats);
} else { } else {
timer[FORCE] = computeForceLJ(&param, &atom, &neighbor, &stats); timer[FORCE] = computeForceLJ(&param, &atom, &neighbor, &stats);
} }
*/
timer[NEIGH] = 0.0; timer[NEIGH] = 0.0;
timer[TOTAL] = getTimeStamp(); timer[TOTAL] = getTimeStamp();
@ -273,7 +274,7 @@ int main(int argc, char** argv)
} }
for(int n = 0; n < param.ntimes; n++) { for(int n = 0; n < param.ntimes; n++) {
//initialIntegrate(&param, &atom); initialIntegrate(&param, &atom);
if((n + 1) % param.every) { if((n + 1) % param.every) {
updatePbc(&atom, &param); updatePbc(&atom, &param);
@ -285,15 +286,13 @@ int main(int argc, char** argv)
traceAddresses(&param, &atom, &neighbor, n + 1); traceAddresses(&param, &atom, &neighbor, n + 1);
#endif #endif
/*
if(param.force_field == FF_EAM) { if(param.force_field == FF_EAM) {
timer[FORCE] += computeForceEam(&eam, &param, &atom, &neighbor, &stats); timer[FORCE] += computeForceEam(&eam, &param, &atom, &neighbor, &stats);
} else { } else {
timer[FORCE] += computeForceLJ(&param, &atom, &neighbor, &stats); timer[FORCE] += computeForceLJ(&param, &atom, &neighbor, &stats);
} }
*/
//finalIntegrate(&param, &atom); finalIntegrate(&param, &atom);
if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) {
computeThermo(n + 1, &param, &atom); computeThermo(n + 1, &param, &atom);

View File

@ -184,8 +184,8 @@ MD_FLOAT getBoundingBoxDistanceSq(Atom *atom, int ci, int cj) {
} }
int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) { int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) {
MD_FLOAT *ciptr = cluster_ptr(ci); MD_FLOAT *ciptr = cluster_pos_ptr(ci);
MD_FLOAT *cjptr = cluster_ptr(cj); MD_FLOAT *cjptr = cluster_pos_ptr(cj);
for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) { for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
for(int cjj = 0; cjj < atom->clusters[cj].natoms; cjj++) { for(int cjj = 0; cjj < atom->clusters[cj].natoms; cjj++) {
@ -422,7 +422,7 @@ void buildClusters(Atom *atom) {
growClusters(atom); growClusters(atom);
} }
MD_FLOAT *cptr = cluster_ptr(ci); MD_FLOAT *cptr = cluster_pos_ptr(ci);
MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY; MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY;
MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY; MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY;
MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY; MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY;
@ -493,7 +493,7 @@ void binClusters(Atom *atom) {
} }
for(int ci = 0; ci < atom->Nclusters_ghost && !resize; ci++) { for(int ci = 0; ci < atom->Nclusters_ghost && !resize; ci++) {
MD_FLOAT *cptr = cluster_ptr(nlocal + ci); MD_FLOAT *cptr = cluster_pos_ptr(nlocal + ci);
MD_FLOAT xtmp, ytmp; MD_FLOAT xtmp, ytmp;
int ix = -1, iy = -1; int ix = -1, iy = -1;
@ -539,7 +539,7 @@ void updateSingleAtoms(Atom *atom) {
int Natom = 0; int Natom = 0;
for(int ci = 0; ci < atom->Nclusters_local; ci++) { for(int ci = 0; ci < atom->Nclusters_local; ci++) {
MD_FLOAT *cptr = cluster_ptr(ci); MD_FLOAT *cptr = cluster_pos_ptr(ci);
for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) { for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
atom_x(Natom) = cluster_x(cptr, cii); atom_x(Natom) = cluster_x(cptr, cii);

View File

@ -51,8 +51,8 @@ void updatePbc(Atom *atom, Parameter *param) {
MD_FLOAT zprd = param->zprd; MD_FLOAT zprd = param->zprd;
for(int ci = 0; ci < atom->Nclusters_ghost; ci++) { for(int ci = 0; ci < atom->Nclusters_ghost; ci++) {
MD_FLOAT *cptr = cluster_ptr(nlocal + ci); MD_FLOAT *cptr = cluster_pos_ptr(nlocal + ci);
MD_FLOAT *bmap_cptr = cluster_ptr(border_map[ci]); MD_FLOAT *bmap_cptr = cluster_pos_ptr(border_map[ci]);
for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) { for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
cluster_x(cptr, cii) = cluster_x(bmap_cptr, cii) + atom->PBCx[ci] * xprd; cluster_x(cptr, cii) = cluster_x(bmap_cptr, cii) + atom->PBCx[ci] * xprd;

View File

@ -30,7 +30,7 @@ void traceAddresses(Parameter *param, Atom *atom, Neighbor *neighbor, int timest
INDEX_TRACER_INIT; INDEX_TRACER_INIT;
int Nlocal = atom->Nlocal; int Nlocal = atom->Nlocal;
int* neighs; int* neighs;
MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; //MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
INDEX_TRACE_NATOMS(Nlocal, atom->Nghost, neighbor->maxneighs); INDEX_TRACE_NATOMS(Nlocal, atom->Nghost, neighbor->maxneighs);
for(int i = 0; i < Nlocal; i++) { for(int i = 0; i < Nlocal; i++) {
@ -60,12 +60,14 @@ void traceAddresses(Parameter *param, Atom *atom, Neighbor *neighbor, int timest
#endif #endif
} }
/*
MEM_TRACE(fx[i], 'R'); MEM_TRACE(fx[i], 'R');
MEM_TRACE(fx[i], 'W'); MEM_TRACE(fx[i], 'W');
MEM_TRACE(fy[i], 'R'); MEM_TRACE(fy[i], 'R');
MEM_TRACE(fy[i], 'W'); MEM_TRACE(fy[i], 'W');
MEM_TRACE(fz[i], 'R'); MEM_TRACE(fz[i], 'R');
MEM_TRACE(fz[i], 'W'); MEM_TRACE(fz[i], 'W');
*/
} }
INDEX_TRACER_END; INDEX_TRACER_END;