From bb599c9ea8d3aad45323d9caaa6ec035bd704c6e Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Tue, 5 Jul 2022 15:33:31 +0200 Subject: [PATCH] Add first version of DEM kernel Signed-off-by: Rafael Ravedutti --- lammps/atom.c | 13 ++- lammps/force_dem.c | 232 +++++++++++++++++++++++++++++++++++++++++ lammps/includes/atom.h | 4 + 3 files changed, 245 insertions(+), 4 deletions(-) create mode 100644 lammps/force_dem.c diff --git a/lammps/atom.c b/lammps/atom.c index d6fac51..d923b0b 100644 --- a/lammps/atom.c +++ b/lammps/atom.c @@ -41,8 +41,7 @@ #define MAX(a,b) ((a) > (b) ? (a) : (b)) #endif -void initAtom(Atom *atom) -{ +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; @@ -56,10 +55,12 @@ void initAtom(Atom *atom) atom->sigma6 = NULL; atom->cutforcesq = NULL; atom->cutneighsq = NULL; + atom->radius = 1.0; + atom->av = NULL; + atom->r = NULL; } -void createAtom(Atom *atom, Parameter *param) -{ +void createAtom(Atom *atom, Parameter *param) { MD_FLOAT xlo = 0.0; MD_FLOAT xhi = param->xprd; MD_FLOAT ylo = 0.0; MD_FLOAT yhi = param->yprd; MD_FLOAT zlo = 0.0; MD_FLOAT zhi = param->zprd; @@ -453,4 +454,8 @@ void growAtom(Atom *atom) { 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->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int)); + #ifdef DEM + atom->av = (MD_FLOAT*) reallocate(atom->av, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3); + atom->r = (MD_FLOAT*) reallocate(atom->r, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 4, nold * sizeof(MD_FLOAT) * 4); + #endif } diff --git a/lammps/force_dem.c b/lammps/force_dem.c new file mode 100644 index 0000000..8162481 --- /dev/null +++ b/lammps/force_dem.c @@ -0,0 +1,232 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2021 RRZE, University Erlangen-Nuremberg + * + * This file is part of MD-Bench. + * + * MD-Bench is free software: you can redistribute it and/or modify it + * under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A + * PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License along + * with MD-Bench. If not, see . + * ======================================================================================= + */ +#include + +#include +#include +#include +#include +#include + +// TODO: Joint common files for gromacs and lammps variants +#include "../gromacs/includes/simd.h" + +double computeForceDEMFullNeigh_plain_c(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) { + int Nlocal = atom->Nlocal; + int* neighs; +#ifndef EXPLICIT_TYPES + MD_FLOAT cutforcesq = param->cutforce * param->cutforce; +#endif + + for(int i = 0; i < Nlocal; i++) { + atom_fx(i) = 0.0; + atom_fy(i) = 0.0; + atom_fz(i) = 0.0; + } + + double S = getTimeStamp(); + LIKWID_MARKER_START("force"); + +#pragma omp parallel for + for(int i = 0; i < Nlocal; i++) { + neighs = &neighbor->neighbors[i * neighbor->maxneighs]; + int numneighs = neighbor->numneigh[i]; + MD_FLOAT irad = atom->radius[i]; + MD_FLOAT xtmp = atom_x(i); + MD_FLOAT ytmp = atom_y(i); + MD_FLOAT ztmp = atom_z(i); + MD_FLOAT fix = 0; + MD_FLOAT fiy = 0; + MD_FLOAT fiz = 0; + +#ifdef EXPLICIT_TYPES + const int type_i = atom->type[i]; +#endif + + for(int k = 0; k < numneighs; k++) { + int j = neighs[k]; + MD_FLOAT jrad = atom->radius[j]; + MD_FLOAT xj = atom_x(j); + MD_FLOAT yj = atom_y(j); + MD_FLOAT zj = atom_z(j); + MD_FLOAT delx = xtmp - xj; + MD_FLOAT dely = ytmp - yj; + MD_FLOAT delz = ztmp - zj; + 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]; +#endif + + if(rsq < cutforcesq) { + MD_FLOAT sq = sqrt(rsq); + // penetration depth + MD_FLOAT p = irad + jrad - sq; + if(p < 0) { continue; } + + // contact position + MD_FLOAT cterm = jrad / (irad + jrad); + MD_FLOAT cx = xj + cterm * delx; + MD_FLOAT cy = yj + cterm * dely; + MD_FLOAT cz = zj + cterm * delz; + + // delta contact and particle position + MD_FLOAT delcx = cx - xtmp; + MD_FLOAT delcy = cy - ytmp; + MD_FLOAT delcz = cz - ztmp; + + // contact velocity + MD_FLOAT cvx = (atom_vx(i) + atom_avx(i) * delcx) - (atom_vx(j) + atom_avx(j) * (cx - xj)); + MD_FLOAT cvy = (atom_vy(i) + atom_avy(i) * delcy) - (atom_vy(j) + atom_avy(j) * (cy - yj)); + MD_FLOAT cvz = (atom_vz(i) + atom_avz(i) * delcz) - (atom_vz(j) + atom_avz(j) * (cz - zj)); + MD_FLOAT vsq = sqrt(cvx * cvx + cvy * cvy + cvz * cvz); + + // normal distance + MD_FLOAT nx = delx / sq; + MD_FLOAT ny = dely / sq; + MD_FLOAT nz = delz / sq; + + // normal contact velocity + MD_FLOAT vnx = cvx / vsq; + MD_FLOAT vny = cvy / vsq; + MD_FLOAT vnz = cvz / vsq; + + // forces + MD_FLOAT fnx = ks * p * nx - kdn * vn; + MD_FLOAT fny = ks * p * ny - kdn * vn; + MD_FLOAT fnz = ks * p * nz - kdn * vn; + MD_FLOAT ftx = MIN(kdt * vtsq, kf * fnx) * tx; + MD_FLOAT fty = MIN(kdt * vtsq, kf * fny) * ty; + MD_FLOAT ftz = MIN(kdt * vtsq, kf * fnz) * tz; + fix += fnx + ftx; + fiy += fny + fty; + fiz += fnz + ftz; + + // torque + //MD_FLOAT taux = delcx * ftx; + //MD_FLOAT tauy = delcy * fty; + //MD_FLOAT tauz = delcz * ftz; + +#ifdef USE_REFERENCE_VERSION + addStat(stats->atoms_within_cutoff, 1); + } else { + addStat(stats->atoms_outside_cutoff, 1); +#endif + } + } + + atom_fx(i) += fix; + atom_fy(i) += fiy; + atom_fz(i) += fiz; + + addStat(stats->total_force_neighs, numneighs); + addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH); + } + + LIKWID_MARKER_STOP("force"); + double E = getTimeStamp(); + return E-S; +} + +double computeForceDEMHalfNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) { + int Nlocal = atom->Nlocal; + int* neighs; +#ifndef EXPLICIT_TYPES + MD_FLOAT cutforcesq = param->cutforce * param->cutforce; + MD_FLOAT sigma6 = param->sigma6; + MD_FLOAT epsilon = param->epsilon; +#endif + + for(int i = 0; i < Nlocal; i++) { + atom_fx(i) = 0.0; + atom_fy(i) = 0.0; + atom_fz(i) = 0.0; + } + + double S = getTimeStamp(); + LIKWID_MARKER_START("forceLJ-halfneigh"); + + for(int i = 0; i < Nlocal; i++) { + neighs = &neighbor->neighbors[i * neighbor->maxneighs]; + int numneighs = neighbor->numneigh[i]; + MD_FLOAT xtmp = atom_x(i); + MD_FLOAT ytmp = atom_y(i); + MD_FLOAT ztmp = atom_z(i); + MD_FLOAT fix = 0; + MD_FLOAT fiy = 0; + MD_FLOAT fiz = 0; + +#ifdef EXPLICIT_TYPES + const int type_i = atom->type[i]; +#endif + + // Pragma required to vectorize the inner loop +#ifdef ENABLE_OMP_SIMD + #pragma omp simd reduction(+: fix,fiy,fiz) +#endif + 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; + +#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) { + MD_FLOAT sr2 = 1.0 / rsq; + MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6; + MD_FLOAT force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon; + fix += delx * force; + fiy += dely * force; + fiz += delz * force; + + // We do not need to update forces for ghost atoms + if(j < Nlocal) { + atom_fx(j) -= delx * force; + atom_fy(j) -= dely * force; + atom_fz(j) -= delz * force; + } + } + } + + atom_fx(i) += fix; + atom_fy(i) += fiy; + atom_fz(i) += fiz; + + addStat(stats->total_force_neighs, numneighs); + addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH); + } + + LIKWID_MARKER_STOP("forceLJ-halfneigh"); + double E = getTimeStamp(); + return E-S; +} diff --git a/lammps/includes/atom.h b/lammps/includes/atom.h index cf161c9..50b15a1 100644 --- a/lammps/includes/atom.h +++ b/lammps/includes/atom.h @@ -37,6 +37,10 @@ typedef struct { MD_FLOAT *sigma6; MD_FLOAT *cutforcesq; MD_FLOAT *cutneighsq; + // DEM + MD_FLOAT radius; + MD_FLOAT *av; + MD_FLOAT *r; } Atom; extern void initAtom(Atom*);