tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

first version of fix ttm/grid

Browse Source
This commit is contained in:
Steve Plimpton
2021-06-09 09:03:50 -06:00
parent d74d7cfd5f
commit 456b81417d
7 changed files with 1062 additions and 225 deletions
Download Patch File Download Diff File Expand all files Collapse all files

372
src/MISC/fix_ttm.cpp
View File

@ -31,7 +31,6 @@
#include "memory.h"
#include "error.h"
#include "tokenizer.h"
using namespace LAMMPS_NS;
@ -108,45 +107,36 @@ FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) :
gfactor1 = new double[atom->ntypes+1];
gfactor2 = new double[atom->ntypes+1];
// allocate 3d grid variables
// check for allowed maxium number of total grid nodes
total_nnodes = (bigint)nxnodes * (bigint)nynodes * (bigint)nznodes;
if (total_nnodes > MAXSMALLINT)
error->all(FLERR,"Too many nodes in fix ttm");
memory->create(nsum,nxnodes,nynodes,nznodes,"ttm:nsum");
memory->create(nsum_all,nxnodes,nynodes,nznodes,"ttm:nsum_all");
memory->create(sum_vsq,nxnodes,nynodes,nznodes,"ttm:sum_vsq");
memory->create(sum_mass_vsq,nxnodes,nynodes,nznodes,"ttm:sum_mass_vsq");
memory->create(sum_vsq_all,nxnodes,nynodes,nznodes,"ttm:sum_vsq_all");
memory->create(sum_mass_vsq_all,nxnodes,nynodes,nznodes,
"ttm:sum_mass_vsq_all");
memory->create(T_electron_old,nxnodes,nynodes,nznodes,"ttm:T_electron_old");
memory->create(T_electron,nxnodes,nynodes,nznodes,"ttm:T_electron");
memory->create(net_energy_transfer,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer");
memory->create(net_energy_transfer_all,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer_all");
// allocate 3d grid variables
allocate_grid();
// allocate per-atom flangevin and zero it
// NOTE: is init to zero necessary?
flangevin = nullptr;
grow_arrays(atom->nmax);
// zero out the flangevin array
for (int i = 0; i < atom->nmax; i++) {
flangevin[i][0] = 0;
flangevin[i][1] = 0;
flangevin[i][2] = 0;
flangevin[i][0] = 0.0;
flangevin[i][1] = 0.0;
flangevin[i][2] = 0.0;
}
// set 2 callbacks
atom->add_callback(Atom::GROW);
atom->add_callback(Atom::RESTART);
// set initial electron temperatures from user input file
if (comm->me == 0) read_initial_electron_temperatures(arg[13]);
MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
read_initial_electron_temperatures(arg[13]);
}
/* ---------------------------------------------------------------------- */
@ -160,17 +150,9 @@ FixTTM::~FixTTM()
delete [] gfactor1;
delete [] gfactor2;
memory->destroy(nsum);
memory->destroy(nsum_all);
memory->destroy(sum_vsq);
memory->destroy(sum_mass_vsq);
memory->destroy(sum_vsq_all);
memory->destroy(sum_mass_vsq_all);
memory->destroy(T_electron_old);
memory->destroy(T_electron);
memory->destroy(flangevin);
memory->destroy(net_energy_transfer);
memory->destroy(net_energy_transfer_all);
deallocate_grid();
}
/* ---------------------------------------------------------------------- */
@ -203,13 +185,12 @@ void FixTTM::init()
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer_all[ixnode][iynode][iznode] = 0;
if (utils::strmatch(update->integrate_style,"^respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
// initialize grid quantities
init_grid();
}
/* ---------------------------------------------------------------------- */
@ -227,6 +208,25 @@ void FixTTM::setup(int vflag)
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_setup(int /*vflag*/)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// apply langevin forces that have been stored from previous run
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force(int /*vflag*/)
{
double **x = atom->x;
@ -279,32 +279,6 @@ void FixTTM::post_force(int /*vflag*/)
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_setup(int /*vflag*/)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// apply langevin forces that have been stored from previous run
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_respa_setup(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == nlevels_respa-1) post_force_setup(vflag);
@ -312,68 +286,9 @@ void FixTTM::post_force_respa_setup(int vflag, int ilevel, int /*iloop*/)
/* ---------------------------------------------------------------------- */
void FixTTM::reset_dt()
void FixTTM::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
for (int i = 1; i <= atom->ntypes; i++)
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
/* ----------------------------------------------------------------------
read in initial electron temperatures from a user-specified file
only called by proc 0
------------------------------------------------------------------------- */
void FixTTM::read_initial_electron_temperatures(const char *filename)
{
int ***T_initial_set;
memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm:T_initial_set");
memset(&T_initial_set[0][0][0],0,total_nnodes*sizeof(int));
std::string name = utils::get_potential_file_path(filename);
if (name.empty())
error->one(FLERR,"Cannot open input file: {}",
filename);
FILE *fpr = fopen(name.c_str(),"r");
// read initial electron temperature values from file
char line[MAXLINE];
int ixnode,iynode,iznode;
double T_tmp;
while (1) {
if (fgets(line,MAXLINE,fpr) == nullptr) break;
ValueTokenizer values(line);
if (values.has_next()) ixnode = values.next_int();
if (values.has_next()) iynode = values.next_int();
if (values.has_next()) iznode = values.next_int();
if (values.has_next()) T_tmp = values.next_double();
else error->one(FLERR,"Incorrect format in fix ttm input file");
// check correctness of input data
if ((ixnode < 0) || (ixnode >= nxnodes)
|| (iynode < 0) || (iynode >= nynodes)
|| (iznode < 0) || (iznode >= nznodes))
error->one(FLERR,"Fix ttm invalide node index in fix ttm input");
if (T_tmp < 0.0)
error->one(FLERR,"Fix ttm electron temperatures must be > 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
fclose(fpr);
// check completeness of input data
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
error->one(FLERR,"Initial temperatures not all set in fix ttm");
memory->destroy(T_initial_set);
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
@ -552,55 +467,20 @@ void FixTTM::end_of_step()
}
}
/* ----------------------------------------------------------------------
memory usage of 3d grid
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
double FixTTM::memory_usage()
void FixTTM::reset_dt()
{
double bytes = 0.0;
bytes += (double)5*total_nnodes * sizeof(int);
bytes += (double)14*total_nnodes * sizeof(double);
return bytes;
for (int i = 1; i <= atom->ntypes; i++)
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
/* ---------------------------------------------------------------------- */
void FixTTM::grow_arrays(int ngrow)
{
memory->grow(flangevin,ngrow,3,"TTM:flangevin");
}
/* ----------------------------------------------------------------------
return the energy of the electronic subsystem or the net_energy transfer
between the subsystems
------------------------------------------------------------------------- */
double FixTTM::compute_vector(int n)
{
double e_energy = 0.0;
double transfer_energy = 0.0;
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
e_energy +=
T_electron[ixnode][iynode][iznode]*electronic_specific_heat*
electronic_density*del_vol;
transfer_energy +=
net_energy_transfer_all[ixnode][iynode][iznode]*update->dt;
}
if (n == 0) return e_energy;
if (n == 1) return transfer_energy;
return 0.0;
memory->grow(flangevin,ngrow,3,"TTM:flangevin");
}
/* ----------------------------------------------------------------------
@ -685,6 +565,15 @@ void FixTTM::unpack_restart(int nlocal, int nth)
flangevin[nlocal][2] = extra[nlocal][m++];
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixTTM::size_restart(int /*nlocal*/)
{
return 4;
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
@ -695,10 +584,157 @@ int FixTTM::maxsize_restart()
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
return the energy of the electronic subsystem or the net_energy transfer
between the subsystems
------------------------------------------------------------------------- */
int FixTTM::size_restart(int /*nlocal*/)
double FixTTM::compute_vector(int n)
{
return 4;
double e_energy = 0.0;
double transfer_energy = 0.0;
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
e_energy +=
T_electron[ixnode][iynode][iznode]*electronic_specific_heat*
electronic_density*del_vol;
transfer_energy +=
net_energy_transfer_all[ixnode][iynode][iznode]*update->dt;
}
if (n == 0) return e_energy;
if (n == 1) return transfer_energy;
return 0.0;
}
/* ----------------------------------------------------------------------
memory usage for flangevin and 3d grid
------------------------------------------------------------------------- */
double FixTTM::memory_usage()
{
double bytes = 0.0;
bytes += (double)atom->nmax * 3 * sizeof(double);
bytes += (double)5*total_nnodes * sizeof(int);
bytes += (double)14*total_nnodes * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate 3d grid quantities
------------------------------------------------------------------------- */
void FixTTM::allocate_grid()
{
memory->create(nsum,nxnodes,nynodes,nznodes,"ttm:nsum");
memory->create(nsum_all,nxnodes,nynodes,nznodes,"ttm:nsum_all");
memory->create(sum_vsq,nxnodes,nynodes,nznodes,"ttm:sum_vsq");
memory->create(sum_mass_vsq,nxnodes,nynodes,nznodes,"ttm:sum_mass_vsq");
memory->create(sum_vsq_all,nxnodes,nynodes,nznodes,"ttm:sum_vsq_all");
memory->create(sum_mass_vsq_all,nxnodes,nynodes,nznodes,
"ttm:sum_mass_vsq_all");
memory->create(T_electron_old,nxnodes,nynodes,nznodes,"ttm:T_electron_old");
memory->create(T_electron,nxnodes,nynodes,nznodes,"ttm:T_electron");
memory->create(net_energy_transfer,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer");
memory->create(net_energy_transfer_all,nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer_all");
}
/* ----------------------------------------------------------------------
initialize 3d grid quantities
------------------------------------------------------------------------- */
void FixTTM::init_grid()
{
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer_all[ixnode][iynode][iznode] = 0;
}
/* ----------------------------------------------------------------------
deallocate 3d grid quantities
------------------------------------------------------------------------- */
void FixTTM::deallocate_grid()
{
memory->destroy(nsum);
memory->destroy(nsum_all);
memory->destroy(sum_vsq);
memory->destroy(sum_mass_vsq);
memory->destroy(sum_vsq_all);
memory->destroy(sum_mass_vsq_all);
memory->destroy(T_electron_old);
memory->destroy(T_electron);
memory->destroy(flangevin);
memory->destroy(net_energy_transfer);
memory->destroy(net_energy_transfer_all);
}
/* ----------------------------------------------------------------------
read in initial electron temperatures from a user-specified file
only called by proc 0
------------------------------------------------------------------------- */
void FixTTM::read_initial_electron_temperatures(const char *filename)
{
int ***T_initial_set;
memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm:T_initial_set");
memset(&T_initial_set[0][0][0],0,total_nnodes*sizeof(int));
std::string name = utils::get_potential_file_path(filename);
if (name.empty())
error->one(FLERR,"Cannot open input file: {}",
filename);
FILE *fpr = fopen(name.c_str(),"r");
// read initial electron temperature values from file
char line[MAXLINE];
int ixnode,iynode,iznode;
double T_tmp;
while (1) {
if (fgets(line,MAXLINE,fpr) == nullptr) break;
ValueTokenizer values(line);
if (values.has_next()) ixnode = values.next_int();
if (values.has_next()) iynode = values.next_int();
if (values.has_next()) iznode = values.next_int();
if (values.has_next()) T_tmp = values.next_double();
else error->one(FLERR,"Incorrect format in fix ttm input file");
// check correctness of input data
if ((ixnode < 0) || (ixnode >= nxnodes)
|| (iynode < 0) || (iynode >= nynodes)
|| (iznode < 0) || (iznode >= nznodes))
error->one(FLERR,"Fix ttm invalide node index in fix ttm input");
if (T_tmp < 0.0)
error->one(FLERR,"Fix ttm electron temperatures must be > 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
fclose(fpr);
// check completeness of input data
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
error->one(FLERR,"Initial temperatures not all set in fix ttm");
memory->destroy(T_initial_set);
MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
}

25
src/MISC/fix_ttm.h
View File

@ -27,27 +27,27 @@ namespace LAMMPS_NS {
class FixTTM : public Fix {
public:
FixTTM(class LAMMPS *, int, char **);
~FixTTM();
virtual ~FixTTM();
int setmask();
void init();
void setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void post_force_setup(int);
virtual void post_force(int);
void post_force_respa_setup(int, int, int);
void end_of_step();
void post_force_respa(int, int, int);
virtual void end_of_step();
void reset_dt();
void write_restart(FILE *);
void restart(char *);
void grow_arrays(int);
virtual void write_restart(FILE *);
virtual void restart(char *);
int pack_restart(int, double *);
void unpack_restart(int, int);
int size_restart(int);
int maxsize_restart();
double memory_usage();
void grow_arrays(int);
double compute_vector(int);
virtual double compute_vector(int);
virtual double memory_usage();
private:
protected:
int nfileevery;
int nlevels_respa;
int seed;
@ -65,7 +65,10 @@ class FixTTM : public Fix {
double electronic_thermal_conductivity;
double gamma_p, gamma_s, v_0, v_0_sq;
void read_initial_electron_temperatures(const char *);
virtual void allocate_grid();
virtual void init_grid();
virtual void deallocate_grid();
virtual void read_initial_electron_temperatures(const char *);
};
} // namespace LAMMPS_NS

609
src/MISC/fix_ttm_grid.cpp Normal file
View File

@ -0,0 +1,609 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Paul Crozier (SNL)
Carolyn Phillips (University of Michigan)
------------------------------------------------------------------------- */
#include "fix_ttm_grid.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "update.h"
#include "domain.h"
#include "comm.h"
#include "gridcomm.h"
#include "neighbor.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "tokenizer.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define MAXLINE 1024
#define OFFSET 16384
enum{T_ELECTRON,NET_ENERGY_TRANSFER,DSUM,SUM_VSQ,SUM_MASS_VSQ};
/* ---------------------------------------------------------------------- */
FixTTMGrid::FixTTMGrid(LAMMPS *lmp, int narg, char **arg) :
FixTTM(lmp, narg, arg) {}
/* ---------------------------------------------------------------------- */
void FixTTMGrid::post_force(int /*vflag*/)
{
int ix,iy,iz;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double gamma1,gamma2;
double *boxlo = domain->boxlo;
// apply damping and thermostat to all atoms in fix group
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ix = static_cast<int> ((x[i][0]-boxlo[0])*delxinv + shift) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*delyinv + shift) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*delzinv + shift) - OFFSET;
if (T_electron[ix][iy][iz] < 0)
error->all(FLERR,"Electronic temperature dropped below zero");
// check that ix,iy,iz is within my ghost cell range
if (ix < nxlo_out || ix > nxhi_out || iy < nylo_out || iy > nyhi_out ||
iz < nzlo_out || iz > nzhi_out) flag = 1;
double tsqrt = sqrt(T_electron[ix][iy][iz]);
gamma1 = gfactor1[type[i]];
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p;
gamma2 = gfactor2[type[i]] * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
if (flag) error->one(FLERR,"Out of range fix ttm/grid atoms");
}
/* ---------------------------------------------------------------------- */
void FixTTMGrid::end_of_step()
{
int ix,iy,iz;
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *boxlo = domain->boxlo;
memset(&net_energy_transfer[nzlo_out][nylo_out][nxlo_out],0,
ngridout*sizeof(double));
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
ix = static_cast<int> ((x[i][0]-boxlo[0])*delxinv+shift) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*delyinv+shift) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*delzinv+shift) - OFFSET;
net_energy_transfer[ix][iy][iz] +=
(flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2]);
}
gc->reverse_comm(1,this,1,sizeof(double),NET_ENERGY_TRANSFER,
gc_buf1,gc_buf2,MPI_DOUBLE);
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
// num_inner_timesteps = # of inner steps (thermal solves)
// required this MD step to maintain a stable explicit solve
int num_inner_timesteps = 1;
double inner_dt = update->dt;
double stability_criterion = 1.0 -
2.0*inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
if (stability_criterion < 0.0) {
inner_dt = 0.5*(electronic_specific_heat*electronic_density) /
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
num_inner_timesteps = static_cast<int>(update->dt/inner_dt) + 1;
inner_dt = update->dt/double(num_inner_timesteps);
if (num_inner_timesteps > 1000000)
error->warning(FLERR,"Too many inner timesteps in fix ttm/grid");
}
for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
ith_inner_timestep++) {
memcpy(&T_electron_old[nzlo_out][nylo_out][nxlo_out],
&T_electron[nzlo_out][nylo_out][nxlo_out],ngridout*sizeof(double));
// compute new electron T profile
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
T_electron[ix][iy][iz] =
T_electron_old[ix][iy][iz] +
inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity *
((T_electron_old[ix-1][iy][iz] + T_electron_old[ix+1][iy][iz] -
2*T_electron_old[ix][iy][iz])/dx/dx +
(T_electron_old[ix][iy-1][iz] + T_electron_old[ix][iy+1][iz] -
2*T_electron_old[ix][iy][iz])/dy/dy +
(T_electron_old[ix][iy][iz-1] + T_electron_old[ix][iy][iz+1] -
2*T_electron_old[ix][iy][iz])/dz/dz) -
net_energy_transfer_all[ix][iy][iz]/del_vol);
}
}
// comm new T_electron to ghost grid points
gc->forward_comm(1,this,1,sizeof(double),T_ELECTRON,
gc_buf1,gc_buf2,MPI_DOUBLE);
// output nodal temperatures for current timestep
if ((nfileevery) && !(update->ntimestep % nfileevery)) {
// compute atomic Ta for each grid point
memset(&dsum[nzlo_out][nylo_out][nxlo_out],0,
ngridout*sizeof(double));
memset(&sum_vsq[nzlo_out][nylo_out][nxlo_out],0,
ngridout*sizeof(double));
memset(&sum_mass_vsq[nzlo_out][nylo_out][nxlo_out],0,
ngridout*sizeof(double));
double massone;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
ix = static_cast<int> ((x[i][0]-boxlo[0])*delxinv+shift) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*delyinv+shift) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*delzinv+shift) - OFFSET;
dsum[ix][iy][iz] += 1.0;
sum_vsq[ix][iy][iz] += vsq;
sum_mass_vsq[ix][iy][iz] += massone*vsq;
}
// can GridComm work on a larger struct of several values?
gc->reverse_comm(1,this,1,sizeof(double),DSUM,
gc_buf1,gc_buf2,MPI_DOUBLE);
gc->reverse_comm(1,this,1,sizeof(double),SUM_VSQ,
gc_buf1,gc_buf2,MPI_DOUBLE);
gc->reverse_comm(1,this,1,sizeof(double),SUM_MASS_VSQ,
gc_buf1,gc_buf2,MPI_DOUBLE);
// NOTE: should this be a write function, work with parallel pings
if (comm->me == 0) {
fmt::print(fp,"{}",update->ntimestep);
double T_a;
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
T_a = 0;
if (dsum[ix][iy][iz] > 0)
T_a = sum_mass_vsq[ix][iy][iz] /
(3.0*force->boltz*dsum[ix][iy][iz]/force->mvv2e);
fmt::print(fp," {}",T_a);
}
fputs("\t",fp);
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
fmt::print(fp," {}",T_electron[ix][iy][iz]);
fputs("\n",fp);
}
}
}
/* ----------------------------------------------------------------------
pack own values to buf to send to another proc
------------------------------------------------------------------------- */
void FixTTMGrid::pack_forward_grid(int flag, void *vbuf, int nlist, int *list)
{
double *buf = (double *) vbuf;
double *src;
if (flag == T_ELECTRON)
src = &T_electron[nzlo_out][nylo_out][nxlo_out];
for (int i = 0; i < nlist; i++)
buf[i] = src[list[i]];
}
/* ----------------------------------------------------------------------
unpack another proc's own values from buf and set own ghost values
------------------------------------------------------------------------- */
void FixTTMGrid::unpack_forward_grid(int flag, void *vbuf, int nlist, int *list)
{
double *buf = (double *) vbuf;
double *dest;
if (flag == T_ELECTRON)
dest = &T_electron[nzlo_out][nylo_out][nxlo_out];
for (int i = 0; i < nlist; i++)
dest[list[i]] = buf[i];
}
/* ----------------------------------------------------------------------
pack ghost values into buf to send to another proc
------------------------------------------------------------------------- */
void FixTTMGrid::pack_reverse_grid(int flag, void *vbuf, int nlist, int *list)
{
double *buf = (double *) vbuf;
double *src;
if (flag == NET_ENERGY_TRANSFER)
src = &net_energy_transfer[nzlo_out][nylo_out][nxlo_out];
else if (flag == DSUM)
src = &dsum[nzlo_out][nylo_out][nxlo_out];
else if (flag == SUM_VSQ)
src = &sum_vsq[nzlo_out][nylo_out][nxlo_out];
else if (flag == SUM_MASS_VSQ)
src = &sum_mass_vsq[nzlo_out][nylo_out][nxlo_out];
for (int i = 0; i < nlist; i++)
buf[i] = src[list[i]];
}
/* ----------------------------------------------------------------------
unpack another proc's ghost values from buf and add to own values
------------------------------------------------------------------------- */
void FixTTMGrid::unpack_reverse_grid(int flag, void *vbuf, int nlist, int *list)
{
double *buf = (double *) vbuf;
double *dest;
if (flag == NET_ENERGY_TRANSFER)
dest = &net_energy_transfer[nzlo_out][nylo_out][nxlo_out];
else if (flag == DSUM)
dest = &dsum[nzlo_out][nylo_out][nxlo_out];
else if (flag == SUM_VSQ)
dest = &sum_vsq[nzlo_out][nylo_out][nxlo_out];
else if (flag == SUM_MASS_VSQ)
dest = &sum_mass_vsq[nzlo_out][nylo_out][nxlo_out];
for (int i = 0; i < nlist; i++)
dest[list[i]] += buf[i];
}
/* ----------------------------------------------------------------------
return the energy of the electronic subsystem or the net_energy transfer
between the subsystems
------------------------------------------------------------------------- */
double FixTTMGrid::compute_vector(int n)
{
int ix,iy,iz;
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
double e_energy_me = 0.0;
double transfer_energy_me = 0.0;
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
e_energy_me +=
T_electron[ix][iy][iz]*electronic_specific_heat*
electronic_density*del_vol;
transfer_energy_me +=
net_energy_transfer_all[ix][iy][iz]*update->dt;
}
// NOTE: only do allreduce once ?
double e_energy,transfer_energy;
MPI_Allreduce(&e_energy_me,&e_energy,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&transfer_energy_me,&transfer_energy,1,MPI_DOUBLE,MPI_SUM,world);
if (n == 0) return e_energy;
if (n == 1) return transfer_energy;
return 0.0;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixTTMGrid::write_restart(FILE *fp)
{
int ix,iy,iz;
double *rlist;
memory->create(rlist,nxnodes*nynodes*nznodes+1,"TTM:rlist");
int n = 0;
rlist[n++] = seed;
// NOTE: this is a parallel grid now, not a serial grid
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
rlist[n++] = T_electron[ix][iy][iz];
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
fwrite(rlist,sizeof(double),n,fp);
}
memory->destroy(rlist);
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTTMGrid::restart(char *buf)
{
int ix,iy,iz;
int n = 0;
double *rlist = (double *) buf;
// the seed must be changed from the initial seed
// NOTE: 0.5 is whacky, could go to zero
// NOTE: maybe GridComm should provide a method to pack a grid with bounds
seed = static_cast<int> (0.5*rlist[n++]);
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
T_electron[ix][iy][iz] = rlist[n++];
delete random;
random = new RanMars(lmp,seed+comm->me);
}
/* ----------------------------------------------------------------------
memory usage for flangevin and 3d grid
------------------------------------------------------------------------- */
double FixTTMGrid::memory_usage()
{
double bytes = 0.0;
bytes += (double)3*atom->nmax * sizeof(double);
bytes += (double)6*ngridout * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate 3d grid quantities
------------------------------------------------------------------------- */
void FixTTMGrid::allocate_grid()
{
// global indices of grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global grid that I own without ghost cells
// both non-tiled and tiled proc layouts use 0-1 fractional subdomain info
if (comm->layout != Comm::LAYOUT_TILED) {
nxlo_in = static_cast<int> (comm->xsplit[comm->myloc[0]] * nxnodes);
nxhi_in = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nxnodes) - 1;
nylo_in = static_cast<int> (comm->ysplit[comm->myloc[1]] * nynodes);
nyhi_in = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * nynodes) - 1;
nzlo_in = static_cast<int> (comm->zsplit[comm->myloc[2]] * nznodes);
nzhi_in = static_cast<int> (comm->zsplit[comm->myloc[2]+1] * nznodes) - 1;
} else {
nxlo_in = static_cast<int> (comm->mysplit[0][0] * nxnodes);
nxhi_in = static_cast<int> (comm->mysplit[0][1] * nxnodes) - 1;
nylo_in = static_cast<int> (comm->mysplit[1][0] * nynodes);
nyhi_in = static_cast<int> (comm->mysplit[1][1] * nynodes) - 1;
nzlo_in = static_cast<int> (comm->mysplit[1][0] * nznodes);
nzhi_in = static_cast<int> (comm->mysplit[1][1] * nznodes) - 1;
}
// nlo,nhi = min/max index of global grid pt my owned atoms can be mapped to
// finite difference stencil requires extra grid pt around my owned grid pts
// max of these 2 quantities is the ghost cells needed in each di
// nlo_out,nhi_out = nlo_in,nhi_in + ghost cells
// NOTE: need to wait until init() when neighskin is defined?
double *boxlo = domain->boxlo;
double *sublo = domain->sublo;
double *subhi = domain->subhi;
shift = OFFSET + 0.0; // change this to 0.5 for nearest grid pt
delxinv = nxnodes/domain->xprd;
delyinv = nxnodes/domain->yprd;
delzinv = nxnodes/domain->zprd;
int nlo,nhi;
double cuthalf = 0.5*neighbor->skin;
nlo = static_cast<int> ((sublo[0]-cuthalf-boxlo[0])*delxinv + shift) - OFFSET;
nhi = static_cast<int> ((subhi[0]+cuthalf-boxlo[0])*delxinv + shift) - OFFSET;
nxlo_out = MIN(nlo,nxlo_in-1);
nxhi_out = MAX(nhi,nxhi_in+1);
nlo = static_cast<int> ((sublo[1]-cuthalf-boxlo[1])*delyinv + shift) - OFFSET;
nhi = static_cast<int> ((subhi[1]+cuthalf-boxlo[1])*delyinv + shift) - OFFSET;
nylo_out = MIN(nlo,nylo_in-1);
nyhi_out = MAX(nhi,nyhi_in+1);
nlo = static_cast<int> ((sublo[2]-cuthalf-boxlo[2])*delzinv + shift) - OFFSET;
nhi = static_cast<int> ((subhi[2]+cuthalf-boxlo[2])*delzinv + shift) - OFFSET;
nzlo_out = MIN(nlo,nzlo_in-1);
nzhi_out = MAX(nhi,nzhi_in+1);
ngridout = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
(nzhi_out-nzlo_out+1);
gc = new GridComm(lmp,world,nxnodes,nynodes,nznodes,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc->setup(ngc_buf1,ngc_buf2);
memory->create(gc_buf1,ngc_buf1,"ttm/grid:gc_buf1");
memory->create(gc_buf2,ngc_buf2,"ttm/grid:gc_buf2");
memory->create3d_offset(nsum,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"ttm:nsum");
memory->create3d_offset(sum_vsq,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"ttm:sum_vsq");
memory->create3d_offset(sum_mass_vsq,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"ttm:sum_mass_vsq");
memory->create3d_offset(T_electron_old,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"ttm:T_electron_old");
memory->create3d_offset(T_electron,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"ttm:T_electron");
memory->create3d_offset(net_energy_transfer,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"ttm:net_energy_transfer");
}
/* ----------------------------------------------------------------------
initialize 3d grid quantities
------------------------------------------------------------------------- */
void FixTTMGrid::init_grid()
{
memset(&net_energy_transfer[nzlo_out][nylo_out][nxlo_out],0,
ngridout*sizeof(double));
}
/* ----------------------------------------------------------------------
deallocate 3d grid quantities
------------------------------------------------------------------------- */
void FixTTMGrid::deallocate_grid()
{
memory->destroy3d_offset(nsum,nzlo_out,nylo_out,nxlo_out);
memory->destroy3d_offset(sum_vsq,nzlo_out,nylo_out,nxlo_out);
memory->destroy3d_offset(sum_mass_vsq,nzlo_out,nylo_out,nxlo_out);
memory->destroy3d_offset(T_electron_old,nzlo_out,nylo_out,nxlo_out);
memory->destroy3d_offset(T_electron,nzlo_out,nylo_out,nxlo_out);
memory->destroy3d_offset(net_energy_transfer,nzlo_out,nylo_out,nxlo_out);
}
/* ----------------------------------------------------------------------
read in initial electron temperatures from a user-specified file
only called by proc 0
------------------------------------------------------------------------- */
void FixTTMGrid::read_initial_electron_temperatures(const char *filename)
{
int ***T_initial_set;
memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm:T_initial_set");
memset(&T_initial_set[0][0][0],0,total_nnodes*sizeof(int));
std::string name = utils::get_potential_file_path(filename);
if (name.empty())
error->one(FLERR,"Cannot open input file: {}",filename);
FILE *fpr = fopen(name.c_str(),"r");
// read initial electron temperature values from file
char line[MAXLINE];
int ixnode,iynode,iznode;
double T_tmp;
while (1) {
if (fgets(line,MAXLINE,fpr) == nullptr) break;
ValueTokenizer values(line);
if (values.has_next()) ixnode = values.next_int();
if (values.has_next()) iynode = values.next_int();
if (values.has_next()) iznode = values.next_int();
if (values.has_next()) T_tmp = values.next_double();
else error->one(FLERR,"Incorrect format in fix ttm/grid input file");
// check correctness of input data
if ((ixnode < 0) || (ixnode >= nxnodes)
|| (iynode < 0) || (iynode >= nynodes)
|| (iznode < 0) || (iznode >= nznodes))
error->one(FLERR,"Fix ttm/grid invalid node index in input");
if (T_tmp < 0.0)
error->one(FLERR,"Fix ttm/grid electron temperatures must be > 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
fclose(fpr);
// check completeness of input data
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
error->one(FLERR,"Initial temperatures not all set in fix ttm/grid");
memory->destroy(T_initial_set);
}

148
src/MISC/fix_ttm_grid.h Normal file
View File

@ -0,0 +1,148 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(ttm/grid,FixTTMGrid);
// clang-format on
#else
#ifndef LMP_FIX_TTM_GRID_H
#define LMP_FIX_TTM_GRID_H
#include "fix_ttm.h"
namespace LAMMPS_NS {
class FixTTMGrid : public FixTTM {
public:
FixTTMGrid(class LAMMPS *, int, char **);
~FixTTMGrid() {}
void post_force(int);
void end_of_step();
double compute_vector(int);
void write_restart(FILE *);
void restart(char *);
double memory_usage();
// grid communication
void pack_forward_grid(int, void *, int, int *);
void unpack_forward_grid(int, void *, int, int *);
void pack_reverse_grid(int, void *, int, int *);
void unpack_reverse_grid(int, void *, int, int *);
private:
int ngridout;
int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in;
int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out;
double shift;
double delxinv,delyinv,delzinv;
double ***dsum;
class GridComm *gc;
int ngc_buf1,ngc_buf2;
double *gc_buf1,*gc_buf2;
void allocate_grid();
void init_grid();
void deallocate_grid();
void read_initial_electron_temperatures(const char *);
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot open file %s
The specified file cannot be opened. Check that the path and name are
correct. If the file is a compressed file, also check that the gzip
executable can be found and run.
E: Cannot open fix ttm file %s
The output file for the fix ttm command cannot be opened. Check that
the path and name are correct.
E: Invalid random number seed in fix ttm command
Random number seed must be > 0.
E: Fix ttm electronic_specific_heat must be > 0.0
Self-explanatory.
E: Fix ttm electronic_density must be > 0.0
Self-explanatory.
E: Fix ttm electronic_thermal_conductivity must be >= 0.0
Self-explanatory.
E: Fix ttm gamma_p must be > 0.0
Self-explanatory.
E: Fix ttm gamma_s must be >= 0.0
Self-explanatory.
E: Fix ttm v_0 must be >= 0.0
Self-explanatory.
E: Fix ttm number of nodes must be > 0
Self-explanatory.
E: Cannot use fix ttm with 2d simulation
This is a current restriction of this fix due to the grid it creates.
E: Cannot use non-periodic boundares with fix ttm
This fix requires a fully periodic simulation box.
E: Cannot use fix ttm with triclinic box
This is a current restriction of this fix due to the grid it creates.
E: Electronic temperature dropped below zero
Something has gone wrong with the fix ttm electron temperature model.
E: Fix ttm electron temperatures must be > 0.0
Self-explanatory.
E: Initial temperatures not all set in fix ttm
Self-explanatory.
W: Too many inner timesteps in fix ttm
Self-explanatory.
*/

5
src/fix.h
View File

@ -205,6 +205,11 @@ class Fix : protected Pointers {
virtual int pack_reverse_comm(int, int, double *) { return 0; }
virtual void unpack_reverse_comm(int, int *, double *) {}
virtual void pack_forward_grid(int, void *, int, int *) {};
virtual void unpack_forward_grid(int, void *, int, int *) {};
virtual void pack_reverse_grid(int, void *, int, int *) {};
virtual void unpack_reverse_grid(int, void *, int, int *) {};
virtual double compute_scalar() { return 0.0; }
virtual double compute_vector(int) { return 0.0; }
virtual double compute_array(int, int) { return 0.0; }

112
src/KSPACE/gridcomm.cpp → src/gridcomm.cpp
View File

@ -18,11 +18,13 @@
#include "error.h"
#include "irregular.h"
#include "kspace.h"
#include "fix.h"
#include "memory.h"
using namespace LAMMPS_NS;
enum{REGULAR,TILED};
enum{KSPACE,FIX};
#define DELTA 16
@ -164,10 +166,14 @@ GridComm::~GridComm()
void GridComm::initialize(MPI_Comm gcomm,
int gnx, int gny, int gnz,
int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi,
int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi,
int fxlo, int fxhi, int fylo, int fyhi, int fzlo, int fzhi,
int pxlo, int pxhi, int pylo, int pyhi, int pzlo, int pzhi)
int ixlo, int ixhi, int iylo, int iyhi,
int izlo, int izhi,
int oxlo, int oxhi, int oylo, int oyhi,
int ozlo, int ozhi,
int fxlo, int fxhi, int fylo, int fyhi,
int fzlo, int fzhi,
int pxlo, int pxhi, int pylo, int pyhi,
int pzlo, int pzhi)
{
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
@ -926,31 +932,43 @@ int GridComm::ghost_adjacent_tiled()
forward comm of my owned cells to other's ghost cells
------------------------------------------------------------------------- */
void GridComm::forward_comm_kspace(KSpace *kspace, int nper, int nbyte, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
void GridComm::forward_comm(int caller, void *ptr, int nper, int nbyte, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
{
if (layout == REGULAR)
forward_comm_kspace_regular(kspace,nper,nbyte,which,buf1,buf2,datatype);
else
forward_comm_kspace_tiled(kspace,nper,nbyte,which,buf1,buf2,datatype);
if (layout == REGULAR) {
if (caller == KSPACE)
forward_comm_regular<KSpace>((KSpace *)ptr,nper,nbyte,which,
buf1,buf2,datatype);
else if (caller == FIX)
forward_comm_regular<Fix>((Fix *)ptr,nper,nbyte,which,
buf1,buf2,datatype);
} else {
if (caller == KSPACE)
forward_comm_tiled<KSpace>((KSpace *)ptr,nper,nbyte,which,
buf1,buf2,datatype);
else if (caller == FIX)
forward_comm_tiled<Fix>((Fix *)ptr,nper,nbyte,
which,buf1,buf2,datatype);
}
}
/* ----------------------------------------------------------------------
forward comm on regular grid of procs via list of swaps with 6 neighbor procs
------------------------------------------------------------------------- */
template < class T >
void GridComm::
forward_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
forward_comm_regular(T *ptr, int nper, int /*nbyte*/, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
{
int m;
MPI_Request request;
for (m = 0; m < nswap; m++) {
if (swap[m].sendproc == me)
kspace->pack_forward_grid(which,buf2,swap[m].npack,swap[m].packlist);
ptr->pack_forward_grid(which,buf2,swap[m].npack,swap[m].packlist);
else
kspace->pack_forward_grid(which,buf1,swap[m].npack,swap[m].packlist);
ptr->pack_forward_grid(which,buf1,swap[m].npack,swap[m].packlist);
if (swap[m].sendproc != me) {
if (swap[m].nunpack) MPI_Irecv(buf2,nper*swap[m].nunpack,datatype,
@ -960,7 +978,7 @@ forward_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which,
if (swap[m].nunpack) MPI_Wait(&request,MPI_STATUS_IGNORE);
}
kspace->unpack_forward_grid(which,buf2,swap[m].nunpack,swap[m].unpacklist);
ptr->unpack_forward_grid(which,buf2,swap[m].nunpack,swap[m].unpacklist);
}
}
@ -968,9 +986,10 @@ forward_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which,
forward comm on tiled grid decomp via Send/Recv lists of each neighbor proc
------------------------------------------------------------------------- */
template < class T >
void GridComm::
forward_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which,
void *buf1, void *vbuf2, MPI_Datatype datatype)
forward_comm_tiled(T *ptr, int nper, int nbyte, int which,
void *buf1, void *vbuf2, MPI_Datatype datatype)
{
int i,m,offset;
@ -987,15 +1006,15 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which,
// perform all sends to other procs
for (m = 0; m < nsend; m++) {
kspace->pack_forward_grid(which,buf1,send[m].npack,send[m].packlist);
ptr->pack_forward_grid(which,buf1,send[m].npack,send[m].packlist);
MPI_Send(buf1,nper*send[m].npack,datatype,send[m].proc,0,gridcomm);
}
// perform all copies to self
for (m = 0; m < ncopy; m++) {
kspace->pack_forward_grid(which,buf1,copy[m].npack,copy[m].packlist);
kspace->unpack_forward_grid(which,buf1,copy[m].nunpack,copy[m].unpacklist);
ptr->pack_forward_grid(which,buf1,copy[m].npack,copy[m].packlist);
ptr->unpack_forward_grid(which,buf1,copy[m].nunpack,copy[m].unpacklist);
}
// unpack all received data
@ -1003,8 +1022,8 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which,
for (i = 0; i < nrecv; i++) {
MPI_Waitany(nrecv,requests,&m,MPI_STATUS_IGNORE);
offset = nper * recv[m].offset * nbyte;
kspace->unpack_forward_grid(which,(void *) &buf2[offset],
recv[m].nunpack,recv[m].unpacklist);
ptr->unpack_forward_grid(which,(void *) &buf2[offset],
recv[m].nunpack,recv[m].unpacklist);
}
}
@ -1012,31 +1031,43 @@ forward_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which,
reverse comm of my ghost cells to sum to owner cells
------------------------------------------------------------------------- */
void GridComm::reverse_comm_kspace(KSpace *kspace, int nper, int nbyte, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
void GridComm::reverse_comm(int caller, void *ptr, int nper, int nbyte, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
{
if (layout == REGULAR)
reverse_comm_kspace_regular(kspace,nper,nbyte,which,buf1,buf2,datatype);
else
reverse_comm_kspace_tiled(kspace,nper,nbyte,which,buf1,buf2,datatype);
if (layout == REGULAR) {
if (caller == KSPACE)
reverse_comm_regular<KSpace>((KSpace *)ptr,nper,nbyte,which,
buf1,buf2,datatype);
else if (caller == FIX)
reverse_comm_regular<Fix>((Fix *)ptr,nper,nbyte,which,
buf1,buf2,datatype);
} else {
if (caller == KSPACE)
reverse_comm_tiled<KSpace>((KSpace *)ptr,nper,nbyte,which,
buf1,buf2,datatype);
else if (caller == FIX)
reverse_comm_tiled<Fix>((Fix *)ptr,nper,nbyte,which,
buf1,buf2,datatype);
}
}
/* ----------------------------------------------------------------------
reverse comm on regular grid of procs via list of swaps with 6 neighbor procs
------------------------------------------------------------------------- */
template < class T >
void GridComm::
reverse_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
reverse_comm_regular(T *ptr, int nper, int /*nbyte*/, int which,
void *buf1, void *buf2, MPI_Datatype datatype)
{
int m;
MPI_Request request;
for (m = nswap-1; m >= 0; m--) {
if (swap[m].recvproc == me)
kspace->pack_reverse_grid(which,buf2,swap[m].nunpack,swap[m].unpacklist);
ptr->pack_reverse_grid(which,buf2,swap[m].nunpack,swap[m].unpacklist);
else
kspace->pack_reverse_grid(which,buf1,swap[m].nunpack,swap[m].unpacklist);
ptr->pack_reverse_grid(which,buf1,swap[m].nunpack,swap[m].unpacklist);
if (swap[m].recvproc != me) {
if (swap[m].npack) MPI_Irecv(buf2,nper*swap[m].npack,datatype,
@ -1046,7 +1077,7 @@ reverse_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which,
if (swap[m].npack) MPI_Wait(&request,MPI_STATUS_IGNORE);
}
kspace->unpack_reverse_grid(which,buf2,swap[m].npack,swap[m].packlist);
ptr->unpack_reverse_grid(which,buf2,swap[m].npack,swap[m].packlist);
}
}
@ -1054,9 +1085,10 @@ reverse_comm_kspace_regular(KSpace *kspace, int nper, int /*nbyte*/, int which,
reverse comm on tiled grid decomp via Send/Recv lists of each neighbor proc
------------------------------------------------------------------------- */
template < class T >
void GridComm::
reverse_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which,
void *buf1, void *vbuf2, MPI_Datatype datatype)
reverse_comm_tiled(T *ptr, int nper, int nbyte, int which,
void *buf1, void *vbuf2, MPI_Datatype datatype)
{
int i,m,offset;
@ -1073,15 +1105,15 @@ reverse_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which,
// perform all sends to other procs
for (m = 0; m < nrecv; m++) {
kspace->pack_reverse_grid(which,buf1,recv[m].nunpack,recv[m].unpacklist);
ptr->pack_reverse_grid(which,buf1,recv[m].nunpack,recv[m].unpacklist);
MPI_Send(buf1,nper*recv[m].nunpack,datatype,recv[m].proc,0,gridcomm);
}
// perform all copies to self
for (m = 0; m < ncopy; m++) {
kspace->pack_reverse_grid(which,buf1,copy[m].nunpack,copy[m].unpacklist);
kspace->unpack_reverse_grid(which,buf1,copy[m].npack,copy[m].packlist);
ptr->pack_reverse_grid(which,buf1,copy[m].nunpack,copy[m].unpacklist);
ptr->unpack_reverse_grid(which,buf1,copy[m].npack,copy[m].packlist);
}
// unpack all received data
@ -1089,8 +1121,8 @@ reverse_comm_kspace_tiled(KSpace *kspace, int nper, int nbyte, int which,
for (i = 0; i < nsend; i++) {
MPI_Waitany(nsend,requests,&m,MPI_STATUS_IGNORE);
offset = nper * send[m].offset * nbyte;
kspace->unpack_reverse_grid(which,(void *) &buf2[offset],
send[m].npack,send[m].packlist);
ptr->unpack_reverse_grid(which,(void *) &buf2[offset],
send[m].npack,send[m].packlist);
}
}

16
src/KSPACE/gridcomm.h → src/gridcomm.h
View File

@ -27,8 +27,8 @@ class GridComm : protected Pointers {
virtual ~GridComm();
void setup(int &, int &);
int ghost_adjacent();
void forward_comm_kspace(class KSpace *, int, int, int, void *, void *, MPI_Datatype);
void reverse_comm_kspace(class KSpace *, int, int, int, void *, void *, MPI_Datatype);
void forward_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
void reverse_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
protected:
int me, nprocs;
@ -181,10 +181,14 @@ class GridComm : protected Pointers {
int ghost_adjacent_regular();
int ghost_adjacent_tiled();
void forward_comm_kspace_regular(class KSpace *, int, int, int, void *, void *, MPI_Datatype);
void forward_comm_kspace_tiled(class KSpace *, int, int, int, void *, void *, MPI_Datatype);
void reverse_comm_kspace_regular(class KSpace *, int, int, int, void *, void *, MPI_Datatype);
void reverse_comm_kspace_tiled(class KSpace *, int, int, int, void *, void *, MPI_Datatype);
template <class T>
void forward_comm_regular(T *, int, int, int, void *, void *, MPI_Datatype);
template <class T>
void forward_comm_tiled(T *, int, int, int, void *, void *, MPI_Datatype);
template <class T>
void reverse_comm_regular(T *, int, int, int, void *, void *, MPI_Datatype);
template <class T>
void reverse_comm_tiled(T *, int, int, int, void *, void *, MPI_Datatype);
virtual void grow_swap();
void grow_overlap();