ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

more refactoring, modernizing, removing of debug code

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-31 06:24:44 -04:00
parent b64e90dfe0
commit b97c30e6e7
7 changed files with 1920 additions and 2250 deletions
Download Patch File Download Diff File Expand all files Collapse all files

720
src/MACHDYN/fix_smd_tlsph_reference_configuration.cpp
View File

@ -27,7 +27,6 @@
#include "fix_smd_tlsph_reference_configuration.h"
#include <Eigen/Eigen>
#include "atom.h"
#include "comm.h"
#include "neigh_list.h"
@ -41,12 +40,14 @@
#include "smd_kernels.h"
#include "smd_math.h"
#include <Eigen/Eigen>
using namespace Eigen;
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace SMD_Kernels;
using namespace std;
using namespace SMD_Math;
#define DELTA 16384
#define INSERT_PREDEFINED_CRACKS false
@ -54,141 +55,134 @@ using namespace SMD_Math;
/* ---------------------------------------------------------------------- */
FixSMD_TLSPH_ReferenceConfiguration::FixSMD_TLSPH_ReferenceConfiguration(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
Fix(lmp, narg, arg) {
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR, "Pair tlsph with partner list requires an atom map, see atom_modify");
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR, "Pair tlsph with partner list requires an atom map, see atom_modify");
maxpartner = 1;
npartner = nullptr;
partner = nullptr;
wfd_list = nullptr;
wf_list = nullptr;
energy_per_bond = nullptr;
degradation_ij = nullptr;
grow_arrays(atom->nmax);
atom->add_callback(Atom::GROW);
maxpartner = 1;
npartner = nullptr;
partner = nullptr;
wfd_list = nullptr;
wf_list = nullptr;
energy_per_bond = nullptr;
degradation_ij = nullptr;
grow_arrays(atom->nmax);
atom->add_callback(Atom::GROW);
// initialize npartner to 0 so neighbor list creation is OK the 1st time
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
npartner[i] = 0;
}
// initialize npartner to 0 so neighbor list creation is OK the 1st time
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
npartner[i] = 0;
}
comm_forward = 14;
updateFlag = 1;
comm_forward = 14;
updateFlag = 1;
}
/* ---------------------------------------------------------------------- */
FixSMD_TLSPH_ReferenceConfiguration::~FixSMD_TLSPH_ReferenceConfiguration() {
// unregister this fix so atom class doesn't invoke it any more
// unregister this fix so atom class doesn't invoke it any more
atom->delete_callback(id,Atom::GROW);
atom->delete_callback(id,Atom::GROW);
// delete locally stored arrays
memory->destroy(npartner);
memory->destroy(partner);
memory->destroy(wfd_list);
memory->destroy(wf_list);
memory->destroy(degradation_ij);
memory->destroy(energy_per_bond);
memory->destroy(npartner);
memory->destroy(partner);
memory->destroy(wfd_list);
memory->destroy(wf_list);
memory->destroy(degradation_ij);
memory->destroy(energy_per_bond);
}
/* ---------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::setmask() {
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::init() {
if (atom->tag_enable == 0)
error->all(FLERR, "Pair style tlsph requires atoms have IDs");
if (atom->tag_enable == 0) error->all(FLERR, "Pair style tlsph requires atoms have IDs");
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::pre_exchange() {
//return;
//printf("in FixSMD_TLSPH_ReferenceConfiguration::pre_exchange()\n");
double **defgrad = atom->smd_data_9;
double *radius = atom->radius;
double *rho = atom->rho;
double *vfrac = atom->vfrac;
double **x = atom->x;
double **x0 = atom->x0;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int i, itmp;
int *mask = atom->mask;
if (igroup == atom->firstgroup) {
nlocal = atom->nfirst;
double **defgrad = atom->smd_data_9;
double *radius = atom->radius;
double *rho = atom->rho;
double *vfrac = atom->vfrac;
double **x = atom->x;
double **x0 = atom->x0;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int i, itmp;
int *mask = atom->mask;
if (igroup == atom->firstgroup) {
nlocal = atom->nfirst;
}
int *updateFlag_ptr = (int *) force->pair->extract("smd/tlsph/updateFlag_ptr", itmp);
if (updateFlag_ptr == nullptr) {
error->one(FLERR,
"fix FixSMD_TLSPH_ReferenceConfiguration failed to access updateFlag pointer. Check if a pair style exist which calculates this quantity.");
}
int *nn = (int *) force->pair->extract("smd/tlsph/numNeighsRefConfig_ptr", itmp);
if (nn == nullptr) {
error->all(FLERR, "FixSMDIntegrateTlsph::updateReferenceConfiguration() failed to access numNeighsRefConfig_ptr array");
}
// sum all update flag across processors
MPI_Allreduce(updateFlag_ptr, &updateFlag, 1, MPI_INT, MPI_MAX, world);
if (updateFlag > 0) {
if (comm->me == 0) utils::logmesg(lmp, "**** updating ref config at step: {}\n", update->ntimestep);
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
// re-set x0 coordinates
x0[i][0] = x[i][0];
x0[i][1] = x[i][1];
x0[i][2] = x[i][2];
// re-set deformation gradient
defgrad[i][0] = 1.0;
defgrad[i][1] = 0.0;
defgrad[i][2] = 0.0;
defgrad[i][3] = 0.0;
defgrad[i][4] = 1.0;
defgrad[i][5] = 0.0;
defgrad[i][6] = 0.0;
defgrad[i][7] = 0.0;
defgrad[i][8] = 1.0;
/*
* Adjust particle volume as the reference configuration is changed.
* We safeguard against excessive deformations by limiting the adjustment range
* to the intervale J \in [0.9..1.1]
*/
vfrac[i] = rmass[i] / rho[i];
if (nn[i] < 15) {
radius[i] *= 1.2;
}
}
}
int *updateFlag_ptr = (int *) force->pair->extract("smd/tlsph/updateFlag_ptr", itmp);
if (updateFlag_ptr == nullptr) {
error->one(FLERR,
"fix FixSMD_TLSPH_ReferenceConfiguration failed to access updateFlag pointer. Check if a pair style exist which calculates this quantity.");
}
// update of reference config could have changed x0, vfrac, radius
// communicate these quantities now to ghosts: x0, vfrac, radius
comm->forward_comm(this);
int *nn = (int *) force->pair->extract("smd/tlsph/numNeighsRefConfig_ptr", itmp);
if (nn == nullptr) {
error->all(FLERR, "FixSMDIntegrateTlsph::updateReferenceConfiguration() failed to access numNeighsRefConfig_ptr array");
}
// sum all update flag across processors
MPI_Allreduce(updateFlag_ptr, &updateFlag, 1, MPI_INT, MPI_MAX, world);
if (updateFlag > 0) {
if (comm->me == 0) {
printf("**** updating ref config at step: " BIGINT_FORMAT "\n", update->ntimestep);
}
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
// re-set x0 coordinates
x0[i][0] = x[i][0];
x0[i][1] = x[i][1];
x0[i][2] = x[i][2];
// re-set deformation gradient
defgrad[i][0] = 1.0;
defgrad[i][1] = 0.0;
defgrad[i][2] = 0.0;
defgrad[i][3] = 0.0;
defgrad[i][4] = 1.0;
defgrad[i][5] = 0.0;
defgrad[i][6] = 0.0;
defgrad[i][7] = 0.0;
defgrad[i][8] = 1.0;
/*
* Adjust particle volume as the reference configuration is changed.
* We safeguard against excessive deformations by limiting the adjustment range
* to the intervale J \in [0.9..1.1]
*/
vfrac[i] = rmass[i] / rho[i];
//
if (nn[i] < 15) {
radius[i] *= 1.2;
} // else //{
// radius[i] *= pow(J, 1.0 / domain->dimension);
//}
}
}
// update of reference config could have changed x0, vfrac, radius
// communicate these quantities now to ghosts: x0, vfrac, radius
comm->forward_comm(this);
setup(0);
}
setup(0);
}
}
/* ----------------------------------------------------------------------
@ -197,155 +191,155 @@ void FixSMD_TLSPH_ReferenceConfiguration::pre_exchange() {
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::setup(int /*vflag*/) {
int i, j, ii, jj, n, inum, jnum;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, h, wf, wfd;
Vector3d dx;
int i, j, ii, jj, n, inum, jnum;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, h, wf, wfd;
Vector3d dx;
if (updateFlag == 0)
return;
if (updateFlag == 0)
return;
int nlocal = atom->nlocal;
nmax = atom->nmax;
grow_arrays(nmax);
int nlocal = atom->nlocal;
nmax = atom->nmax;
grow_arrays(nmax);
// 1st loop over neighbor list
// calculate npartner for each owned atom
// nlocal_neigh = nlocal when neigh list was built, may be smaller than nlocal
double **x0 = atom->x;
double *radius = atom->radius;
int *mask = atom->mask;
tagint *tag = atom->tag;
NeighList *list = pair->list;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double **x0 = atom->x;
double *radius = atom->radius;
int *mask = atom->mask;
tagint *tag = atom->tag;
NeighList *list = pair->list;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// zero npartner for all current atoms
for (i = 0; i < nlocal; i++)
npartner[i] = 0;
// zero npartner for all current atoms
for (i = 0; i < nlocal; i++)
npartner[i] = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
h = radius[i] + radius[j];
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
h = radius[i] + radius[j];
if (r <= h) {
npartner[i]++;
if (j < nlocal) {
npartner[j]++;
}
}
}
if (r <= h) {
npartner[i]++;
if (j < nlocal) {
npartner[j]++;
}
}
}
}
maxpartner = 0;
for (i = 0; i < nlocal; i++)
maxpartner = MAX(maxpartner, npartner[i]);
int maxall;
MPI_Allreduce(&maxpartner, &maxall, 1, MPI_INT, MPI_MAX, world);
maxpartner = maxall;
maxpartner = 0;
for (i = 0; i < nlocal; i++)
maxpartner = MAX(maxpartner, npartner[i]);
int maxall;
MPI_Allreduce(&maxpartner, &maxall, 1, MPI_INT, MPI_MAX, world);
maxpartner = maxall;
grow_arrays(nmax);
grow_arrays(nmax);
for (i = 0; i < nlocal; i++) {
npartner[i] = 0;
for (jj = 0; jj < maxpartner; jj++) {
wfd_list[i][jj] = 0.0;
wf_list[i][jj] = 0.0;
degradation_ij[i][jj] = 0.0;
energy_per_bond[i][jj] = 0.0;
}
for (i = 0; i < nlocal; i++) {
npartner[i] = 0;
for (jj = 0; jj < maxpartner; jj++) {
wfd_list[i][jj] = 0.0;
wf_list[i][jj] = 0.0;
degradation_ij[i][jj] = 0.0;
energy_per_bond[i][jj] = 0.0;
}
}
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
h = radius[i] + radius[j];
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
if (r < h) {
spiky_kernel_and_derivative(h, r, domain->dimension, wf, wfd);
partner[i][npartner[i]] = tag[j];
wfd_list[i][npartner[i]] = wfd;
wf_list[i][npartner[i]] = wf;
npartner[i]++;
if (j < nlocal) {
partner[j][npartner[j]] = tag[i];
wfd_list[j][npartner[j]] = wfd;
wf_list[j][npartner[j]] = wf;
npartner[j]++;
}
}
}
}
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
// count number of particles for which this group is active
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
// bond statistics
if (update->ntimestep > -1) {
n = 0;
int count = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
n += npartner[i];
count += 1;
}
}
int nall, countall;
MPI_Allreduce(&n, &nall, 1, MPI_INT, MPI_SUM, world);
MPI_Allreduce(&count, &countall, 1, MPI_INT, MPI_SUM, world);
if (countall < 1) countall = 1;
dx(0) = x0[i][0] - x0[j][0];
dx(1) = x0[i][1] - x0[j][1];
dx(2) = x0[i][2] - x0[j][2];
r = dx.norm();
h = radius[i] + radius[j];
if (comm->me == 0) {
if (screen) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
fprintf(screen, "TLSPH neighbors:\n");
fprintf(screen, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(screen, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
}
if (logfile) {
fprintf(logfile, "\nTLSPH neighbors:\n");
fprintf(logfile, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(logfile, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
}
}
}
if (INSERT_PREDEFINED_CRACKS) {
if (!crack_exclude(i, j))
continue;
}
if (r < h) {
spiky_kernel_and_derivative(h, r, domain->dimension, wf, wfd);
partner[i][npartner[i]] = tag[j];
wfd_list[i][npartner[i]] = wfd;
wf_list[i][npartner[i]] = wf;
npartner[i]++;
if (j < nlocal) {
partner[j][npartner[j]] = tag[i];
wfd_list[j][npartner[j]] = wfd;
wf_list[j][npartner[j]] = wf;
npartner[j]++;
}
}
}
}
// count number of particles for which this group is active
// bond statistics
if (update->ntimestep > -1) {
n = 0;
int count = 0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
n += npartner[i];
count += 1;
}
}
int nall, countall;
MPI_Allreduce(&n, &nall, 1, MPI_INT, MPI_SUM, world);
MPI_Allreduce(&count, &countall, 1, MPI_INT, MPI_SUM, world);
if (countall < 1) countall = 1;
if (comm->me == 0) {
if (screen) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
fprintf(screen, "TLSPH neighbors:\n");
fprintf(screen, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(screen, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
}
if (logfile) {
fprintf(logfile, "\nTLSPH neighbors:\n");
fprintf(logfile, " max # of neighbors for a single particle = %d\n", maxpartner);
fprintf(logfile, " average # of neighbors/particle in group tlsph = %g\n", (double) nall / countall);
}
}
}
updateFlag = 0; // set update flag to zero after the update
updateFlag = 0; // set update flag to zero after the update
}
@ -354,14 +348,14 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int /*vflag*/) {
------------------------------------------------------------------------- */
double FixSMD_TLSPH_ReferenceConfiguration::memory_usage() {
int nmax = atom->nmax;
int bytes = nmax * sizeof(int);
bytes += (double)nmax * maxpartner * sizeof(tagint); // partner array
bytes += (double)nmax * maxpartner * sizeof(float); // wf_list
bytes += (double)nmax * maxpartner * sizeof(float); // wfd_list
bytes += (double)nmax * maxpartner * sizeof(float); // damage_per_interaction array
bytes += (double)nmax * sizeof(int); // npartner array
return bytes;
int nmax = atom->nmax;
int bytes = nmax * sizeof(int);
bytes += (double)nmax * maxpartner * sizeof(tagint); // partner array
bytes += (double)nmax * maxpartner * sizeof(float); // wf_list
bytes += (double)nmax * maxpartner * sizeof(float); // wfd_list
bytes += (double)nmax * maxpartner * sizeof(float); // damage_per_interaction array
bytes += (double)nmax * sizeof(int); // npartner array
return bytes;
}
@ -370,13 +364,12 @@ double FixSMD_TLSPH_ReferenceConfiguration::memory_usage() {
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::grow_arrays(int nmax) {
//printf("in FixSMD_TLSPH_ReferenceConfiguration::grow_arrays\n");
memory->grow(npartner, nmax, "tlsph_refconfig_neigh:npartner");
memory->grow(partner, nmax, maxpartner, "tlsph_refconfig_neigh:partner");
memory->grow(wfd_list, nmax, maxpartner, "tlsph_refconfig_neigh:wfd");
memory->grow(wf_list, nmax, maxpartner, "tlsph_refconfig_neigh:wf");
memory->grow(degradation_ij, nmax, maxpartner, "tlsph_refconfig_neigh:degradation_ij");
memory->grow(energy_per_bond, nmax, maxpartner, "tlsph_refconfig_neigh:damage_onset_strain");
memory->grow(npartner, nmax, "tlsph_refconfig_neigh:npartner");
memory->grow(partner, nmax, maxpartner, "tlsph_refconfig_neigh:partner");
memory->grow(wfd_list, nmax, maxpartner, "tlsph_refconfig_neigh:wfd");
memory->grow(wf_list, nmax, maxpartner, "tlsph_refconfig_neigh:wf");
memory->grow(degradation_ij, nmax, maxpartner, "tlsph_refconfig_neigh:degradation_ij");
memory->grow(energy_per_bond, nmax, maxpartner, "tlsph_refconfig_neigh:damage_onset_strain");
}
/* ----------------------------------------------------------------------
@ -384,14 +377,14 @@ void FixSMD_TLSPH_ReferenceConfiguration::grow_arrays(int nmax) {
------------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::copy_arrays(int i, int j, int /*delflag*/) {
npartner[j] = npartner[i];
for (int m = 0; m < npartner[j]; m++) {
partner[j][m] = partner[i][m];
wfd_list[j][m] = wfd_list[i][m];
wf_list[j][m] = wf_list[i][m];
degradation_ij[j][m] = degradation_ij[i][m];
energy_per_bond[j][m] = energy_per_bond[i][m];
}
npartner[j] = npartner[i];
for (int m = 0; m < npartner[j]; m++) {
partner[j][m] = partner[i][m];
wfd_list[j][m] = wfd_list[i][m];
wf_list[j][m] = wf_list[i][m];
degradation_ij[j][m] = degradation_ij[i][m];
energy_per_bond[j][m] = energy_per_bond[i][m];
}
}
/* ----------------------------------------------------------------------
@ -402,18 +395,16 @@ int FixSMD_TLSPH_ReferenceConfiguration::pack_exchange(int i, double *buf) {
// NOTE: how do I know comm buf is big enough if extreme # of touching neighs
// Comm::BUFEXTRA may need to be increased
//printf("pack_exchange ...\n");
int m = 0;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
int m = 0;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
}
@ -422,27 +413,25 @@ int FixSMD_TLSPH_ReferenceConfiguration::pack_exchange(int i, double *buf) {
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::unpack_exchange(int nlocal, double *buf) {
if (nlocal == nmax) {
//printf("nlocal=%d, nmax=%d\n", nlocal, nmax);
nmax = nmax / DELTA * DELTA;
nmax += DELTA;
grow_arrays(nmax);
if (nlocal == nmax) {
nmax = nmax / DELTA * DELTA;
nmax += DELTA;
grow_arrays(nmax);
error->message(FLERR,
"in Fixtlsph_refconfigNeighGCG::unpack_exchange: local arrays too small for receiving partner information; growing arrays");
}
//printf("nlocal=%d, nmax=%d\n", nlocal, nmax);
if (comm->me == 0)
error->message(FLERR, "in Fixtlsph_refconfigNeighGCG::unpack_exchange: local arrays too small for receiving partner information; growing arrays");
}
int m = 0;
npartner[nlocal] = static_cast<int>(buf[m++]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<tagint>(buf[m++]);
wfd_list[nlocal][n] = static_cast<float>(buf[m++]);
wf_list[nlocal][n] = static_cast<float>(buf[m++]);
degradation_ij[nlocal][n] = static_cast<float>(buf[m++]);
energy_per_bond[nlocal][n] = static_cast<float>(buf[m++]);
}
return m;
int m = 0;
npartner[nlocal] = static_cast<int>(buf[m++]);
for (int n = 0; n < npartner[nlocal]; n++) {
partner[nlocal][n] = static_cast<tagint>(buf[m++]);
wfd_list[nlocal][n] = static_cast<float>(buf[m++]);
wf_list[nlocal][n] = static_cast<float>(buf[m++]);
degradation_ij[nlocal][n] = static_cast<float>(buf[m++]);
energy_per_bond[nlocal][n] = static_cast<float>(buf[m++]);
}
return m;
}
/* ----------------------------------------------------------------------
@ -450,18 +439,18 @@ int FixSMD_TLSPH_ReferenceConfiguration::unpack_exchange(int nlocal, double *buf
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_restart(int i, double *buf) {
int m = 0;
// pack buf[0] this way because other fixes unpack it
buf[m++] = 4 * npartner[i] + 2;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
int m = 0;
// pack buf[0] this way because other fixes unpack it
buf[m++] = 4 * npartner[i] + 2;
buf[m++] = npartner[i];
for (int n = 0; n < npartner[i]; n++) {
buf[m++] = partner[i][n];
buf[m++] = wfd_list[i][n];
buf[m++] = wf_list[i][n];
buf[m++] = degradation_ij[i][n];
buf[m++] = energy_per_bond[i][n];
}
return m;
}
/* ----------------------------------------------------------------------
@ -496,9 +485,9 @@ void FixSMD_TLSPH_ReferenceConfiguration::unpack_restart(int /*nlocal*/, int /*n
int FixSMD_TLSPH_ReferenceConfiguration::maxsize_restart() {
// maxtouch_all = max # of touching partners across all procs
int maxtouch_all;
MPI_Allreduce(&maxpartner, &maxtouch_all, 1, MPI_INT, MPI_MAX, world);
return 4 * maxtouch_all + 2;
int maxtouch_all;
MPI_Allreduce(&maxpartner, &maxtouch_all, 1, MPI_INT, MPI_MAX, world);
return 4 * maxtouch_all + 2;
}
/* ----------------------------------------------------------------------
@ -506,72 +495,71 @@ int FixSMD_TLSPH_ReferenceConfiguration::maxsize_restart() {
------------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::size_restart(int nlocal) {
return 4 * npartner[nlocal] + 2;
return 4 * npartner[nlocal] + 2;
}
/* ---------------------------------------------------------------------- */
int FixSMD_TLSPH_ReferenceConfiguration::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) {
int i, j, m;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
int i, j, m;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
//printf("FixSMD_TLSPH_ReferenceConfiguration:::pack_forward_comm\n");
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x0[j][0];
buf[m++] = x0[j][1];
buf[m++] = x0[j][2];
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x0[j][0];
buf[m++] = x0[j][1];
buf[m++] = x0[j][2];
buf[m++] = vfrac[j];
buf[m++] = radius[j];
buf[m++] = vfrac[j];
buf[m++] = radius[j];
buf[m++] = defgrad0[i][0];
buf[m++] = defgrad0[i][1];
buf[m++] = defgrad0[i][2];
buf[m++] = defgrad0[i][3];
buf[m++] = defgrad0[i][4];
buf[m++] = defgrad0[i][5];
buf[m++] = defgrad0[i][6];
buf[m++] = defgrad0[i][7];
buf[m++] = defgrad0[i][8];
buf[m++] = defgrad0[i][0];
buf[m++] = defgrad0[i][1];
buf[m++] = defgrad0[i][2];
buf[m++] = defgrad0[i][3];
buf[m++] = defgrad0[i][4];
buf[m++] = defgrad0[i][5];
buf[m++] = defgrad0[i][6];
buf[m++] = defgrad0[i][7];
buf[m++] = defgrad0[i][8];
}
return m;
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixSMD_TLSPH_ReferenceConfiguration::unpack_forward_comm(int n, int first, double *buf) {
int i, m, last;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
int i, m, last;
double *radius = atom->radius;
double *vfrac = atom->vfrac;
double **x0 = atom->x0;
double **defgrad0 = atom->smd_data_9;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x0[i][0] = buf[m++];
x0[i][1] = buf[m++];
x0[i][2] = buf[m++];
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x0[i][0] = buf[m++];
x0[i][1] = buf[m++];
x0[i][2] = buf[m++];
vfrac[i] = buf[m++];
radius[i] = buf[m++];
vfrac[i] = buf[m++];
radius[i] = buf[m++];
defgrad0[i][0] = buf[m++];
defgrad0[i][1] = buf[m++];
defgrad0[i][2] = buf[m++];
defgrad0[i][3] = buf[m++];
defgrad0[i][4] = buf[m++];
defgrad0[i][5] = buf[m++];
defgrad0[i][6] = buf[m++];
defgrad0[i][7] = buf[m++];
defgrad0[i][8] = buf[m++];
}
defgrad0[i][0] = buf[m++];
defgrad0[i][1] = buf[m++];
defgrad0[i][2] = buf[m++];
defgrad0[i][3] = buf[m++];
defgrad0[i][4] = buf[m++];
defgrad0[i][5] = buf[m++];
defgrad0[i][6] = buf[m++];
defgrad0[i][7] = buf[m++];
defgrad0[i][8] = buf[m++];
}
}
/* ----------------------------------------------------------------------
@ -582,26 +570,24 @@ void FixSMD_TLSPH_ReferenceConfiguration::unpack_forward_comm(int n, int first,
bool FixSMD_TLSPH_ReferenceConfiguration::crack_exclude(int i, int j) {
double **x = atom->x;
double l0 = domain->lattice->xlattice;
double **x = atom->x;
double l0 = domain->lattice->xlattice;
// line between pair of atoms i,j
double x1 = x[i][0] / l0;
double y1 = x[i][1] / l0;
// line between pair of atoms i,j
double x1 = x[i][0] / l0;
double y1 = x[i][1] / l0;
double x2 = x[j][0] / l0;
double y2 = x[j][1] / l0;
double x2 = x[j][0] / l0;
double y2 = x[j][1] / l0;
// hardcoded crack line
double x3 = -0.1 / l0;
double y3 = ((int) 1.0 / l0) + 0.5;
//printf("y3 = %f\n", y3);
double x4 = 0.1 / l0 - 1.0 + 0.1;
double y4 = y3;
// hardcoded crack line
double x3 = -0.1 / l0;
double y3 = ((int) 1.0 / l0) + 0.5;
double x4 = 0.1 / l0 - 1.0 + 0.1;
double y4 = y3;
bool retVal = DoLineSegmentsIntersect(x1, y1, x2, y2, x3, y3, x4, y4);
bool retVal = DoLineSegmentsIntersect(x1, y1, x2, y2, x3, y3, x4, y4);
return !retVal;
//return 1;
return !retVal;
}

3296
src/MACHDYN/pair_smd_tlsph.cpp
View File

File diff suppressed because it is too large Load Diff

20
src/MANIFOLD/fix_nve_manifold_rattle.cpp
View File

@ -537,13 +537,9 @@ void FixNVEManifoldRattle::rattle_manifold_x(double *x, double *v,
// gg = ptr_m->g(x);
}
if (iters >= max_iter && res > tolerance) {
char msg[2048];
sprintf(msg,"Failed to constrain atom " TAGINT_FORMAT
" (x = (%f, %f, %f)! res = %e, iters = %d\n",
tagi, x[0], x[1], x[2], res, iters);
error->one(FLERR,msg);
}
if (iters >= max_iter && res > tolerance)
error->one(FLERR, "Failed to constrain atom {} (x = ({}, {}, {})! res = {}, iters = {}\n",
tagi, x[0], x[1], x[2], res, iters);
// "sync" x and v:
v[0] = vt[0] - l*no_dt[0];
@ -631,13 +627,9 @@ void FixNVEManifoldRattle::rattle_manifold_v(double *v, double *f,
++iters;
} while ((res > tolerance) && (iters < max_iter));
if (iters >= max_iter && res >= tolerance) {
char msg[2048];
sprintf(msg,"Failed to constrain atom " TAGINT_FORMAT
" (x = (%f, %f, %f)! res = %e, iters = %d\n",
tagi, x[0], x[1], x[2], res, iters);
error->all(FLERR,msg);
}
if (iters >= max_iter && res >= tolerance)
error->all(FLERR,"Failed to constrain atom {} (x = ({}, {}, {})! res = {}, iters = {}\n",
tagi, x[0], x[1], x[2], res, iters);
stats.v_iters += iters;
}

65
src/MC/fix_bond_break.cpp
View File

@ -14,16 +14,17 @@
#include "fix_bond_break.h"
#include <cstring>
#include "update.h"
#include "respa.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "respa.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -157,15 +158,11 @@ void FixBondBreak::init()
else improperflag = 0;
if (force->improper) {
if (force->improper_match("class2") || force->improper_match("ring"))
error->all(FLERR,"Cannot yet use fix bond/break with this "
"improper style");
if (force->improper_match("^class2") || force->improper_match("^ring"))
error->all(FLERR,"Cannot yet use fix bond/break with this improper style");
}
lastcheck = -1;
// DEBUG
//print_bb();
}
/* ---------------------------------------------------------------------- */
@ -790,50 +787,6 @@ void FixBondBreak::unpack_reverse_comm(int n, int *list, double *buf)
}
}
/* ---------------------------------------------------------------------- */
void FixBondBreak::print_bb()
{
for (int i = 0; i < atom->nlocal; i++) {
printf("TAG " TAGINT_FORMAT ": %d nbonds: ",atom->tag[i],atom->num_bond[i]);
for (int j = 0; j < atom->num_bond[i]; j++) {
printf(" " TAGINT_FORMAT, atom->bond_atom[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d nangles: ",atom->tag[i],atom->num_angle[i]);
for (int j = 0; j < atom->num_angle[i]; j++) {
printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT ",",
atom->angle_atom1[i][j],atom->angle_atom2[i][j],
atom->angle_atom3[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d ndihedrals: ",atom->tag[i],atom->num_dihedral[i]);
for (int j = 0; j < atom->num_dihedral[i]; j++) {
printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT ",",
atom->dihedral_atom1[i][j],atom->dihedral_atom2[i][j],
atom->dihedral_atom3[i][j],atom->dihedral_atom4[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d %d %d nspecial: ",atom->tag[i],
atom->nspecial[i][0],atom->nspecial[i][1],atom->nspecial[i][2]);
for (int j = 0; j < atom->nspecial[i][2]; j++) {
printf(" " TAGINT_FORMAT, atom->special[i][j]);
}
printf("\n");
}
}
/* ---------------------------------------------------------------------- */
void FixBondBreak::print_copy(const char *str, tagint m,
int n1, int n2, int n3, int *v)
{
printf("%s " TAGINT_FORMAT ": %d %d %d nspecial: ",str,m,n1,n2,n3);
for (int j = 0; j < n3; j++) printf(" %d",v[j]);
printf("\n");
}
/* ---------------------------------------------------------------------- */
double FixBondBreak::compute_vector(int n)

5
src/MC/fix_bond_break.h
View File

@ -71,11 +71,6 @@ class FixBondBreak : public Fix {
void break_impropers(int, tagint, tagint);
void rebuild_special_one(int);
int dedup(int, int, tagint *);
// DEBUG
void print_bb();
void print_copy(const char *, tagint, int, int, int, int *);
};
} // namespace LAMMPS_NS

59
src/MC/fix_bond_create.cpp
View File

@ -621,9 +621,6 @@ void FixBondCreate::post_integrate()
// also add angles/dihedrals/impropers induced by created bonds
update_topology();
// DEBUG
//print_bb();
}
/* ----------------------------------------------------------------------
@ -685,11 +682,6 @@ void FixBondCreate::update_topology()
nimpropers = 0;
overflow = 0;
//printf("NCREATE %d: ",ncreate);
//for (i = 0; i < ncreate; i++)
// printf(" %d %d,",created[i][0],created[i][1]);
//printf("\n");
for (i = 0; i < nlocal; i++) {
influenced = 0;
slist = special[i];
@ -1439,54 +1431,3 @@ double FixBondCreate::memory_usage()
return bytes;
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::print_bb()
{
for (int i = 0; i < atom->nlocal; i++) {
printf("TAG " TAGINT_FORMAT ": %d nbonds: ",atom->tag[i],atom->num_bond[i]);
for (int j = 0; j < atom->num_bond[i]; j++) {
printf(" " TAGINT_FORMAT,atom->bond_atom[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d nangles: ",atom->tag[i],atom->num_angle[i]);
for (int j = 0; j < atom->num_angle[i]; j++) {
printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT ",",
atom->angle_atom1[i][j], atom->angle_atom2[i][j],
atom->angle_atom3[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d ndihedrals: ",atom->tag[i],atom->num_dihedral[i]);
for (int j = 0; j < atom->num_dihedral[i]; j++) {
printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT ",", atom->dihedral_atom1[i][j],
atom->dihedral_atom2[i][j],atom->dihedral_atom3[i][j],
atom->dihedral_atom4[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d nimpropers: ",atom->tag[i],atom->num_improper[i]);
for (int j = 0; j < atom->num_improper[i]; j++) {
printf(" " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT ",",atom->improper_atom1[i][j],
atom->improper_atom2[i][j],atom->improper_atom3[i][j],
atom->improper_atom4[i][j]);
}
printf("\n");
printf("TAG " TAGINT_FORMAT ": %d %d %d nspecial: ",atom->tag[i],
atom->nspecial[i][0],atom->nspecial[i][1],atom->nspecial[i][2]);
for (int j = 0; j < atom->nspecial[i][2]; j++) {
printf(" " TAGINT_FORMAT,atom->special[i][j]);
}
printf("\n");
}
}
/* ---------------------------------------------------------------------- */
void FixBondCreate::print_copy(const char *str, tagint m,
int n1, int n2, int n3, int *v)
{
printf("%s " TAGINT_FORMAT ": %d %d %d nspecial: ",str,m,n1,n2,n3);
for (int j = 0; j < n3; j++) printf(" %d",v[j]);
printf("\n");
}

5
src/MC/fix_bond_create.h
View File

@ -88,11 +88,6 @@ class FixBondCreate : public Fix {
int dedup(int, int, tagint *);
virtual int constrain(int, int, double, double) { return 1; }
// DEBUG
void print_bb();
void print_copy(const char *, tagint, int, int, int, int *);
};
} // namespace LAMMPS_NS