Unverified Commit 4acdd161 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #1937 from jibril-b-coulibaly/fixadapt

fix adapt: implement scale keyword for diameter/charge and 2d compatibility
parents b74aabf0 832cb57d
Loading
Loading
Loading
Loading
+11 −14
Original line number Diff line number Diff line
@@ -14,7 +14,7 @@ Syntax
* adapt = style name of this fix command
* N = adapt simulation settings every this many timesteps
* one or more attribute/arg pairs may be appended
* attribute = *pair* or *kspace* or *atom*
* attribute = *pair* or *bond* or *kspace* or *atom*

  .. parsed-literal::

@@ -86,8 +86,8 @@ the end of a simulation. Even if *reset* is specified as *yes*\ , a
restart file written during a simulation will contain the modified
settings.

If the *scale* keyword is set to *no*\ , then the value the parameter is
set to will be whatever the variable generates.  If the *scale*
If the *scale* keyword is set to *no*\ , then the value of the altered
parameter will be whatever the variable generates.  If the *scale*
keyword is set to *yes*\ , then the value of the altered parameter will
be the initial value of that parameter multiplied by whatever the
variable generates.  I.e. the variable is now a "scale factor" applied
@@ -319,26 +319,23 @@ The *atom* keyword enables various atom properties to be changed. The
current list of atom parameters that can be varied by this fix:

* charge = charge on particle
* diameter = diameter of particle
* diameter, or, diameter/disc = diameter of particle

The *v_name* argument of the *atom* keyword is the name of an
:doc:`equal-style variable <variable>` which will be evaluated each time
this fix is invoked to set the parameter to a new value.  It should be
specified as v_name, where name is the variable name.  See the
this fix is invoked to set, or scale the parameter to a new value.
It should be specified as v_name, where name is the variable name.  See the
discussion above describing the formulas associated with equal-style
variables.  The new value is assigned to the corresponding attribute
for all atoms in the fix group.

.. note::

   The *atom* keyword works this way whether the *scale* keyword is
   set to *no* or *yes*\ .  I.e. the use of scale yes is not yet supported
   by the *atom* keyword.

If the atom parameter is *diameter* and per-atom density and per-atom
mass are defined for particles (e.g. :doc:`atom_style granular <atom_style>`), then the mass of each particle is also
changed when the diameter changes (density is assumed to stay
constant).
changed when the diameter changes. The mass is set from the particle volume
for 3d systems (density is assumed to stay constant). For 2d, the default is
for LAMMPS to model particles with a radius attribute as spheres.
However, if the atom parameter is *diameter/disc*, then the mass is
set from the particle area (the density is assumed to be in mass/distance^2 units).

