scale diameter and charge incrementally without using fix store. Corrects bug of fix store array accessed out of bounds

This commit is contained in:
Jibril B. Coulibaly
2020-05-19 18:28:18 -05:00
parent e1362e9a61
commit 03419b398a
2 changed files with 12 additions and 17 deletions

View File

@ -238,6 +238,7 @@ int FixAdapt::setmask()
void FixAdapt::post_constructor() void FixAdapt::post_constructor()
{ {
if (!resetflag) return;
if (!diamflag && !chgflag) return; if (!diamflag && !chgflag) return;
// new id = fix-ID + FIX_STORE_ATTRIBUTE // new id = fix-ID + FIX_STORE_ATTRIBUTE
@ -268,9 +269,8 @@ void FixAdapt::post_constructor()
double *radius = atom->radius; double *radius = atom->radius;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
for (int i = 0; i < nall; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) vec[i] = radius[i]; if (mask[i] & groupbit) vec[i] = radius[i];
else vec[i] = 0.0; else vec[i] = 0.0;
} }
@ -292,9 +292,8 @@ void FixAdapt::post_constructor()
double *q = atom->q; double *q = atom->q;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
for (int i = 0; i < nall; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) vec[i] = q[i]; if (mask[i] & groupbit) vec[i] = q[i];
else vec[i] = 0.0; else vec[i] = 0.0;
} }
@ -436,10 +435,12 @@ void FixAdapt::init()
error->all(FLERR,"Fix adapt requires atom attribute mass"); error->all(FLERR,"Fix adapt requires atom attribute mass");
if (discflag && domain->dimension!=2) if (discflag && domain->dimension!=2)
error->all(FLERR,"Fix adapt requires 2d simulation"); error->all(FLERR,"Fix adapt requires 2d simulation");
if(scaleflag) diam_scale = 1.0;
} }
if (ad->aparam == CHARGE) { if (ad->aparam == CHARGE) {
if (!atom->q_flag) if (!atom->q_flag)
error->all(FLERR,"Fix adapt requires atom attribute charge"); error->all(FLERR,"Fix adapt requires atom attribute charge");
if(scaleflag) chg_scale = 1.0;
} }
} }
} }
@ -578,9 +579,6 @@ void FixAdapt::change_settings()
if (ad->aparam == DIAMETER) { if (ad->aparam == DIAMETER) {
double density; double density;
// Get initial diameter if `scale` keyword is used
double *vec = fix_diam->vstore;
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
@ -592,18 +590,15 @@ void FixAdapt::change_settings()
if (discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]); if (discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]);
else density = rmass[i] / (4.0*MY_PI/3.0 * else density = rmass[i] / (4.0*MY_PI/3.0 *
radius[i]*radius[i]*radius[i]); radius[i]*radius[i]*radius[i]);
if (scaleflag) radius[i] = value * vec[i]; if (scaleflag) radius[i] = value * radius[i] / diam_scale;
else radius[i] = 0.5*value; else radius[i] = 0.5*value;
if (discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density; if (discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density;
else rmass[i] = 4.0*MY_PI/3.0 * else rmass[i] = 4.0*MY_PI/3.0 *
radius[i]*radius[i]*radius[i] * density; radius[i]*radius[i]*radius[i] * density;
} }
if (scaleflag) diam_scale = value;
} else if (ad->aparam == CHARGE) { } else if (ad->aparam == CHARGE) {
// Get initial charge if `scale` keyword is used
double *vec = fix_chg->vstore;
double *q = atom->q; double *q = atom->q;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -611,9 +606,10 @@ void FixAdapt::change_settings()
for (i = 0; i < nall; i++) for (i = 0; i < nall; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (scaleflag) q[i] = value * vec[i]; if (scaleflag) q[i] = value * q[i] / chg_scale;
else q[i] = value; else q[i] = value;
} }
if (scaleflag) chg_scale = value;
} }
} }
} }
@ -681,9 +677,8 @@ void FixAdapt::restore_settings()
double *rmass = atom->rmass; double *rmass = atom->rmass;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
for (int i = 0; i < nall; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if(discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]); if(discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]);
else density = rmass[i] / (4.0*MY_PI/3.0 * else density = rmass[i] / (4.0*MY_PI/3.0 *
@ -699,9 +694,8 @@ void FixAdapt::restore_settings()
double *q = atom->q; double *q = atom->q;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
for (int i = 0; i < nall; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) q[i] = vec[i]; if (mask[i] & groupbit) q[i] = vec[i];
} }
} }

View File

@ -47,6 +47,7 @@ class FixAdapt : public Fix {
int nlevels_respa; int nlevels_respa;
char *id_fix_diam,*id_fix_chg; char *id_fix_diam,*id_fix_chg;
class FixStore *fix_diam,*fix_chg; class FixStore *fix_diam,*fix_chg;
double diam_scale,chg_scale;
int discflag; int discflag;
struct Adapt { struct Adapt {