From 79483a446eb5546e1de8f08b82b521fdf8ae2382 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Wed, 6 Jul 2022 01:07:39 +0200 Subject: [PATCH] Adjust code with DEM to be compilable Signed-off-by: Rafael Ravedutti --- lammps/atom.c | 34 ++++++++------- lammps/force_dem.c | 85 +++++++++++++++++++------------------ lammps/includes/atom.h | 8 +++- lammps/includes/parameter.h | 4 ++ lammps/main-stub.c | 6 +-- lammps/main.c | 23 ++++------ lammps/neighbor.c | 19 ++++++--- lammps/parameter.c | 6 +++ lammps/thermo.c | 27 +++++------- 9 files changed, 116 insertions(+), 96 deletions(-) diff --git a/lammps/atom.c b/lammps/atom.c index d923b0b..46e4110 100644 --- a/lammps/atom.c +++ b/lammps/atom.c @@ -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 diff --git a/lammps/force_dem.c b/lammps/force_dem.c index 8162481..6ba1c70 100644 --- a/lammps/force_dem.c +++ b/lammps/force_dem.c @@ -20,6 +20,7 @@ * with MD-Bench. If not, see . * ======================================================================================= */ +#include #include #include @@ -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; - // 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; - // 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 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); - // 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 / r; + MD_FLOAT ny = dely / r; + MD_FLOAT nz = delz / r; - // normal distance - MD_FLOAT nx = delx / sq; - MD_FLOAT ny = dely / sq; - MD_FLOAT nz = delz / sq; + // normal contact velocity + MD_FLOAT nvx = delvx / vr; + MD_FLOAT nvy = delvy / vr; + MD_FLOAT nvz = delvz / vr; - // 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; + // 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; + // 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 { diff --git a/lammps/includes/atom.h b/lammps/includes/atom.h index 50b15a1..4c473d8 100644 --- a/lammps/includes/atom.h +++ b/lammps/includes/atom.h @@ -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] diff --git a/lammps/includes/parameter.h b/lammps/includes/parameter.h index 5c8eba7..497f4ce 100644 --- a/lammps/includes/parameter.h +++ b/lammps/includes/parameter.h @@ -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*); diff --git a/lammps/main-stub.c b/lammps/main-stub.c index 2e8f1ff..eb6aeee 100644 --- a/lammps/main-stub.c +++ b/lammps/main-stub.c @@ -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++; } diff --git a/lammps/main.c b/lammps/main.c index f770fa1..3551ce1 100644 --- a/lammps/main.c +++ b/lammps/main.c @@ -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); } } diff --git a/lammps/neighbor.c b/lammps/neighbor.c index 8900a04..826a00c 100644 --- a/lammps/neighbor.c +++ b/lammps/neighbor.c @@ -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; } diff --git a/lammps/parameter.c b/lammps/parameter.c index 90fc3cf..e16f538 100644 --- a/lammps/parameter.c +++ b/lammps/parameter.c @@ -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) { diff --git a/lammps/thermo.c b/lammps/thermo.c index 73e94ae..0fbbf41 100644 --- a/lammps/thermo.c +++ b/lammps/thermo.c @@ -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; } }