From 577955dfb7acfc741a6a62de7042f069594d2b42 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Wed, 13 Jul 2022 02:34:33 +0200 Subject: [PATCH] Apply first changes to DEM kernel Signed-off-by: Rafael Ravedutti --- lammps/force_dem.c | 140 ++++++++------------------------------------- lammps/main.c | 4 +- lammps/parameter.c | 3 + 3 files changed, 28 insertions(+), 119 deletions(-) diff --git a/lammps/force_dem.c b/lammps/force_dem.c index 698c587..47b5564 100644 --- a/lammps/force_dem.c +++ b/lammps/force_dem.c @@ -37,9 +37,7 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor int* neighs; MD_FLOAT k_s = param->k_s; MD_FLOAT k_dn = param->k_dn; -#ifndef EXPLICIT_TYPES MD_FLOAT cutforcesq = param->cutforce * param->cutforce; -#endif for(int i = 0; i < Nlocal; i++) { atom_fx(i) = 0.0; @@ -50,7 +48,6 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor 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]; @@ -62,10 +59,6 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor 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]; @@ -77,32 +70,11 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor 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 r = sqrt(rsq); - // penetration depth MD_FLOAT p = irad + jrad - r; - if(p >= 0) { - // 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)); + if(p > 0) { MD_FLOAT delvx = atom_vx(i) - atom_vx(j); MD_FLOAT delvy = atom_vy(i) - atom_vy(j); MD_FLOAT delvz = atom_vz(i) - atom_vz(j); @@ -119,9 +91,28 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor MD_FLOAT nvz = delvz / vr; // forces - fix += k_s * p * nx - k_dn * nvx; - fiy += k_s * p * ny - k_dn * nvy; - fiz += k_s * p * nz - k_dn * nvz; + atom_fx(i) += k_s * p * nx - k_dn * nvx; + atom_fy(i) += k_s * p * ny - k_dn * nvy; + atom_fz(i) += k_s * p * nz - k_dn * nvz; + atom_fx(j) += -k_s * p * nx - k_dn * nvx; + atom_fy(j) += -k_s * p * ny - k_dn * nvy; + atom_fz(j) += -k_s * p * nz - k_dn * nvz; + + // 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)); // tangential force //fix += MIN(kdt * vtsq, kf * fnx) * tx; @@ -140,10 +131,6 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor } } - 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); } @@ -152,84 +139,3 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor 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/main.c b/lammps/main.c index 03760bf..c372470 100644 --- a/lammps/main.c +++ b/lammps/main.c @@ -50,7 +50,6 @@ extern double computeForceLJFullNeigh_simd(Parameter*, Atom*, Neighbor*, Stats*) extern double computeForceLJHalfNeigh(Parameter*, Atom*, Neighbor*, Stats*); extern double computeForceEam(Eam*, Parameter*, Atom*, Neighbor*, Stats*); extern double computeForceDemFullNeigh(Parameter*, Atom*, Neighbor*, Stats*); -extern double computeForceDemHalfNeigh(Parameter*, Atom*, Neighbor*, Stats*); #ifdef USE_SIMD_KERNEL # define KERNEL_NAME "SIMD" @@ -135,7 +134,8 @@ double computeForce(Eam *eam, Parameter *param, Atom *atom, Neighbor *neighbor, return computeForceEam(eam, param, atom, neighbor, stats); } else if(param->force_field == FF_DEM) { if(param->half_neigh) { - return computeForceDemHalfNeigh(param, atom, neighbor, stats); + fprintf(stderr, "Error: DEM cannot use half neighbor-lists!\n"); + return 0.0; } else { return computeForceDemFullNeigh(param, atom, neighbor, stats); } diff --git a/lammps/parameter.c b/lammps/parameter.c index faeb7ba..5606bb2 100644 --- a/lammps/parameter.c +++ b/lammps/parameter.c @@ -130,6 +130,9 @@ void readParameter(Parameter *param, const char *filename) { } } + // Update dtforce + param->dtforce = 0.5 * param->dt; + // Update sigma6 parameter MD_FLOAT s2 = param->sigma * param->sigma; param->sigma6 = s2 * s2 * s2;