Commit a20c38af authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

complete implementation of using an instance of fix STORE for persistent storage

parent cf91a938
Loading
Loading
Loading
Loading
+47 −26
Original line number Diff line number Diff line
@@ -71,6 +71,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
  // for CauchyStat

  usePK = 1;
  restartPK = 0;
  restart_global = 1;
  restart_stored = 0;

@@ -79,7 +80,6 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
  pcouple = NONE;
  drag = 0.0;
  allremap = 1;
  id_dilate = NULL;
  mtchain = mpchain = 3;
  nc_tchain = nc_pchain = 1;
  mtk_flag = 1;
@@ -94,8 +94,6 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :

  tcomputeflag = 0;
  pcomputeflag = 0;
  id_temp = NULL;
  id_press = NULL;

  // turn on tilt factor scaling, whenever applicable

@@ -348,14 +346,15 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
      } else error->all(FLERR,"Illegal fix nvt/npt/nph command");
      iarg += 2;
    } else if (strcmp(arg[iarg],"cauchystat") == 0) {
      usePK = 0;
      if (iarg+3 > narg) error->all(FLERR,
				    "Illegal fix npt cauchystat command:"
				    " wrong number of arguments");
      if (strcmp(arg[iarg+2],"yes") != 0 && strcmp(arg[iarg+2],"no") != 0)
	error->all(FLERR,"Illegal cauchystat continue value.  Must be 'yes' or 'no'");
        error->all(FLERR,"Illegal cauchystat continue value.  "
                   "Must be 'yes' or 'no'");
      alpha = force->numeric(FLERR,arg[iarg+1]);
      restartPK = !strcmp(arg[iarg+2],"yes");
      CauchyStat_init();
      iarg += 3;
    } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
      if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
@@ -594,6 +593,14 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :

/* ---------------------------------------------------------------------- */

void FixNH::post_constructor()
{
  if (!usePK) CauchyStat_init();
  else CauchyStat_cleanup();
}

/* ---------------------------------------------------------------------- */

FixNH::~FixNH()
{
  if (copymode) return;
@@ -2439,12 +2446,11 @@ void FixNH::CauchyStat_init()
  }

  if (!id_store) {
    int n = strlen(id) + 11;
    int n = strlen(id) + 14;
    id_store = new char[n];
    strcpy(id_store,id);
    strcat(id_store,"_FIX_STORE");
    strcat(id_store,"_FIX_NH_STORE");
  }

  restart_stored = modify->find_fix(id_store);

  if (restartPK==1 && restart_stored < 0)
@@ -2455,8 +2461,6 @@ void FixNH::CauchyStat_init()
    error->all(FLERR,"Illegal cauchystat command: "
               " Alpha cannot be zero or negative.");


  if (restartPK == 0) {
  if (restart_stored < 0) {
    char **newarg = new char *[6];
    newarg[0] = id_store;
@@ -2467,14 +2471,12 @@ void FixNH::CauchyStat_init()
    newarg[5] = (char *) "6";
    modify->add_fix(6,newarg);
    delete[] newarg;
    restart_stored = modify->find_fix(id_store);
  }
    int n = modify->find_fix(id_store);
    init_store = (FixStore *)modify->fix[n];
  }
  init_store = (FixStore *)modify->fix[restart_stored];

  initRUN = 0;
  initPK = 1;
  usePK = 0;

#define H0(row,col) (H0[(row-1)][(col-1)])
#define invH0(row,col) (invH0[(row-1)][(col-1)])
@@ -2498,7 +2500,26 @@ void FixNH::CauchyStat_init()
#undef invH0
}

/* ----------------------------------------------------------------------
   cleanup when not using Cauchystat
------------------------------------------------------------------------- */

void FixNH::CauchyStat_cleanup()
{
  if (id_store) {
    modify->delete_fix(id_store);
    delete[] id_store;
  } else {
    int n = strlen(id) + 14;
    id_store = new char[n];
    strcpy(id_store,id);
    strcat(id_store,"_FIX_NH_STORE");
    modify->delete_fix(id_store);
    delete[] id_store;
  }
  id_store = NULL;
  restart_stored = -1;
}

/* ----------------------------------------------------------------------
   Cauchystat
@@ -2589,6 +2610,7 @@ void FixNH::CauchyStat()
  //initialize:
  if (initPK==1) {
    if (restartPK==1) {
      double *setPKinit = init_store->astore[0];
      setPK(1,1)=setPKinit[0]; setPK(1,2)=setPKinit[1]; setPK(1,3)=setPKinit[2];
      setPK(2,1)=setPKinit[1]; setPK(2,2)=setPKinit[3]; setPK(2,3)=setPKinit[4];
      setPK(3,1)=setPKinit[2]; setPK(3,2)=setPKinit[4]; setPK(3,3)=setPKinit[5];
@@ -2623,14 +2645,13 @@ void FixNH::CauchyStat()
  p_hydro /= pdim;

  // save information for Cauchystat restart

  double *setPKinit = init_store->astore[0];
  setPKinit[0] = setcauchy(1,1);
  setPKinit[1] = setcauchy(1,2);
  setPKinit[2] = setcauchy(1,3);
  setPKinit[3] = setcauchy(2,2);
  setPKinit[4] = setcauchy(2,3);
  setPKinit[5] = setcauchy(3,3);
  restart_stored=1;

#undef H
#undef Hdot
+2 −2
Original line number Diff line number Diff line
@@ -21,6 +21,7 @@ namespace LAMMPS_NS {
class FixNH : public Fix {
 public:
  FixNH(class LAMMPS *, int, char **);
  virtual void post_constructor();
  virtual ~FixNH();
  int setmask();
  virtual void init();
@@ -151,8 +152,6 @@ class FixNH : public Fix {
  double setPK[3][3];   // current set values of the PK stress
                        // (this is modified until the cauchy
                        // stress converges)
  double setPKinit[6];  // initialization value of setPK for
                        // continuation runs
  double alpha;         // integration parameter for the cauchystat
  int initPK;           // 1 if setPK needs to be initialized either
                        // from cauchy or restart, else 0
@@ -164,6 +163,7 @@ class FixNH : public Fix {
                        // else 1 (pressure->vector available)

  void CauchyStat_init();
  void CauchyStat_cleanup();
  void CauchyStat();
  void CauchyStat_Step(double (&Fi)[3][3], double (&Fdot)[3][3],
                       double (&cauchy)[3][3], double (&setcauchy)[3][3],