avoid division by zero in peri pair styles (358915d1) · Commits · 郑智淋 / lammps

src/PERI/pair_peri_eps.cpp

+18 −19

Original line number	Diff line number	Diff line
		@@ -46,7 +46,7 @@ PairPeriEPS::PairPeriEPS(LAMMPS *lmp) : Pair(lmp)

		ifix_peri = -1;

		nmax = 0;
		nmax = -1;
		s0_new = NULL;
		theta = NULL;

		@@ -209,7 +209,7 @@ void PairPeriEPS::compute(int eflag, int vflag)
		for (i = 0; i < nlocal; i++) maxpartner = MAX(maxpartner,npartner[i]);


		if (atom->nmax > nmax) {
		if (nlocal > nmax) {
		memory->destroy(s0_new);
		memory->destroy(theta);
		nmax = atom->nmax;
		@@ -220,13 +220,13 @@ void PairPeriEPS::compute(int eflag, int vflag)

		// ****** temp array to store Plastic extension ********* ///
		// create on heap to reduce stack use and to allow for faster zeroing
		double **deviatorPlasticExtTemp;
		double **deviatorPlasticExtTemp = NULL;
		if (nlocal*maxpartner > 0) {
		memory->create(deviatorPlasticExtTemp,nlocal,maxpartner,"pair:plastext");
		memset(&(deviatorPlasticExtTemp[0][0]),0,sizeof(double)nlocalmaxpartner);
		}
		// ****** temp array to store Plastic extension ********* ///



		// compute the dilatation on each particle
		compute_dilatation();

		@@ -280,12 +280,10 @@ void PairPeriEPS::compute(int eflag, int vflag)
		double fsurf = (tdnorm * tdnorm)/2 - pointwiseYieldvalue;
		bool elastic = true;

		double alphavalue = (15 * shearmodulus[itype][itype]) /wvolume[i];


		if (fsurf > 0) {
		elastic = false;
		deltalambda = ((tdnorm /sqrt(2.0 * pointwiseYieldvalue)) - 1.0) / alphavalue;
		deltalambda = ((tdnorm /sqrt(2.0 * pointwiseYieldvalue)) - 1.0) * wvolume[i]
		/ (15 * shearmodulus[itype][itype]);
		double templambda = lambdaValue[i];
		lambdaValue[i] = templambda + deltalambda;
		}
		@@ -350,8 +348,7 @@ void PairPeriEPS::compute(int eflag, int vflag)

		if (elastic) {
		rkNew = tdtrialValue;
		}
		else {
		} else {
		rkNew = (sqrt(2.0pointwiseYieldvalue) tdtrialValue) / tdnorm;
		deviatorPlasticExtTemp[i][jj] = edpNp1 + rkNew * deltalambda;
		}
		@@ -402,11 +399,13 @@ void PairPeriEPS::compute(int eflag, int vflag)

		memcpy(s0,s0_new,sizeof(double)*nlocal);

		if (nlocal*maxpartner > 0) {
		memcpy(&(deviatorPlasticextension[0][0]),
		&(deviatorPlasticExtTemp[0][0]),
		sizeof(double)nlocalmaxpartner);
		memory->destroy(deviatorPlasticExtTemp);
		}
		}

		/* ----------------------------------------------------------------------
		allocate all arrays

src/PERI/pair_peri_lps.cpp

+13 −9

Original line number	Diff line number	Diff line
		@@ -33,6 +33,7 @@
		#include "memory.h"
		#include "error.h"
		#include "update.h"
		#include "math_const.h"

		using namespace LAMMPS_NS;

		@@ -174,7 +175,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
		// of the bond-based theory used in PMB model

		double kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) /
		(3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype]);
		(MathConst::MY_PI * cutsq[itype][jtype] * cutsq[itype][jtype]);
		rk = (kshort * vfrac[j]) * (dr / cut[itype][jtype]);

		if (r > 0.0) fpair = -(rk/r);
		@@ -285,13 +286,14 @@ void PairPeriLPS::compute(int eflag, int vflag)

		omega_plus = influence_function(-1.0delx0,-1.0dely0,-1.0*delz0);
		omega_minus = influence_function(delx0,dely0,delz0);

		if ((wvolume[i] > 0.0) && (wvolume[j] > 0.0)) {
		rk = ( (3.0 * bulkmodulus[itype][itype]) -
		(5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale *
		( (omega_plus * theta[i] / wvolume[i]) +
		( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj];
		rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) *
		( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr;
		} else rk = 0.0;

		if (r > 0.0) fbond = -(rk/r);
		else fbond = 0.0;
		@@ -305,9 +307,11 @@ void PairPeriLPS::compute(int eflag, int vflag)
		double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);


		if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
		if (eflag && (wvolume[i] > 0.0))
		evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
		omega_plus(deviatoric_extension deviatoric_extension) *
		vfrac[j] * vfrac_scale;
		else evdwl = 0.0;
		if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,
		0.5fbondvfrac[i],delx,dely,delz);

src/USER-OMP/pair_peri_lps_omp.cpp

+11 −7

Original line number	Diff line number	Diff line
		@@ -310,12 +310,14 @@ void PairPeriLPSOMP::eval(int iifrom, int iito, ThrData * const thr)

		omega_plus = influence_function(-1.0delx0,-1.0dely0,-1.0*delz0);
		omega_minus = influence_function(delx0,dely0,delz0);
		if ((wvolume[i] > 0.0) && (wvolume[j] > 0.0)) {
		rk = ( (3.0 * bulkmodulus[itype][itype]) -
		(5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale *
		( (omega_plus * theta[i] / wvolume[i]) +
		( omega_minus * theta[j] / wvolume[j] ) ) * r0[i][jj];
		rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) *
		( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr;
		} else rk = 0.0;

		if (r > 0.0) fbond = -(rk/r);
		else fbond = 0.0;
		@@ -327,9 +329,11 @@ void PairPeriLPSOMP::eval(int iifrom, int iito, ThrData * const thr)
		// since I-J is double counted, set newton off & use 1/2 factor and I,I

		double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);
		if (EFLAG) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
		if (EFLAG && (wvolume[i] > 0.0))
		evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
		omega_plus(deviatoric_extension deviatoric_extension) *
		vfrac[j] * vfrac_scale;
		else evdwl = 0.0;
		if (EVFLAG) ev_tally_thr(this,i,i,nlocal,0,0.5*evdwl,0.0,
		0.5fbondvfrac[i],delx,dely,delz,thr);

Admin message