git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4906 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
192
src/min_fire.cpp
Normal file
192
src/min_fire.cpp
Normal file
@ -0,0 +1,192 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "min_fire.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "output.h"
|
||||
#include "timer.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
// same as in other min classes
|
||||
|
||||
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
|
||||
|
||||
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
|
||||
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
|
||||
|
||||
#define DELAYSTEP 5
|
||||
#define DT_GROW 1.1
|
||||
#define DT_SHRINK 0.5
|
||||
#define ALPHA0 0.1
|
||||
#define ALPHA_SHRINK 0.99
|
||||
#define TMAX 10.0
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
MinFire::MinFire(LAMMPS *lmp) : Min(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void MinFire::init_style()
|
||||
{
|
||||
dt = update->dt;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void MinFire::setup_style()
|
||||
{
|
||||
double **v = atom->v;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
v[i][0] = v[i][1] = v[i][2] = 0.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set current vector lengths and pointers
|
||||
called after atoms have migrated
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MinFire::reset_vectors()
|
||||
{
|
||||
// atomic dof
|
||||
|
||||
nvec = 3 * atom->nlocal;
|
||||
if (nvec) xvec = atom->x[0];
|
||||
if (nvec) fvec = atom->f[0];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int MinFire::iterate(int maxiter)
|
||||
{
|
||||
int ntimestep;
|
||||
double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall;
|
||||
double scale1,scale2;
|
||||
double dtvone,dtv,dtfm;
|
||||
|
||||
alpha_final = 0.0;
|
||||
double alpha = ALPHA0;
|
||||
double dtmax = TMAX * dt;
|
||||
int last_negative = update->ntimestep;
|
||||
|
||||
for (int iter = 0; iter < maxiter; iter++) {
|
||||
ntimestep = ++update->ntimestep;
|
||||
niter++;
|
||||
|
||||
// vdotfall = v dot f
|
||||
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
vdotf = 0.0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2];
|
||||
MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// if (v dot f) > 0:
|
||||
// v = (1-alpha) v + alpha |v| Fhat
|
||||
// |v| = length of v, Fhat = unit f
|
||||
// if more than DELAYSTEP since v dot f was negative:
|
||||
// increase timestep and decrease alpha
|
||||
|
||||
if (vdotfall > 0.0) {
|
||||
vdotv = 0.0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
vdotv += v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
|
||||
MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
fdotf = 0.0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
fdotf += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
|
||||
MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
scale1 = 1.0 - alpha;
|
||||
if (fdotfall == 0.0) scale2 = 0.0;
|
||||
else scale2 = alpha * sqrt(vdotvall/fdotfall);
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
|
||||
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
|
||||
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
|
||||
}
|
||||
|
||||
if (ntimestep - last_negative > DELAYSTEP) {
|
||||
dt = MIN(dt*DT_GROW,dtmax);
|
||||
alpha *= ALPHA_SHRINK;
|
||||
}
|
||||
|
||||
// else (v dot f) <= 0:
|
||||
// decrease timestep, reset alpha, set v = 0
|
||||
|
||||
} else {
|
||||
last_negative = ntimestep;
|
||||
dt *= DT_SHRINK;
|
||||
alpha = ALPHA0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
v[i][0] = v[i][1] = v[i][2] = 0.0;
|
||||
}
|
||||
|
||||
// limit timestep so no particle moves further than dmax
|
||||
|
||||
double *mass = atom->mass;
|
||||
int *type = atom->type;
|
||||
|
||||
dtvone = dt;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
vmax = MAX(fabs(v[i][0]),fabs(v[i][1]));
|
||||
vmax = MAX(vmax,fabs(v[i][2]));
|
||||
if (dtvone*vmax > dmax) dtvone = dmax/vmax;
|
||||
}
|
||||
MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world);
|
||||
|
||||
// Euler integration step
|
||||
|
||||
double **x = atom->x;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
dtfm = dtv / mass[type[i]];
|
||||
x[i][0] += dtv * v[i][0];
|
||||
x[i][1] += dtv * v[i][1];
|
||||
x[i][2] += dtv * v[i][2];
|
||||
v[i][0] += dtfm * f[i][0];
|
||||
v[i][1] += dtfm * f[i][1];
|
||||
v[i][2] += dtfm * f[i][2];
|
||||
}
|
||||
|
||||
eprevious = ecurrent;
|
||||
ecurrent = energy_force(0);
|
||||
neval++;
|
||||
|
||||
// force tolerance criterion
|
||||
|
||||
//fdotf = fnorm_sqr();
|
||||
//if (fdotf < update->ftol*update->ftol) return FTOL;
|
||||
|
||||
// output for thermo, dump, restart files
|
||||
|
||||
if (output->next == ntimestep) {
|
||||
timer->stamp();
|
||||
output->write(ntimestep);
|
||||
timer->stamp(TIME_OUTPUT);
|
||||
}
|
||||
}
|
||||
|
||||
return MAXITER;
|
||||
}
|
||||
Reference in New Issue
Block a user