moebiusband/MD-Bench
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Adjust code with DEM to be compilable

Browse Source 
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2022-07-06 01:07:39 +02:00
parent bb599c9ea8
commit 79483a446e
9 changed files with 116 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
lammps/atom.c
View File

@ -55,7 +55,7 @@ void initAtom(Atom *atom) {
atom->sigma6 = NULL;
atom->cutforcesq = NULL;
atom->cutneighsq = NULL;
atom->radius = 1.0;
atom->radius = NULL;
atom->av = NULL;
atom->r = NULL;
}
@ -144,9 +144,9 @@ void createAtom(Atom *atom, Parameter *param) {
atom_x(atom->Nlocal) = xtmp;
atom_y(atom->Nlocal) = ytmp;
atom_z(atom->Nlocal) = ztmp;
atom->vx[atom->Nlocal] = vxtmp;
atom->vy[atom->Nlocal] = vytmp;
atom->vz[atom->Nlocal] = vztmp;
atom_vx(atom->Nlocal) = vxtmp;
atom_vy(atom->Nlocal) = vytmp;
atom_vz(atom->Nlocal) = vztmp;
atom->type[atom->Nlocal] = rand() % atom->ntypes;
atom->Nlocal++;
}
@ -220,9 +220,9 @@ int readAtom_pdb(Atom* atom, Parameter* param) {
atom_x(atom_id) = atof(strtok(NULL, " "));
atom_y(atom_id) = atof(strtok(NULL, " "));
atom_z(atom_id) = atof(strtok(NULL, " "));
atom->vx[atom_id] = 0.0;
atom->vy[atom_id] = 0.0;
atom->vz[atom_id] = 0.0;
atom_vx(atom_id) = 0.0;
atom_vy(atom_id) = 0.0;
atom_vz(atom_id) = 0.0;
occupancy = atof(strtok(NULL, " "));
charge = atof(strtok(NULL, " "));
atom->ntypes = MAX(atom->type[atom_id] + 1, atom->ntypes);
@ -299,9 +299,9 @@ int readAtom_gro(Atom* atom, Parameter* param) {
atom_x(atom_id) = atof(strtok(NULL, " "));
atom_y(atom_id) = atof(strtok(NULL, " "));
atom_z(atom_id) = atof(strtok(NULL, " "));
atom->vx[atom_id] = atof(strtok(NULL, " "));
atom->vy[atom_id] = atof(strtok(NULL, " "));
atom->vz[atom_id] = atof(strtok(NULL, " "));
atom_vx(atom_id) = atof(strtok(NULL, " "));
atom_vy(atom_id) = atof(strtok(NULL, " "));
atom_vz(atom_id) = atof(strtok(NULL, " "));
atom->ntypes = MAX(atom->type[atom_id] + 1, atom->ntypes);
atom->Natoms++;
atom->Nlocal++;
@ -396,9 +396,9 @@ int readAtom_dmp(Atom* atom, Parameter* param) {
atom_x(atom_id) = atof(strtok(NULL, " "));
atom_y(atom_id) = atof(strtok(NULL, " "));
atom_z(atom_id) = atof(strtok(NULL, " "));
atom->vx[atom_id] = atof(strtok(NULL, " "));
atom->vy[atom_id] = atof(strtok(NULL, " "));
atom->vz[atom_id] = atof(strtok(NULL, " "));
atom_vx(atom_id) = atof(strtok(NULL, " "));
atom_vy(atom_id) = atof(strtok(NULL, " "));
atom_vz(atom_id) = atof(strtok(NULL, " "));
atom->ntypes = MAX(atom->type[atom_id], atom->ntypes);
read_atoms++;
}
@ -441,20 +441,22 @@ void growAtom(Atom *atom) {
#ifdef AOS
atom->x = (MD_FLOAT*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
#else
atom->x = (MD_FLOAT*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->y = (MD_FLOAT*) reallocate(atom->y, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->z = (MD_FLOAT*) reallocate(atom->z, 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->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));
#endif
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->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->radius = (int *) reallocate(atom->radius, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
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

65
lammps/force_dem.c
View File

@ -20,6 +20,7 @@
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
* =======================================================================================
*/
#include <math.h>
#include <likwid-marker.h>
#include <timing.h>
@ -34,6 +35,8 @@
double computeForceDEMFullNeigh_plain_c(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
int Nlocal = atom->Nlocal;
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
@ -81,54 +84,54 @@ double computeForceDEMFullNeigh_plain_c(Parameter *param, Atom *atom, Neighbor *
#endif
if(rsq < cutforcesq) {
MD_FLOAT sq = sqrt(rsq);
MD_FLOAT r = sqrt(rsq);
// penetration depth
MD_FLOAT p = irad + jrad - sq;
if(p < 0) { continue; }
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;
//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;
//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);
//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 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);
MD_FLOAT vr = sqrt(delvx * delvx + delvy * delvy + delvz * delvz);
// normal distance
MD_FLOAT nx = delx / sq;
MD_FLOAT ny = dely / sq;
MD_FLOAT nz = delz / sq;
MD_FLOAT nx = delx / r;
MD_FLOAT ny = dely / r;
MD_FLOAT nz = delz / r;
// normal contact velocity
MD_FLOAT vnx = cvx / vsq;
MD_FLOAT vny = cvy / vsq;
MD_FLOAT vnz = cvz / vsq;
MD_FLOAT nvx = delvx / vr;
MD_FLOAT nvy = delvy / vr;
MD_FLOAT nvz = delvz / vr;
// 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;
fix += k_s * p * nx - k_dn * nvx;
fiy += k_s * p * ny - k_dn * nvy;
fiz += k_s * p * nz - k_dn * nvz;
// tangential force
//fix += MIN(kdt * vtsq, kf * fnx) * tx;
//fiy += MIN(kdt * vtsq, kf * fny) * ty;
//fiz += MIN(kdt * vtsq, kf * fnz) * tz;
// 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 {

8
lammps/includes/atom.h
View File

@ -38,7 +38,7 @@ typedef struct {
MD_FLOAT *cutforcesq;
MD_FLOAT *cutneighsq;
// DEM
MD_FLOAT radius;
MD_FLOAT *radius;
MD_FLOAT *av;
MD_FLOAT *r;
} Atom;
@ -56,6 +56,9 @@ extern void growAtom(Atom*);
#define atom_x(i) atom->x[(i) * 3 + 0]
#define atom_y(i) atom->x[(i) * 3 + 1]
#define atom_z(i) atom->x[(i) * 3 + 2]
#define atom_vx(i) atom->vx[(i) * 3 + 0]
#define atom_vy(i) atom->vx[(i) * 3 + 1]
#define atom_vz(i) atom->vx[(i) * 3 + 2]
#define atom_fx(i) atom->fx[(i) * 3 + 0]
#define atom_fy(i) atom->fx[(i) * 3 + 1]
#define atom_fz(i) atom->fx[(i) * 3 + 2]
@ -64,6 +67,9 @@ extern void growAtom(Atom*);
#define atom_x(i) atom->x[i]
#define atom_y(i) atom->y[i]
#define atom_z(i) atom->z[i]
#define atom_vx(i) atom->vx[i]
#define atom_vy(i) atom->vy[i]
#define atom_vz(i) atom->vz[i]
#define atom_fx(i) atom->fx[i]
#define atom_fy(i) atom->fy[i]
#define atom_fz(i) atom->fz[i]

4
lammps/includes/parameter.h
View File

@ -58,6 +58,10 @@ typedef struct {
MD_FLOAT xprd, yprd, zprd;
double proc_freq;
char* eam_file;
// DEM
MD_FLOAT k_s;
MD_FLOAT k_dn;
MD_FLOAT gx, gy, gz;
} Parameter;
void initParameter(Parameter*);

6
lammps/main-stub.c
View File

@ -227,9 +227,9 @@ int main(int argc, const char *argv[]) {
atom_x(atom->Nlocal) = (MD_FLOAT)(i) * 0.00001;
atom_y(atom->Nlocal) = (MD_FLOAT)(i) * 0.00001;
atom_z(atom->Nlocal) = (MD_FLOAT)(i) * 0.00001;
atom->vx[atom->Nlocal] = 0.0;
atom->vy[atom->Nlocal] = 0.0;
atom->vz[atom->Nlocal] = 0.0;
atom_vx(atom->Nlocal) = 0.0;
atom_vy(atom->Nlocal) = 0.0;
atom_vz(atom->Nlocal) = 0.0;
atom->Nlocal++;
}

23
lammps/main.c
View File

@ -102,26 +102,21 @@ double reneighbour(Parameter *param, Atom *atom, Neighbor *neighbor) {
}
void initialIntegrate(Parameter *param, Atom *atom) {
MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
for(int i = 0; i < atom->Nlocal; i++) {
vx[i] += param->dtforce * atom_fx(i);
vy[i] += param->dtforce * atom_fy(i);
vz[i] += param->dtforce * atom_fz(i);
atom_x(i) = atom_x(i) + param->dt * vx[i];
atom_y(i) = atom_y(i) + param->dt * vy[i];
atom_z(i) = atom_z(i) + param->dt * vz[i];
atom_vx(i) += param->dtforce * atom_fx(i);
atom_vy(i) += param->dtforce * atom_fy(i);
atom_vz(i) += param->dtforce * atom_fz(i);
atom_x(i) = atom_x(i) + param->dt * atom_vx(i);
atom_y(i) = atom_y(i) + param->dt * atom_vy(i);
atom_z(i) = atom_z(i) + param->dt * atom_vz(i);
}
}
void finalIntegrate(Parameter *param, Atom *atom) {
MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
for(int i = 0; i < atom->Nlocal; i++) {
vx[i] += param->dtforce * atom_fx(i);
vy[i] += param->dtforce * atom_fy(i);
vz[i] += param->dtforce * atom_fz(i);
atom_vx(i) += param->dtforce * atom_fx(i);
atom_vy(i) += param->dtforce * atom_fy(i);
atom_vz(i) += param->dtforce * atom_fz(i);
}
}

19
lammps/neighbor.c
View File

@ -349,14 +349,15 @@ void sortAtom(Atom* atom) {
#ifdef AOS
MD_FLOAT* new_x = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT) * 3);
MD_FLOAT* new_vx = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT) * 3);
#else
MD_FLOAT* new_x = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_y = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_z = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
#endif
MD_FLOAT* new_vx = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_vy = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_vz = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
#endif
MD_FLOAT* old_x = atom->x; MD_FLOAT* old_y = atom->y; MD_FLOAT* old_z = atom->z;
MD_FLOAT* old_vx = atom->vx; MD_FLOAT* old_vy = atom->vy; MD_FLOAT* old_vz = atom->vz;
@ -370,24 +371,32 @@ void sortAtom(Atom* atom) {
new_x[new_i * 3 + 0] = old_x[old_i * 3 + 0];
new_x[new_i * 3 + 1] = old_x[old_i * 3 + 1];
new_x[new_i * 3 + 2] = old_x[old_i * 3 + 2];
new_vx[new_i * 3 + 0] = old_vx[old_i * 3 + 0];
new_vx[new_i * 3 + 1] = old_vx[old_i * 3 + 1];
new_vx[new_i * 3 + 2] = old_vx[old_i * 3 + 2];
#else
new_x[new_i] = old_x[old_i];
new_y[new_i] = old_y[old_i];
new_z[new_i] = old_z[old_i];
#endif
new_vx[new_i] = old_vx[old_i];
new_vy[new_i] = old_vy[old_i];
new_vz[new_i] = old_vz[old_i];
#endif
}
}
free(atom->x);
free(atom->vx);
atom->x = new_x;
atom->vx = new_vx;
#ifndef AOS
free(atom->y);
free(atom->z);
atom->y = new_y; atom->z = new_z;
free(atom->vy);
free(atom->vz);
atom->y = new_y;
atom->z = new_z;
atom->vy = new_vy;
atom->vz = new_vz;
#endif
free(atom->vx); free(atom->vy); free(atom->vz);
atom->vx = new_vx; atom->vy = new_vy; atom->vz = new_vz;
}

6
lammps/parameter.c
View File

@ -54,6 +54,12 @@ void initParameter(Parameter *param) {
param->v_out_every = 5;
param->half_neigh = 0;
param->proc_freq = 2.4;
// DEM
param->k_s = 1.0;
param->k_dn = 1.0;
param->gx = 0.0;
param->gy = 0.0;
param->gz = 0.0;
}
void readParameter(Parameter *param, const char *filename) {

27
lammps/thermo.c
View File

@ -70,12 +70,8 @@ void setupThermo(Parameter *param, int natoms)
void computeThermo(int iflag, Parameter *param, Atom *atom)
{
MD_FLOAT t = 0.0, p;
MD_FLOAT* vx = atom->vx;
MD_FLOAT* vy = atom->vy;
MD_FLOAT* vz = atom->vz;
for(int i = 0; i < atom->Nlocal; i++) {
t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass;
t += (atom_vx(i) * atom_vx(i) + atom_vy(i) * atom_vy(i) + atom_vz(i) * atom_vz(i)) * param->mass;
}
t = t * t_scale;
@ -100,12 +96,11 @@ void adjustThermo(Parameter *param, Atom *atom)
{
/* zero center-of-mass motion */
MD_FLOAT vxtot = 0.0; MD_FLOAT vytot = 0.0; MD_FLOAT vztot = 0.0;
MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
for(int i = 0; i < atom->Nlocal; i++) {
vxtot += vx[i];
vytot += vy[i];
vztot += vz[i];
vxtot += atom_vx(i);
vytot += atom_vy(i);
vztot += atom_vz(i);
}
vxtot = vxtot / atom->Natoms;
@ -113,24 +108,24 @@ void adjustThermo(Parameter *param, Atom *atom)
vztot = vztot / atom->Natoms;
for(int i = 0; i < atom->Nlocal; i++) {
vx[i] -= vxtot;
vy[i] -= vytot;
vz[i] -= vztot;
atom_vx(i) -= vxtot;
atom_vy(i) -= vytot;
atom_vz(i) -= vztot;
}
t_act = 0;
MD_FLOAT t = 0.0;
for(int i = 0; i < atom->Nlocal; i++) {
t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass;
t += (atom_vx(i) * atom_vx(i) + atom_vy(i) * atom_vy(i) + atom_vz(i) * atom_vz(i)) * param->mass;
}
t *= t_scale;
MD_FLOAT factor = sqrt(param->temp / t);
for(int i = 0; i < atom->Nlocal; i++) {
vx[i] *= factor;
vy[i] *= factor;
vz[i] *= factor;
atom_vx(i) *= factor;
atom_vy(i) *= factor;
atom_vz(i) *= factor;
}
}