For example, these commands would shrink the diameter of all granular
particles in the "center" group from 1.0 to 0.1 in a linear fashion
+40 −25
Original line number Diff line number Diff line
@@ -15,6 +15,7 @@
#include <cstring>
#include "atom.h"
#include "bond.h"
#include "domain.h"
#include "update.h"
#include "group.h"
#include "modify.h"
@@ -139,9 +140,11 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
    } else if (strcmp(arg[iarg],"atom") == 0) {
      if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command");
      adapt[nadapt].which = ATOM;
      if (strcmp(arg[iarg+1],"diameter") == 0) {
      if (strcmp(arg[iarg+1],"diameter") == 0 || strcmp(arg[iarg+1],"diameter/disc") == 0) {
        adapt[nadapt].aparam = DIAMETER;
        diamflag = 1;
        discflag = 0;
        if (strcmp(arg[iarg+1],"diameter/disc") == 0) discflag = 1;
      } else if (strcmp(arg[iarg+1],"charge") == 0) {
        adapt[nadapt].aparam = CHARGE;
        chgflag = 1;
@@ -235,7 +238,6 @@ int FixAdapt::setmask()

void FixAdapt::post_constructor()
{
  if (!resetflag) return;
  if (!diamflag && !chgflag) return;

  // new id = fix-ID + FIX_STORE_ATTRIBUTE
@@ -251,7 +253,7 @@ void FixAdapt::post_constructor()
  newarg[4] = (char *) "1";
  newarg[5] = (char *) "1";

  if (diamflag) {
  if (diamflag && atom->radius_flag) {
    int n = strlen(id) + strlen("_FIX_STORE_DIAM") + 1;
    id_fix_diam = new char[n];
    strcpy(id_fix_diam,id);
@@ -274,7 +276,7 @@ void FixAdapt::post_constructor()
    }
  }

  if (chgflag) {
  if (chgflag && atom->q_flag) {
    int n = strlen(id) + strlen("_FIX_STORE_CHG") + 1;
    id_fix_chg = new char[n];
    strcpy(id_fix_chg,id);
@@ -428,6 +430,10 @@ void FixAdapt::init()
      if (ad->aparam == DIAMETER) {
        if (!atom->radius_flag)
          error->all(FLERR,"Fix adapt requires atom attribute diameter");
        if (!atom->rmass_flag)
          error->all(FLERR,"Fix adapt requires atom attribute mass");
        if (discflag && domain->dimension!=2)
          error->all(FLERR,"Fix adapt requires 2d simulation");
      }
      if (ad->aparam == CHARGE) {
        if (!atom->q_flag)
@@ -568,38 +574,44 @@ void FixAdapt::change_settings()
      // also scale rmass to new value

      if (ad->aparam == DIAMETER) {
        int mflag = 0;
        if (atom->rmass_flag) mflag = 1;
        double density;

        // Get initial diameter if `scale` keyword is used

        double *vec = fix_diam->vstore;
        double *radius = atom->radius;
        double *rmass = atom->rmass;
        int *mask = atom->mask;
        int nlocal = atom->nlocal;
        int nall = nlocal + atom->nghost;

        if (mflag == 0) {
          for (i = 0; i < nall; i++)
            if (mask[i] & groupbit)
              radius[i] = 0.5*value;
        } else {
        for (i = 0; i < nall; i++)
          if (mask[i] & groupbit) {
              density = rmass[i] / (4.0*MY_PI/3.0 *
            if (discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]);
            else density = rmass[i] / (4.0*MY_PI/3.0 *
                                       radius[i]*radius[i]*radius[i]);
              radius[i] = 0.5*value;
              rmass[i] = 4.0*MY_PI/3.0 *
            if (scaleflag) radius[i] = value * vec[i];
            else radius[i] = 0.5*value;
            if (discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density;
            else rmass[i] = 4.0*MY_PI/3.0 *
                            radius[i]*radius[i]*radius[i] * density;
          }
        }

      } else if (ad->aparam == CHARGE) {

        // Get initial charge if `scale` keyword is used

        double *vec = fix_chg->vstore;
        double *q = atom->q;
        int *mask = atom->mask;
        int nlocal = atom->nlocal;
        int nall = nlocal + atom->nghost;

        for (i = 0; i < nall; i++)
          if (mask[i] & groupbit) q[i] = value;
          if (mask[i] & groupbit) {
            if (scaleflag) q[i] = value * vec[i];
            else q[i] = value;
          }
      }
    }
  }
@@ -607,7 +619,7 @@ void FixAdapt::change_settings()
  modify->addstep_compute(update->ntimestep + nevery);

  // re-initialize pair styles if any PAIR settings were changed
  // ditto for bond styles if any BOND setitings were changes
  // ditto for bond styles if any BOND settings were changed
  // this resets other coeffs that may depend on changed values,
  //   and also offset and tail corrections

@@ -670,10 +682,13 @@ void FixAdapt::restore_settings()

        for (int i = 0; i < nlocal; i++)
          if (mask[i] & groupbit) {
            density = rmass[i] / (4.0*MY_PI/3.0 *
            if(discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]);
            else density = rmass[i] / (4.0*MY_PI/3.0 *
                                       radius[i]*radius[i]*radius[i]);
            radius[i] = vec[i];
            rmass[i] = 4.0*MY_PI/3.0 * radius[i]*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 *
                            radius[i]*radius[i]*radius[i] * density;
          }
      }
      if (chgflag) {
+1 −0
Original line number Diff line number Diff line
@@ -47,6 +47,7 @@ class FixAdapt : public Fix {
  int nlevels_respa;
  char *id_fix_diam,*id_fix_chg;
  class FixStore *fix_diam,*fix_chg;
  int discflag;

  struct Adapt {
    int which,ivar;