forked from lijiext/lammps
2231 lines
85 KiB
C++
2231 lines
85 KiB
C++
// ATC_Transfer headers ?
|
|
#include "ATC_TransferHardy.h"
|
|
#include "ATC_Error.h"
|
|
#include "FE_Engine.h"
|
|
#include "ATC_HardyKernel.h"
|
|
#include "LammpsInterface.h"
|
|
#include "Quadrature.h"
|
|
|
|
// Other Headers
|
|
#include <vector>
|
|
#include <map>
|
|
#include <set>
|
|
#include <utility>
|
|
#include <fstream>
|
|
#include <exception>
|
|
|
|
static int sgn(double x) { return (int) ((x>0) - (x<0)); }
|
|
static int rnd(double x) { return (int) (x+sgn(x)*0.5); }
|
|
|
|
using namespace std;
|
|
|
|
namespace ATC {
|
|
|
|
ATC_TransferHardy::ATC_TransferHardy(string groupName,
|
|
string matParamFile)
|
|
: ATC_Transfer(),
|
|
setRefPE_(false),
|
|
setRefPEvalue_(false),
|
|
xPointer_(NULL),
|
|
kernelFunction_(NULL),
|
|
cauchyBornStress_(NULL),
|
|
bondOnTheFly_(false),
|
|
kernelOnTheFly_(false),
|
|
useAtomicShapeFunctions_(true)
|
|
{
|
|
// assign default "internal" group
|
|
int igroup = lammpsInterface_->find_group(groupName.c_str());
|
|
groupbit_ |= lammpsInterface_->group_bit(igroup);
|
|
igroups_.insert(igroup);
|
|
|
|
nTypes_ = lammpsInterface_->ntypes();
|
|
|
|
|
|
// Defaults
|
|
simTime_ = 0.;
|
|
|
|
fieldFlags_.reset(NUM_HARDY_FIELDS);
|
|
fieldFlags_ = true;
|
|
gradFlags_.reset(NUM_HARDY_FIELDS);
|
|
gradFlags_ = false;
|
|
rateFlags_.reset(NUM_HARDY_FIELDS);
|
|
rateFlags_ = false;
|
|
|
|
fieldSizes_.reset(NUM_HARDY_FIELDS);
|
|
fieldSizes_ = 1;
|
|
fieldSizes_(HARDY_DISPLACEMENT) =3;
|
|
fieldSizes_(HARDY_MOMENTUM) =3;
|
|
fieldSizes_(HARDY_VELOCITY) =3;
|
|
fieldSizes_(HARDY_PROJECTED_VELOCITY) =3;
|
|
fieldSizes_(HARDY_HEAT_FLUX) =3;
|
|
fieldSizes_(HARDY_STRESS) =9; // LAGRANGIAN
|
|
fieldSizes_(HARDY_ESHELBY_STRESS) =9;
|
|
fieldSizes_(HARDY_CAUCHY_BORN_STRESS) =6;
|
|
fieldSizes_(HARDY_TRANSFORMED_STRESS) =9;
|
|
fieldSizes_(HARDY_TYPE_CONCENTRATION) =nTypes_;
|
|
|
|
atomicOutputMask_.insert("displacement");
|
|
atomicOutputMask_.insert("velocity");
|
|
atomicOutputMask_.insert("temperature");
|
|
atomicOutputMask_.insert("potential_energy");
|
|
atomicOutputMask_.insert("force");
|
|
atomicOutputMask_.insert("virial");
|
|
atomicOutputMask_.insert("centrosymmetry");
|
|
|
|
set_line_quadrature(line_ngauss,line_xg,line_wg);
|
|
|
|
int computeID;
|
|
// create compute for pe/atom
|
|
computeID = lammpsInterface_->atomPE_create();
|
|
if (lammpsInterface_->comm_rank() == 0 ) {
|
|
cout << "atomPE compute created with ID: " << computeID << "\n" << flush;
|
|
}
|
|
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
ATC_TransferHardy::~ATC_TransferHardy()
|
|
{
|
|
// Deallocate
|
|
if (kernelFunction_) delete kernelFunction_;
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// called before the beginning of a "run"
|
|
void ATC_TransferHardy::initialize()
|
|
{
|
|
// Base class initalizations
|
|
if (!initialized_) {
|
|
feEngine_->initialize();
|
|
nNodes_ = feEngine_->get_nNodes();
|
|
nsd_ = feEngine_->get_nsd();
|
|
// set ref positions
|
|
set_xref();
|
|
|
|
// set nlocal & internal to atom map
|
|
reset_nlocal();
|
|
if (useAtomicShapeFunctions_) {
|
|
try {
|
|
reset_shape_functions();
|
|
}
|
|
catch(bad_alloc&) {
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << " ATC:: insufficient memory for reset_shape_functions() to execute\n";
|
|
throw;
|
|
}
|
|
}
|
|
set_atomic_weights();
|
|
}
|
|
|
|
NodeVolumes_.reset(nNodes_,nNodes_);
|
|
invNodeVolumes_.reset(nNodes_,nNodes_);
|
|
feEngine_->compute_lumped_mass_matrix(NodeVolumes_);
|
|
invNodeVolumes_ = inv(NodeVolumes_);
|
|
}
|
|
|
|
// positions
|
|
set_xPointer();
|
|
|
|
// set up workspaces
|
|
if (!initialized_) {
|
|
|
|
nNodesGlobal_ = feEngine_->get_feMesh()->get_nNodes();
|
|
|
|
// get periodicity and box bounds/lengths
|
|
lammpsInterface_->get_box_periodicity(periodicity[0],
|
|
periodicity[1],periodicity[2]);
|
|
lammpsInterface_->get_box_bounds(box_bounds[0][0],box_bounds[1][0],
|
|
box_bounds[0][1],box_bounds[1][1],
|
|
box_bounds[0][2],box_bounds[1][2]);
|
|
|
|
for (int k = 0; k < 3; k++) {
|
|
box_length[k] = box_bounds[1][k] - box_bounds[0][k];
|
|
}
|
|
|
|
if (atomToElementMapType_ == EULERIAN) {
|
|
fieldSizes_(HARDY_STRESS) =6;
|
|
}
|
|
|
|
// PE compute
|
|
lammpsInterface_->atomPE_init();
|
|
|
|
// ground state for PE
|
|
nodalRefPotentialEnergy_.reset(nNodes_,1);
|
|
|
|
// copy requested fieldFlags_ settings to outputFlags_
|
|
outputFlags_ = fieldFlags_;
|
|
|
|
// perform consistency check on fieldFlags
|
|
check_fieldFlags_consistency();
|
|
|
|
// check whether single_enable==0 for stress/heat flux calculation
|
|
if (fieldFlags_(HARDY_STRESS) || fieldFlags_(HARDY_HEAT_FLUX)) {
|
|
if (lammpsInterface_->single_enable()==0) {
|
|
throw ATC_Error(0,"Calculation of Hardy stress field not possible with selected pair type.");
|
|
}
|
|
}
|
|
|
|
// size arrays for requested/required Hardy fields
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
string name;
|
|
int size;
|
|
if (fieldFlags_(index)) {
|
|
hardy_field_to_string(index,name);
|
|
size = fieldSizes_(index);
|
|
hardyData_ [name].reset(nNodes_,size);
|
|
hardyDataOld_ [name].reset(nNodes_,size);
|
|
filteredHardyData_[name].reset(nNodes_,size);
|
|
}
|
|
}
|
|
|
|
// compute table of Hardy bond functions for processors that own atoms
|
|
if ((! bondOnTheFly_)
|
|
&& ( ( fieldFlags_(HARDY_STRESS)
|
|
|| fieldFlags_(HARDY_ESHELBY_STRESS)
|
|
|| fieldFlags_(HARDY_HEAT_FLUX) ) ) ) {
|
|
try {
|
|
compute_bond_matrix();
|
|
}
|
|
catch(bad_alloc&) {
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << "\n ATC:: stress/heat_flux will be computed on-the-fly\n";
|
|
bondOnTheFly_ = true;
|
|
}
|
|
}
|
|
}
|
|
|
|
// compute shape function matrix for kernel functions
|
|
if (kernelFunction_ && (! kernelOnTheFly_)) {
|
|
try{
|
|
compute_kernel_matrix();
|
|
}
|
|
catch(bad_alloc&) {
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << "\n ATC:: kernel will be computed on-the-fly\n";
|
|
kernelOnTheFly_ = true;
|
|
}
|
|
}
|
|
}
|
|
|
|
// compute table of gradient of Hardy functions
|
|
if (fieldFlags_(HARDY_ESHELBY_STRESS)
|
|
|| fieldFlags_(HARDY_CAUCHY_BORN_STRESS)
|
|
|| fieldFlags_(HARDY_VACANCY_CONCENTRATION)) {
|
|
gradFlags_(HARDY_DISPLACEMENT) = true;
|
|
}
|
|
gradientTable_ = vector<SPAR_MAT>(nsd_);
|
|
if (gradFlags_.has_member(true)) {
|
|
compute_gradient_matrix();
|
|
}
|
|
|
|
// construct & initialize filter
|
|
init_filter();
|
|
}
|
|
initialized_ = true;
|
|
|
|
// addstep needs to be done _every_ initialize not just first
|
|
lammpsInterface_->atomPE_addstep(lammpsInterface_->ntimestep()+1);
|
|
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << " ATC:: conversion factor for energy/vol -> stress "
|
|
<< lammpsInterface_->nktv2p() << "\n";
|
|
cout << " ATC:: cutoff radius "
|
|
<< lammpsInterface_->pair_cutoff() << "\n";
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// sets initial values of filtered quantities
|
|
void ATC_TransferHardy::init_filter()
|
|
{
|
|
timeFilters_.reset(NUM_HARDY_FIELDS);
|
|
sampleCounter_ = 0;
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
string name;
|
|
hardy_field_to_string(index,name);
|
|
filteredHardyData_[name] = 0.0;
|
|
timeFilters_(index) = timeFilterManager_.construct();
|
|
}
|
|
}
|
|
|
|
|
|
//-------------------------------------------------------------------
|
|
// called after the end of a "run"
|
|
void ATC_TransferHardy::finish()
|
|
{
|
|
// base class
|
|
ATC_Transfer::finish();
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// this is the parser
|
|
bool ATC_TransferHardy::modify(int narg, char **arg)
|
|
{
|
|
bool match = false;
|
|
double val, dval;
|
|
|
|
// check to see if it is a transfer class command
|
|
if (strcmp(arg[0],"transfer")==0) {
|
|
/*! \page man_hardy_fields fix_modify AtC transfer fields
|
|
\section syntax
|
|
fix_modify AtC transfer fields <all | none> \n
|
|
fix_modify AtC transfer fields <add | delete> <list_of_fields> \n
|
|
- all | none (keyword) = output all or no fields \n
|
|
- add | delete (keyword) = add or delete the listed output fields \n
|
|
- fields (keyword) = \n
|
|
density : mass per unit volume \n
|
|
displacement : displacement vector \n
|
|
momentum : momentum per unit volume \n
|
|
velocity : defined by momentum divided by density \n
|
|
projected_velocity : simple kernel estimation of atomic velocities \n
|
|
temperature : temperature derived from the relative atomic kinetic energy (as done by Hardy) \n
|
|
kinetic_temperature : temperature derived from the full kinetic energy \n
|
|
energy : total energy (potential + kinetic) per unit volume \n
|
|
number_density : simple kernel estimation of number of atoms per unit volume \n
|
|
stress :
|
|
Cauchy stress tensor for eulerian analysis (atom_element_map), or
|
|
1st Piola-Kirchhoff stress tensor for lagrangian analysis \n
|
|
transformed_stress :
|
|
1st Piola-Kirchhoff stress tensor for eulerian analysis (atom_element_map), or
|
|
Cauchy stress tensor for lagrangian analysis \n
|
|
heat_flux : spatial heat flux vector for eulerian,
|
|
or referential heat flux vector for lagrangian \n
|
|
\section examples
|
|
<TT> fix_modify AtC transfer fields add velocity temperature </TT>
|
|
\section description
|
|
Allows modification of the fields calculated and output by the Hardy
|
|
transfer class. The commands are cumulative, e.g.\n
|
|
<TT> fix_modify AtC transfer fields none </TT> \n
|
|
followed by \n
|
|
<TT> fix_modify AtC transfer fields add velocity temperature </TT> \n
|
|
will only output the velocity and temperature fields.
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer, see \ref man_fix_atc.
|
|
Currently, the stress and heat flux formulas are only correct for
|
|
central force potentials, e.g. Lennard-Jones and EAM
|
|
but not Stillinger-Weber.
|
|
\section related
|
|
See \ref man_hardy_gradients , \ref man_hardy_rates and \ref man_hardy_computes
|
|
\section default
|
|
All fields are output by default
|
|
*/
|
|
if (strcmp(arg[1],"fields")==0) {
|
|
if (strcmp(arg[2],"all")==0) {
|
|
fieldFlags_ = true;
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"none")==0) {
|
|
fieldFlags_ = false;
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"add")==0) {
|
|
hardyFieldName field_name;
|
|
for (int i = 3; i < narg; ++i) {
|
|
if (string_to_hardy_field(arg[i],field_name)) {
|
|
fieldFlags_(field_name) = true; }
|
|
else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
}
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"delete")==0) {
|
|
hardyFieldName field_name;
|
|
for (int i = 3; i < narg; ++i) {
|
|
if (string_to_hardy_field(arg[i],field_name)) {
|
|
fieldFlags_(field_name) = false; }
|
|
else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
}
|
|
match = true;
|
|
}
|
|
}
|
|
|
|
/*! \page man_hardy_gradients fix_modify AtC transfer gradients
|
|
\section syntax
|
|
fix_modify AtC transfer gradients <add | delete> <list_of_fields> \n
|
|
- add | delete (keyword) = add or delete the calculation of gradients for the listed output fields \n
|
|
- fields (keyword) = \n
|
|
gradients can be calculated for all fields listed in \ref man_hardy_fields
|
|
|
|
\section examples
|
|
<TT> fix_modify AtC transfer gradients add temperature velocity stress </TT> \n
|
|
<TT> fix_modify AtC transfer gradients delete velocity </TT> \n
|
|
\section description
|
|
Requests calculation and ouput of gradients of the fields from the Hardy
|
|
transfer class. These gradients will be with regard to spatial or material
|
|
coordinate for eulerian or lagrangian analysis, respectively, as specified by
|
|
atom_element_map (see \ref man_atom_element_map )
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer
|
|
( see \ref man_fix_atc )
|
|
\section related
|
|
\section default
|
|
No gradients are calculated by default
|
|
*/
|
|
else if (strcmp(arg[1],"gradients")==0) {
|
|
if (strcmp(arg[2],"add")==0) {
|
|
hardyFieldName field_name;
|
|
for (int i = 3; i < narg; ++i) {
|
|
if (string_to_hardy_field(arg[i],field_name)) {
|
|
gradFlags_(field_name) = true; }
|
|
else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
}
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"delete")==0) {
|
|
hardyFieldName field_name;
|
|
for (int i = 3; i < narg; ++i) {
|
|
if (string_to_hardy_field(arg[i],field_name)) {
|
|
gradFlags_(field_name) = false; }
|
|
else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
}
|
|
match = true;
|
|
}
|
|
}
|
|
|
|
/*! \page man_hardy_rates fix_modify AtC transfer rates
|
|
\section syntax
|
|
fix_modify AtC transfer rates <add | delete> <list_of_fields> \n
|
|
- add | delete (keyword) = add or delete the calculation of rates (time derivatives) for the listed output fields \n
|
|
- fields (keyword) = \n
|
|
rates can be calculated for all fields listed in \ref man_hardy_fields
|
|
|
|
\section examples
|
|
<TT> fix_modify AtC transfer rates add temperature velocity stress </TT> \n
|
|
<TT> fix_modify AtC transfer rates delete stress </TT> \n
|
|
\section description
|
|
Requests calculation and ouput of rates (time derivatives) of the fields from the Hardy
|
|
transfer class. For eulerian analysis (see \ref man_atom_element_map ), these rates
|
|
are the partial time derivatives of the nodal fields, not the full (material) time
|
|
derivatives. \n
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer
|
|
( see \ref man_fix_atc )
|
|
\section related
|
|
\section default
|
|
No rates are calculated by default
|
|
*/
|
|
else if (strcmp(arg[1],"rates")==0) {
|
|
if (strcmp(arg[2],"add")==0) {
|
|
hardyFieldName field_name;
|
|
for (int i = 3; i < narg; ++i) {
|
|
if (string_to_hardy_field(arg[i],field_name)) {
|
|
rateFlags_(field_name) = true; }
|
|
else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
}
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"delete")==0) {
|
|
hardyFieldName field_name;
|
|
for (int i = 3; i < narg; ++i) {
|
|
if (string_to_hardy_field(arg[i],field_name)) {
|
|
rateFlags_(field_name) = false; }
|
|
else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
}
|
|
match = true;
|
|
}
|
|
}
|
|
|
|
/*! \page man_hardy_computes fix_modify AtC transfer computes
|
|
\section syntax
|
|
fix_modify AtC transfer computes <add | delete> [per-atom compute id] <volume | number> \n
|
|
- add | delete (keyword) = add or delete the calculation of an equivalent continuum field
|
|
for the specified per-atom compute as volume or number density quantity \n
|
|
- per-atom compute id = name/id for per-atom compute,
|
|
fields can be calculated for all per-atom computes available from LAMMPS \n
|
|
- volume | number (keyword) = field created is a per-unit-volume quantity
|
|
or a per-atom quantity as weighted by kernel functions \n
|
|
|
|
\section examples
|
|
<TT> compute virial all stress/atom </TT> \n
|
|
<TT> fix_modify AtC transfer computes add virial volume </TT> \n
|
|
<TT> fix_modify AtC transfer computes delete virial </TT> \n
|
|
\n
|
|
<TT> compute centrosymmetry all centro/atom </TT> \n
|
|
<TT> fix_modify AtC transfer computes add centrosymmetry number </TT> \n
|
|
\section description
|
|
Calculates continuum fields corresponding to specified per-atom computes created by LAMMPS \n
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer ( see \ref man_fix_atc ) \n
|
|
Per-atom compute must be specified before corresponding continuum field can be requested \n
|
|
\section related
|
|
See manual page for compute
|
|
\section default
|
|
No defaults exist for this command
|
|
*/
|
|
else if (strcmp(arg[1],"computes")==0) {
|
|
if (strcmp(arg[2],"add")==0) {
|
|
string tag(arg[3]);
|
|
int icompute = lammpsInterface_->find_compute(tag.c_str());
|
|
if (icompute < 0)
|
|
throw ATC_Error(0,"Could not find compute "+tag);
|
|
int normalization = NO_NORMALIZATION;
|
|
if (narg > 4) {
|
|
if (strcmp(arg[4],"volume")==0) {
|
|
normalization = VOLUME_NORMALIZATION;
|
|
}
|
|
else if (strcmp(arg[4],"number")==0) {
|
|
normalization = NUMBER_NORMALIZATION;
|
|
}
|
|
}
|
|
computes_[tag] = normalization;
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"delete")==0) {
|
|
string tag(arg[3]);
|
|
if (computes_.find(tag) != computes_.end()) {
|
|
computes_.erase(tag);
|
|
}
|
|
else {
|
|
throw ATC_Error(0,tag+" compute is not in list");
|
|
}
|
|
match = true;
|
|
}
|
|
}
|
|
|
|
/*! \page man_hardy_kernel fix_modify AtC kernel
|
|
\section syntax
|
|
fix_modify AtC transfer kernel <type> <parameters>
|
|
- type (keyword) = mesh, step, cell, cubic_cylinder, cubic_sphere, quartic_cylinder, quartic_sphere \n
|
|
- parameters :\n
|
|
mesh = none\n
|
|
step = radius (double) \n
|
|
cell = hx, hy, hz (double) or h (double) \n
|
|
cubic_cylinder = radius (double) \n
|
|
cubic_sphere = radius (double) \n
|
|
quartic_cylinder = radius (double) \n
|
|
quartic_sphere = radius (double) \n
|
|
\section examples
|
|
fix_modify AtC transfer kernel cell 1.0 1.0 1.0
|
|
\section description
|
|
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer \n
|
|
For cylinder kernel types, cylindrical axis is assumed to be in z-direction \n
|
|
( see \ref man_fix_atc )
|
|
\section related
|
|
\section default
|
|
Default to the mesh based kernel
|
|
*/
|
|
else if (strcmp(arg[1],"kernel")==0) {
|
|
if (kernelFunction_) delete kernelFunction_;
|
|
if (strcmp(arg[2],"mesh")==0) {
|
|
match = true;
|
|
useAtomicShapeFunctions_ = true;
|
|
}
|
|
else if (strcmp(arg[2],"step")==0) {
|
|
double parameters[1] = {atof(arg[3])};
|
|
kernelFunction_ = new ATC_HardyKernelStep(1,parameters);
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"cell")==0) {
|
|
double parameters[3];
|
|
parameters[0] = parameters[1] = parameters[2] = atof(arg[3]);
|
|
if (narg > 5) {
|
|
for (int i = 1; i < 3; i++) { parameters[i] = atof(arg[3+i]); }
|
|
}
|
|
kernelFunction_ = new ATC_HardyKernelCell(2,parameters);
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"cubic_cylinder")==0) {
|
|
double parameters[1] = {atof(arg[3])}; // cutoff radius
|
|
kernelFunction_ = new ATC_HardyKernelCubicCyl(1,parameters);
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"cubic_sphere")==0) {
|
|
double parameters[1] = {atof(arg[3])}; // cutoff radius
|
|
kernelFunction_ = new ATC_HardyKernelCubicSphere(1,parameters);
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"quartic_cylinder")==0) {
|
|
double parameters[1] = {atof(arg[3])}; // cutoff radius
|
|
kernelFunction_ = new ATC_HardyKernelQuarticCyl(1,parameters);
|
|
match = true;
|
|
}
|
|
else if (strcmp(arg[2],"quartic_sphere")==0) {
|
|
double parameters[1] = {atof(arg[3])}; // cutoff radius
|
|
kernelFunction_ = new ATC_HardyKernelQuarticSphere(1,parameters);
|
|
match = true;
|
|
}
|
|
}
|
|
|
|
/*! \page man_hardy_on_the_fly fix_modify AtC transfer on_the_fly
|
|
\section syntax
|
|
fix_modify AtC transfer on_the_fly <bond | kernel> <optional on | off> \n
|
|
- bond | kernel (keyword) = specifies on-the-fly calculation of bond or kernel
|
|
matrix elements \n
|
|
- on | off (keyword) = activate or discontinue on-the-fly mode \n
|
|
|
|
\section examples
|
|
<TT> fix_modify AtC transfer on_the_fly bond on </TT> \n
|
|
<TT> fix_modify AtC transfer on_the_fly kernel </TT> \n
|
|
<TT> fix_modify AtC transfer on_the_fly kernel off </TT> \n
|
|
\section description
|
|
Overrides normal mode of pre-calculating and storing bond pair-to-node and
|
|
kernel atom-to-node matrices. If activated, will calculate elements of these
|
|
matrices during repeated calls of field computations (i.e. "on-the-fly") and not store them for
|
|
future use. \n
|
|
on flag is optional - if omitted, on_the_fly will be activated for the specified
|
|
matrix. Can be deactivated using off flag. \n
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer
|
|
( see \ref man_fix_atc )
|
|
\section related
|
|
\section default
|
|
By default, on-the-fly calculation is not active (i.e. off). However, code does a memory allocation
|
|
check to determine if it can store all needed bond and kernel matrix elements. If this allocation
|
|
fails, on-the-fly is activated. \n
|
|
*/
|
|
else if (strcmp(arg[1],"on_the_fly")==0) {
|
|
if (strcmp(arg[2],"bond")==0) {
|
|
bondOnTheFly_ = true;
|
|
if (narg > 3 && strcmp(arg[3],"off")==0) bondOnTheFly_ = false;
|
|
}
|
|
else if (strcmp(arg[2],"kernel")==0) {
|
|
kernelOnTheFly_ = true;
|
|
if (narg > 3 && strcmp(arg[3],"off")==0) kernelOnTheFly_ = false;
|
|
}
|
|
else { throw ATC_Error(0,"unsupported on_the_fly type"); }
|
|
match = true;
|
|
}
|
|
|
|
/*! \page man_hardy_set fix_modify AtC set
|
|
\section syntax
|
|
fix_modify AtC transfer set reference_potential_energy <value>
|
|
- value (double) : optional user specified zero point for PE
|
|
\section examples
|
|
<TT> fix_modify AtC transfer set reference_potential_energy </TT> \n
|
|
<TT> fix_modify AtC transfer set reference_potential_energy -0.05 </TT> \n
|
|
\section description
|
|
Used to set various quantities for the post-processing algorithms.
|
|
Currently it only
|
|
sets the zero point for the potential energy density using
|
|
the value provided for all nodes, or from the current
|
|
configuration of the lattice if no value is provided
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer
|
|
( see \ref man_fix_atc )
|
|
\section related
|
|
\section default
|
|
Defaults to lammps zero point i.e. isolated atoms
|
|
*/
|
|
else if (strcmp(arg[1],"set")==0) {
|
|
if (strcmp(arg[2],"reference_potential_energy")==0) {
|
|
if (narg > 3) {
|
|
double value = atof(arg[3]);
|
|
nodalRefPEvalue_ = value;
|
|
setRefPEvalue_ = true;
|
|
}
|
|
else { // NOTE make this part of initialize
|
|
setRefPE_ = true;
|
|
}
|
|
match = true;
|
|
}
|
|
} // end "set"
|
|
|
|
/*! \page man_boundary_integral fix_modify AtC boundary_integral
|
|
\section syntax
|
|
fix_modify AtC transfer boundary_integral [field] faceset [name]
|
|
- field (string) : name of hardy field
|
|
- name (string) : name of faceset
|
|
\section examples
|
|
<TT> fix_modify AtC transfer boundary_integral stress faceset loop1 </TT> \n
|
|
\section description
|
|
Calculates a surface integral of the given field dotted with the
|
|
outward normal of the faces and puts output in the "GLOBALS" file
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer
|
|
( see \ref man_fix_atc )
|
|
\section related
|
|
\section default
|
|
none
|
|
*/
|
|
else if (strcmp(arg[1],"boundary_integral")==0) {
|
|
hardyFieldName field;
|
|
// if (string_to_hardy_field(arg[2],field)) { }
|
|
// else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
if(strcmp(arg[3],"faceset")==0) {
|
|
string name(arg[4]);
|
|
string field(arg[2]);
|
|
const set< pair<int,int> > * faceSet
|
|
= & ( (feEngine_->get_feMesh())->get_faceset(name));
|
|
pair <string,string> pair_name(name,field);
|
|
bndyIntegralData_[pair_name] = faceSet;
|
|
match = true;
|
|
}
|
|
} // end "boundary_integral"
|
|
|
|
/*! \page man_contour_integral fix_modify AtC contour_integral
|
|
\section syntax
|
|
fix_modify AtC transfer contour_integral [field] faceset [name] <axis [x | y | z]>
|
|
- field (string) : name of hardy field
|
|
- name (string) : name of faceset
|
|
- axis (string) : x or y or z
|
|
\section examples
|
|
<TT> fix_modify AtC transfer contour_integral stress faceset loop1 </TT> \n
|
|
\section description
|
|
Calculates a surface integral of the given field dotted with the
|
|
outward normal of the faces and puts output in the "GLOBALS" file
|
|
\section restrictions
|
|
Must be used with the hardy AtC transfer
|
|
( see \ref man_fix_atc )
|
|
\section related
|
|
\section default
|
|
none
|
|
*/
|
|
else if (strcmp(arg[1],"contour_integral")==0) {
|
|
hardyFieldName field;
|
|
// if (string_to_hardy_field(arg[2],field)) { }
|
|
// else { throw ATC_Error(0,"unsupported Hardy field"); }
|
|
if(strcmp(arg[3],"faceset")==0) {
|
|
string name(arg[4]);
|
|
string field(arg[2]);
|
|
const set< pair<int,int> > * faceSet
|
|
= & ( (feEngine_->get_feMesh())->get_faceset(name));
|
|
pair <string,string> pair_name(name,field);
|
|
contourIntegralData_[pair_name] = faceSet;
|
|
match = true;
|
|
}
|
|
} // end "contour_integral"
|
|
}
|
|
|
|
// no match, call base class parser
|
|
if (!match) {
|
|
match = ATC_Transfer::modify(narg, arg);
|
|
}
|
|
|
|
return match;
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// called at the beginning of a timestep
|
|
void ATC_TransferHardy::pre_init_integrate()
|
|
{
|
|
// output initial configuration
|
|
if (stepCounter_ == 0 && outputFrequency_ > 0) {
|
|
double dt = lammpsInterface_->dt();
|
|
time_filter_pre (dt);
|
|
compute_fields();
|
|
// initialize filtered data
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
string name;
|
|
hardy_field_to_string(index,name);
|
|
filteredHardyData_[name] = hardyData_[name];
|
|
}
|
|
time_filter_post(dt);
|
|
compute_boundary_integrals();
|
|
output();
|
|
}
|
|
}
|
|
//-------------------------------------------------------------------
|
|
// called at the end of first half of a timestep
|
|
void ATC_TransferHardy::post_init_integrate()
|
|
{
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// called at the begining of second half timestep
|
|
void ATC_TransferHardy::pre_final_integrate()
|
|
{
|
|
// cases to recompute transfer matrices:
|
|
// (1) if atoms have changed processors
|
|
bool atomSwitch = atomSwitch_;
|
|
int localSwitch = atomSwitch_, globalSwitch = 0;
|
|
lammpsInterface_->int_allmax(&localSwitch,&globalSwitch);
|
|
atomSwitch = globalSwitch;
|
|
// (2) if eulerian and we are at the atom_to_element_map reset frequency
|
|
// which is used to track the convection of atoms
|
|
bool eulerianReset = (atomToElementMapType_ == EULERIAN)
|
|
&& (stepCounter_ % atomToElementMapFrequency_ ==0);
|
|
|
|
bool needsKernel = ( kernelFunction_ && (! kernelOnTheFly_ ) );
|
|
|
|
bool needsBond = (! bondOnTheFly_ ) &&
|
|
(fieldFlags_(HARDY_STRESS)
|
|
|| fieldFlags_(HARDY_ESHELBY_STRESS)
|
|
|| fieldFlags_(HARDY_HEAT_FLUX));
|
|
|
|
// (3) if we suspect neighbor lists changing and are
|
|
// at the neighbor reset frequency
|
|
int ago = lammpsInterface_->neighbor_ago();
|
|
int interval = outputFrequency_;
|
|
if (timeFilterManager_.filter_dynamics() && sampleFrequency_ < interval)
|
|
{ interval = sampleFrequency_; }
|
|
bool neighborReset = (ago <= interval);
|
|
if (neighborResetFrequency_ > 0) {
|
|
neighborReset = (neighborReset
|
|
|| (stepCounter_ % neighborResetFrequency_ == 0) );
|
|
}
|
|
|
|
if (atomSwitch || eulerianReset) {
|
|
reset_nlocal();
|
|
if (needsKernel) {
|
|
compute_kernel_matrix();
|
|
}
|
|
if (useAtomicShapeFunctions_) {
|
|
reset_shape_functions();
|
|
set_atomic_weights();
|
|
}
|
|
if (needsBond) {
|
|
compute_bond_matrix();
|
|
}
|
|
if (gradFlags_.has_member(true)) {
|
|
compute_gradient_matrix();
|
|
}
|
|
}
|
|
else if (neighborReset && needsBond) {
|
|
check_pair_map(); // will recompute bond matrix if necessary
|
|
}
|
|
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// called at the end of second half timestep
|
|
void ATC_TransferHardy::post_final_integrate()
|
|
{
|
|
double dt = lammpsInterface_->dt();
|
|
simTime_ += dt;
|
|
if (dt == 0.0) simTime_ = stepCounter_;
|
|
++stepCounter_;
|
|
|
|
bool output_now = ( (outputFrequency_ > 0)
|
|
&& ((stepCounter_ % outputFrequency_ == 0)) );
|
|
bool sample_now = ( (output_now)
|
|
|| (timeFilterManager_.filter_dynamics() && (stepCounter_ % sampleFrequency_ == 0)));
|
|
|
|
// compute spatially smoothed quantities
|
|
if (sample_now) {
|
|
time_filter_pre (dt, output_now); // NOTE drop output_now
|
|
compute_fields();
|
|
time_filter_post(dt, output_now);
|
|
compute_boundary_integrals();
|
|
}
|
|
|
|
// output
|
|
if ( output_now ) output();
|
|
if ( (outputFrequencyAtom_ > 0)
|
|
&& ((stepCounter_ % outputFrequencyAtom_ == 0)) )
|
|
atomic_output();
|
|
|
|
// "1" should be outputFrequency_ or sampleFrequency_
|
|
lammpsInterface_->atomPE_addstep(lammpsInterface_->ntimestep()+1);
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_fields(void)
|
|
{
|
|
if (setRefPE_) {
|
|
compute_potential_energy(nodalRefPotentialEnergy_);
|
|
setRefPE_ = false;
|
|
}
|
|
|
|
if (setRefPEvalue_) {
|
|
nodalRefPotentialEnergy_ = nodalRefPEvalue_;
|
|
setRefPEvalue_ = false;
|
|
}
|
|
|
|
|
|
// NOTE remove this
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ** v = lammpsInterface_->vatom();
|
|
double ** x = lammpsInterface_->xatom();
|
|
int atomIdx;
|
|
|
|
bool needs_velocity = fieldFlags_(HARDY_VELOCITY);
|
|
bool needs_projected_velocity = fieldFlags_(HARDY_PROJECTED_VELOCITY);
|
|
bool needs_momentum = fieldFlags_(HARDY_MOMENTUM);
|
|
bool needs_density = fieldFlags_(HARDY_DENSITY);
|
|
bool needs_displacement = fieldFlags_(HARDY_DISPLACEMENT);
|
|
bool needs_temperature = fieldFlags_(HARDY_TEMPERATURE);
|
|
bool needs_kinetic_temperature = fieldFlags_(HARDY_KINETIC_TEMPERATURE);
|
|
bool needs_stress = fieldFlags_(HARDY_STRESS);
|
|
bool needs_eshelby_stress = fieldFlags_(HARDY_ESHELBY_STRESS);
|
|
bool needs_cauchy_born_stress = fieldFlags_(HARDY_CAUCHY_BORN_STRESS);
|
|
bool needs_heat = fieldFlags_(HARDY_HEAT_FLUX);
|
|
bool needs_energy = fieldFlags_(HARDY_ENERGY);
|
|
bool needs_number_density = fieldFlags_(HARDY_NUMBER_DENSITY);
|
|
bool needs_transformed_stress = fieldFlags_(HARDY_TRANSFORMED_STRESS);
|
|
bool needs_vacancy_concentration = fieldFlags_(HARDY_VACANCY_CONCENTRATION);
|
|
bool needs_type_concentration = fieldFlags_(HARDY_TYPE_CONCENTRATION);
|
|
|
|
// (1) direct quantities
|
|
//compute : kernel norm is a _weighted_ number density
|
|
if (needs_number_density) {
|
|
compute_number_density(hardyData_["number_density"]);
|
|
}
|
|
// compute: density
|
|
if (needs_density) {
|
|
compute_mass_density(hardyData_["density"]);
|
|
}
|
|
// compute: displacement
|
|
if (needs_displacement) {
|
|
compute_displacement(hardyData_["displacement"],hardyData_["density"]);
|
|
}
|
|
// compute: momentum
|
|
if (needs_momentum) {
|
|
compute_momentum(hardyData_["momentum"]);
|
|
}
|
|
//compute: projected velocity
|
|
if (needs_projected_velocity) {
|
|
compute_projected_velocity(hardyData_["projected_velocity"]);
|
|
}
|
|
// compute: velocity
|
|
if (needs_velocity) {
|
|
compute_velocity(hardyData_["velocity"],
|
|
hardyData_["density"], hardyData_["momentum"]);
|
|
// NOTE make this only for mesh
|
|
compute_variation_velocity(uVariationVelocity_, hardyData_["velocity"]);
|
|
}
|
|
|
|
// compute : temperature & (full) kinetic temperature
|
|
if (needs_temperature) {
|
|
compute_temperature(hardyData_["temperature"]);
|
|
}
|
|
if (needs_kinetic_temperature) {
|
|
compute_kinetic_temperature(hardyData_["kinetic_temperature"]);
|
|
}
|
|
// compute: stress
|
|
if (needs_stress) {
|
|
compute_stress(hardyData_["stress"]);
|
|
}
|
|
// compute: heat flux
|
|
if (needs_heat) {
|
|
compute_heatflux(hardyData_["heat_flux"]);
|
|
}
|
|
// compute: energy
|
|
if (needs_energy) {
|
|
compute_total_energy(hardyData_["energy"]);
|
|
}
|
|
|
|
// (2) derived quantities
|
|
// compute: gradients
|
|
if (gradFlags_.has_member(true)) {
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
if (gradFlags_(index)) {
|
|
string field;
|
|
hardy_field_to_string(index,field);
|
|
string grad_field = field + "_gradient";
|
|
if (hardyData_.find(field) == hardyData_.end() ) {
|
|
throw ATC_Error(0,"field " + field + " needs to be defined for gradient");
|
|
}
|
|
gradient_compute(hardyData_[field], hardyData_[grad_field]);
|
|
}
|
|
}
|
|
}
|
|
|
|
// (3) computes
|
|
map <string,int>::const_iterator iter;
|
|
for (iter = computes_.begin(); iter != computes_.end(); iter++) {
|
|
string tag = iter->first;
|
|
int projection = iter->second;
|
|
int ncols = lammpsInterface_->compute_ncols(tag.c_str());;
|
|
DENS_MAT atomicData(nLocal_,ncols);
|
|
if (ncols == 1) {
|
|
double * atomData = lammpsInterface_->compute_scalar_data(tag.c_str());
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
atomicData(i,0) = atomData[atomIdx];
|
|
}
|
|
}
|
|
else {
|
|
double ** atomData = lammpsInterface_->compute_vector_data(tag.c_str());
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
for (int k = 0; k < ncols; k++) {
|
|
atomicData(i,k) = atomData[atomIdx][k];
|
|
}
|
|
}
|
|
}
|
|
if (projection == VOLUME_NORMALIZATION) {
|
|
project_volume_normalized(atomicData, hardyData_[tag]);
|
|
}
|
|
else if (projection == NUMBER_NORMALIZATION) {
|
|
project_count_normalized(atomicData, hardyData_[tag]);
|
|
}
|
|
else {
|
|
ATC_Transfer::restrict_unscaled(atomicData, hardyData_[tag]);
|
|
}
|
|
}
|
|
|
|
|
|
}// end of compute_fields routine
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_boundary_integrals(void)
|
|
{
|
|
if (! bndyIntegralData_.empty()) {
|
|
map < pair<string,string>, const set< PAIR > * >::const_iterator iter;
|
|
DENS_MAT values;
|
|
for (iter = bndyIntegralData_.begin();
|
|
iter != bndyIntegralData_.end(); iter++) {
|
|
string bndyName = (iter->first).first;
|
|
string fieldName = (iter->first).second;
|
|
const set< PAIR > & faceSet = *(iter->second);
|
|
const DENS_MAT & data = (hardyData_.find(fieldName)->second);
|
|
feEngine_->field_surface_flux(data,faceSet,values);
|
|
for (int i = 0; i < values.nRows() ; ++i ) {
|
|
string name = bndyName + "_" + fieldName + "_"
|
|
+ ATC_STRING::tostring(i+1);
|
|
feEngine_->add_global(name, values(i,0));
|
|
}
|
|
}
|
|
}
|
|
if (! contourIntegralData_.empty()) {
|
|
map < pair<string,string>, const set< PAIR > * >::const_iterator iter;
|
|
DENS_MAT values;
|
|
for (iter = contourIntegralData_.begin();
|
|
iter != contourIntegralData_.end(); iter++) {
|
|
string bndyName = (iter->first).first;
|
|
string fieldName = (iter->first).second;
|
|
const set< PAIR > & faceSet = *(iter->second);
|
|
const DENS_MAT & data = (hardyData_.find(fieldName)->second);
|
|
feEngine_->field_surface_flux(data,faceSet,values,true);
|
|
for (int i = 0; i < values.nRows() ; ++i ) {
|
|
string name = "contour_"+bndyName + "_" + fieldName + "_"
|
|
+ ATC_STRING::tostring(i+1);
|
|
feEngine_->add_global(name, values(i,0));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_potential_energy(DENS_MAT & nodalPE)
|
|
{
|
|
DENS_MAT atomicEnergy(nLocal_,1);
|
|
|
|
// compute pair energy per atom
|
|
double * atomPE = lammpsInterface_->atomPE_compute();
|
|
|
|
// add up potential energy per atom
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
atomicEnergy(i,0) += atomPE[atomIdx];
|
|
}
|
|
project_volume_normalized(atomicEnergy, nodalPE );
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_number_density(DENS_MAT & density)
|
|
{
|
|
DENS_MAT atomCnt(nLocal_,1);
|
|
atomCnt = 1;
|
|
project_volume_normalized(atomCnt, density);
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_mass_density(DENS_MAT & density)
|
|
{
|
|
atomicDensity_.reset(nLocal_,1);
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
// density
|
|
atomicDensity_(i,0) = ma;
|
|
}
|
|
project_volume_normalized(atomicDensity_, density);
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_displacement(DENS_MAT & displacement,
|
|
const DENS_MAT & density)
|
|
{
|
|
atomicDisplacement_.reset(nLocal_,nsd_);
|
|
double ** x = lammpsInterface_->xatom();
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
for (int j = 0; j < nsd_; j++) {
|
|
// displacement
|
|
atomicDisplacement_(i,j) = x[atomIdx][j] - xref_[atomIdx][j];
|
|
// the cheap version of unwrapping atomic positions
|
|
if ((bool) periodicity[j]) {
|
|
double u = atomicDisplacement_(i,j);
|
|
if (u >= 0.5*box_length[j]) { u -= box_length[j]; }
|
|
if (u <= -0.5*box_length[j]) { u += box_length[j]; }
|
|
atomicDisplacement_(i,j) = u;
|
|
}
|
|
atomicDisplacement_(i,j) *= ma;
|
|
}
|
|
}
|
|
project_volume_normalized(atomicDisplacement_, displacement);
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
double rho_i = density(i,0);
|
|
if (rho_i > 0.0) {
|
|
for (int j = 0; j < nsd_; j++) {
|
|
displacement(i,j) = displacement(i,j)/rho_i;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_momentum(DENS_MAT & momentum)
|
|
{
|
|
atomicMomentum_.reset(nLocal_,nsd_);
|
|
double ** v = lammpsInterface_->vatom();
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
// momentum
|
|
for (int j = 0; j < nsd_; j++) {
|
|
atomicMomentum_(i,j) = ma*(v[atomIdx][j]);
|
|
}
|
|
}
|
|
project_volume_normalized(atomicMomentum_, momentum);
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_projected_velocity(DENS_MAT & velocity)
|
|
{
|
|
atomicVelocity_.reset(nLocal_,nsd_);
|
|
double ** v = lammpsInterface_->vatom();
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
for (int j = 0; j < nsd_; j++) {
|
|
// velocity
|
|
atomicVelocity_(i,j) += v[atomIdx][j];
|
|
}
|
|
}
|
|
project_count_normalized(atomicVelocity_, velocity);
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_velocity(DENS_MAT & velocity,
|
|
const DENS_MAT & density,
|
|
const DENS_MAT & momentum)
|
|
{
|
|
velocity.reset(nNodes_,nsd_);
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
double rho_i = density(i,0);
|
|
if (rho_i > 0.0) {
|
|
for (int j = 0; j < nsd_; j++) {
|
|
velocity(i,j) = momentum(i,j)/rho_i;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_variation_velocity(DENS_MAT & velocity,
|
|
const DENS_MAT & vI)
|
|
{
|
|
if (nLocal_>0) {
|
|
// interpolate nodal velocities to the atoms
|
|
vbar_.reset(nLocal_,nsd_);
|
|
double ** v = lammpsInterface_->vatom();
|
|
// use of prolong assumes atom system contained within mesh
|
|
if (useAtomicShapeFunctions_) {
|
|
ATC_Transfer::prolong(vI,vbar_);
|
|
}
|
|
// compute and store variation velocities of atoms
|
|
uVariationVelocity_.reset(nLocal_,nsd_);
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
for (int j = 0; j < nsd_; j++) {
|
|
velocity(i,j) = v[atomIdx][j] - vbar_(i,j);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_temperature(DENS_MAT & temperature)
|
|
{
|
|
atomicTemperature_.reset(nLocal_,1);
|
|
double kB = lammpsInterface_->kBoltzmann();
|
|
double Tcoef = 1./(nsd_*kB);
|
|
double ** v = lammpsInterface_->vatom();
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
for (int j = 0; j < nsd_; j++) {
|
|
// variation temperature
|
|
atomicTemperature_(i,0) += Tcoef*ma*(uVariationVelocity_(i,j)*uVariationVelocity_(i,j));
|
|
}
|
|
}
|
|
project_count_normalized(atomicTemperature_, temperature);
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_kinetic_temperature(DENS_MAT & temperature)
|
|
{
|
|
atomicKineticTemperature_.reset(nLocal_,1);
|
|
double kB = lammpsInterface_->kBoltzmann();
|
|
double Tcoef = 1./(nsd_*kB);
|
|
double ** v = lammpsInterface_->vatom();
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
for (int j = 0; j < nsd_; j++) {
|
|
// full temperature
|
|
atomicKineticTemperature_(i,0) += Tcoef*ma*(v[atomIdx][j]*v[atomIdx][j]);
|
|
}
|
|
}
|
|
project_count_normalized(atomicKineticTemperature_, temperature);
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_stress(DENS_MAT & stress)
|
|
{
|
|
// table of bond functions already calculated in initialize function
|
|
// get conversion factor for mvv to e units
|
|
double mvv2e = lammpsInterface_->mvv2e();
|
|
// get conversion factor for nktV to p units
|
|
double nktv2p = lammpsInterface_->nktv2p();
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
|
|
// calculate kinetic energy tensor part of stress for Eulerian analysis
|
|
DENS_MAT & vI = uVariationVelocity_;
|
|
if (atomToElementMapType_ == EULERIAN && nLocal_>0) {
|
|
atomicStress_.reset(nLocal_,6);
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
// convert mass to appropriate units
|
|
ma = mvv2e*ma;
|
|
atomicStress_(i,0) -= ma*vI(i,0)*vI(i,0);
|
|
atomicStress_(i,1) -= ma*vI(i,1)*vI(i,1);
|
|
atomicStress_(i,2) -= ma*vI(i,2)*vI(i,2);
|
|
atomicStress_(i,3) -= ma*vI(i,0)*vI(i,1);
|
|
atomicStress_(i,4) -= ma*vI(i,0)*vI(i,2);
|
|
atomicStress_(i,5) -= ma*vI(i,1)*vI(i,2);
|
|
}
|
|
project_volume_normalized(atomicStress_, stress);
|
|
}
|
|
else {
|
|
// zero stress table for Lagrangian analysis or if nLocal_ = 0
|
|
stress.zero();
|
|
}
|
|
// add-in potential part of stress tensor
|
|
int nrows = stress.nRows();
|
|
int ncols = stress.nCols();
|
|
DENS_MAT local_potential_hardy_stress(nrows,ncols);
|
|
if (nLocal_>0) {
|
|
if (bondOnTheFly_) {
|
|
compute_potential_stress(local_potential_hardy_stress);
|
|
}
|
|
else {
|
|
// compute table of force & position dyad
|
|
compute_force_matrix();
|
|
// calculate force part of stress tensor
|
|
local_potential_hardy_stress = atomicBondTable_*atomicForceTable_;
|
|
local_potential_hardy_stress *= 0.5; // NOTE: work around for matrix multiplication
|
|
}
|
|
}
|
|
// global summation of potential part of stress tensor
|
|
DENS_MAT potential_hardy_stress(nrows,ncols);
|
|
int count = nrows*ncols;
|
|
lammpsInterface_->allsum(local_potential_hardy_stress.get_ptr(),
|
|
potential_hardy_stress.get_ptr(), count);
|
|
stress += potential_hardy_stress;
|
|
stress = nktv2p*stress;
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_heatflux(DENS_MAT & flux)
|
|
{
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
// calculate kinetic part of heat flux
|
|
if (atomToElementMapType_ == EULERIAN && nLocal_>0) {
|
|
// get conversion factor for mvv to e units
|
|
double mvv2e = lammpsInterface_->mvv2e();
|
|
// compute pair energy per atom
|
|
double * atomPE = lammpsInterface_->atomPE_compute();
|
|
double atomKE, atomEnergy;
|
|
atomicHeat_.reset(nLocal_,3);
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
// convert mass to appropriate units
|
|
ma = mvv2e*ma;
|
|
atomKE = 0.0;
|
|
for (int j = 0; j < nsd_; j++) {
|
|
atomKE += 0.5*ma*(uVariationVelocity_(i,j)*uVariationVelocity_(i,j));
|
|
}
|
|
atomEnergy = atomKE + atomPE[atomIdx];
|
|
for (int j = 0; j < nsd_; j++) {
|
|
atomicHeat_(i,j) += atomEnergy*uVariationVelocity_(i,j);
|
|
}
|
|
}
|
|
project_volume_normalized(atomicHeat_,flux);
|
|
}
|
|
else {
|
|
// zero stress table for Lagrangian analysis or if nLocal_ = 0
|
|
flux.zero();
|
|
}
|
|
// add-in potential part of heat flux vector
|
|
int nrows = flux.nRows();
|
|
int ncols = flux.nCols();
|
|
DENS_MAT local_hardy_heat(nrows,ncols);
|
|
if (nLocal_>0) {
|
|
if (bondOnTheFly_) {
|
|
compute_potential_heatflux(local_hardy_heat);
|
|
}
|
|
else {
|
|
// compute table of heat vectors
|
|
compute_heat_matrix();
|
|
// calculate force/potential-derivative part of heat flux
|
|
local_hardy_heat = atomicBondTable_*atomicHeatTable_;
|
|
}
|
|
}
|
|
// global summation of potential part of heat flux vector
|
|
DENS_MAT hardy_heat(nrows,ncols);
|
|
int count = nrows*ncols;
|
|
lammpsInterface_->allsum(local_hardy_heat.get_ptr(),
|
|
hardy_heat.get_ptr(), count);
|
|
flux += hardy_heat;
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_total_energy(DENS_MAT & energy)
|
|
{
|
|
atomicEnergy_.reset(nLocal_,1);
|
|
// compute pair energy per atom
|
|
double * atomPE = lammpsInterface_->atomPE_compute();
|
|
// get conversion factor for mvv to e units
|
|
double mvv2e = lammpsInterface_->mvv2e();
|
|
int * type = lammpsInterface_->atom_type();
|
|
double * mass = lammpsInterface_->atom_mass();
|
|
double * rmass = lammpsInterface_->atom_rmass();
|
|
double ma;
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int atomIdx = internalToAtom_(i);
|
|
if (mass) ma = mass[type[atomIdx]];
|
|
else ma = rmass[atomIdx];
|
|
// convert mass to appropriate units
|
|
ma = mvv2e*ma;
|
|
// compute kinetic energy per atom
|
|
double atomKE = 0.0;
|
|
for (int k = 0; k < nsd_; k++) {
|
|
atomKE += 0.5*ma*(uVariationVelocity_(i,k)*uVariationVelocity_(i,k));
|
|
}
|
|
// add up total energy per atom
|
|
atomicEnergy_(i,0) += atomPE[atomIdx] + atomKE;
|
|
}
|
|
project_volume_normalized(atomicEnergy_, energy);
|
|
// subtract zero point energy
|
|
energy -= nodalRefPotentialEnergy_;
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::time_filter_pre(double dt, bool output_now)
|
|
{
|
|
sampleCounter_++;
|
|
string name;
|
|
double delta_t = dt*sampleFrequency_;
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
hardy_field_to_string(index,name);
|
|
timeFilters_(index)->apply_pre_step1(filteredHardyData_[name],
|
|
hardyData_[name], delta_t);
|
|
}
|
|
// NOTE add computes_ here
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::time_filter_post(double dt, bool output_now)
|
|
{
|
|
sampleCounter_++;
|
|
string name;
|
|
double delta_t = dt*sampleFrequency_;
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
hardy_field_to_string(index,name);
|
|
timeFilters_(index)->apply_post_step2(filteredHardyData_[name],
|
|
hardyData_[name], delta_t);
|
|
}
|
|
if (rateFlags_.has_member(true)) {
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
if (rateFlags_(index)) {
|
|
string field;
|
|
hardy_field_to_string(index,field);
|
|
string rate_field = field + "_rate";
|
|
timeFilters_(index)->rate(hardyData_[rate_field],
|
|
filteredHardyData_[field],
|
|
hardyData_[field], delta_t);
|
|
}
|
|
}
|
|
}
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
hardy_field_to_string(index,name);
|
|
if (rateFlags_(index)) {
|
|
string rate_field = name + "_rate";
|
|
filteredHardyData_[rate_field] = hardyData_[rate_field];
|
|
}
|
|
if (gradFlags_(index)) {
|
|
string grad_field = name + "_gradient";
|
|
filteredHardyData_[grad_field] = hardyData_[grad_field];
|
|
}
|
|
}
|
|
// NOTE add computes_ here
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::output()
|
|
{
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
// data
|
|
OUTPUT_LIST output_data;
|
|
for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
|
|
string name;
|
|
hardy_field_to_string(index,name);
|
|
if (outputFlags_(index)) {
|
|
output_data[name] = & ( filteredHardyData_[name]);
|
|
}
|
|
if (rateFlags_(index)) {
|
|
string rate_name = name + "_rate";
|
|
output_data[rate_name] = & ( filteredHardyData_[rate_name]);
|
|
}
|
|
if (gradFlags_(index)) {
|
|
string grad_name = name + "_gradient";
|
|
output_data[grad_name] = & ( filteredHardyData_[grad_name]);
|
|
}
|
|
}
|
|
map <string,int>::const_iterator iter;
|
|
for (iter = computes_.begin(); iter != computes_.end(); iter++) {
|
|
string tag = iter->first;
|
|
output_data[tag] = & ( hardyData_[tag]); // NOTE should be filtered
|
|
}
|
|
|
|
// output
|
|
feEngine_->write_data(simTime_, & output_data);
|
|
}
|
|
}
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_kernel_matrix(void)
|
|
{
|
|
if (nLocal_>0) {
|
|
set_xPointer();
|
|
kernelShpFcn_.reset(nNodes_,nLocal_);
|
|
DENS_VEC xI(nsd_),xa(nsd_),xaI(nsd_);
|
|
if (lammpsInterface_->comm_rank() == 0 ) {
|
|
cout << "ATC:: computing kernel matrix " << flush;
|
|
}
|
|
int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
if (lammpsInterface_->comm_rank() == 0 && i % heartbeatFreq == 0 ) {
|
|
cout << "." << flush;
|
|
}
|
|
// Hardy point
|
|
xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
|
|
//xI = (feEngine_->get_feMesh())->global_coordinates(i);
|
|
int ii = i;
|
|
//int ii = feEngine_->get_feMesh()->map_global_to_unique(i);
|
|
for (int j = 0; j < nLocal_; j++) {
|
|
// atom location
|
|
int lammps_j = internalToAtom_(j);
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
// difference vector
|
|
xaI = xa - xI;
|
|
periodicity_correction(xaI);
|
|
// compute kernel value & add it to the matrix
|
|
double val = kernelFunction_->value(xaI);
|
|
if (val > 0) kernelShpFcn_.add(ii,j,val);
|
|
}
|
|
}
|
|
kernelShpFcn_.compress();
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << "done\n";
|
|
}
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// count normalized
|
|
void ATC_TransferHardy::project_count_normalized(const DENS_MAT & atomData,
|
|
DENS_MAT & nodeData)
|
|
{
|
|
int ncols_atomData = atomData.nCols();
|
|
nodeData.reset(nNodes_,ncols_atomData);
|
|
nodeData.zero();
|
|
// kernel functions
|
|
if (kernelFunction_ ) {
|
|
project_volume_normalized(atomData,nodeData);
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
double numDens = hardyData_["number_density"](i,0);
|
|
if (numDens > 0) {
|
|
double invNumDens = 1.0/numDens;
|
|
for (int j = 0; j < ncols_atomData; j++) {
|
|
nodeData(i,j) *= invNumDens;
|
|
}
|
|
}
|
|
else {
|
|
for (int j = 0; j < ncols_atomData; j++) {
|
|
nodeData(i,j) = 0;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// mesh-based kernel functions
|
|
else {
|
|
ATC_Transfer::restrict(atomData,nodeData);
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// volume normalized
|
|
void ATC_TransferHardy::project_volume_normalized(const DENS_MAT & atomData,
|
|
DENS_MAT & nodeData)
|
|
{
|
|
int ncols_atomData = atomData.nCols();
|
|
DENS_MAT workNodeArray(nNodes_,ncols_atomData);
|
|
workNodeArray.zero();
|
|
nodeData.reset(workNodeArray.nRows(),workNodeArray.nCols());
|
|
nodeData.zero();
|
|
// kernel functions
|
|
if (kernelFunction_ ) {
|
|
if (nLocal_>0) {
|
|
double invVol = kernelFunction_->inv_vol();
|
|
// on the fly calculation
|
|
if (kernelOnTheFly_) {
|
|
set_xPointer();
|
|
DENS_VEC xI(nsd_),xa(nsd_),xaI(nsd_);
|
|
double val;
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
|
|
for (int j = 0; j < nLocal_; j++) {
|
|
int lammps_j = internalToAtom_(j);
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
xaI = xa - xI;
|
|
periodicity_correction(xaI);
|
|
val = kernelFunction_->value(xaI);
|
|
if (val > 0) {
|
|
for (int k=0; k < ncols_atomData; k++) {
|
|
workNodeArray(i,k) += val*atomData(j,k);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// matrix calculation
|
|
else {
|
|
workNodeArray = kernelShpFcn_*atomData;
|
|
}
|
|
workNodeArray *= invVol;
|
|
}
|
|
}
|
|
// mesh-based kernel functions
|
|
else {
|
|
// computes nodeData = N*atomData w/ N: volume normalized shape functions
|
|
if (nLocal_>0)
|
|
workNodeArray = invNodeVolumes_*shpFcn_.transMat(atomData);
|
|
}
|
|
// accumulate across processors
|
|
int count = workNodeArray.nRows()*workNodeArray.nCols();
|
|
lammpsInterface_->allsum(workNodeArray.get_ptr(),nodeData.get_ptr(),count);
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::gradient_compute(const DENS_MAT & inNodeData,
|
|
DENS_MAT & outNodeData)
|
|
{
|
|
int nrows = inNodeData.nRows();
|
|
int ncols = inNodeData.nCols();
|
|
outNodeData.reset(nrows,ncols*nsd_);
|
|
for (int j = 0; j < nNodes_; j++) {
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
int index = 0;
|
|
for (int n = 0; n < ncols; n++) { //output v1,1 v1,2 v1,3 ...
|
|
for (int m = 0; m < nsd_; m++) {
|
|
outNodeData(j,index) += gradientTable_[m](j,i)*inNodeData(i,n);
|
|
index++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_pair_map()
|
|
{
|
|
// neighbor lists for "internal" group
|
|
// a: a "real" atom in the "internal" group
|
|
// b: any neighbor to a (it should project into the mesh for mesh-based)
|
|
int inum = lammpsInterface_->neighbor_list_inum();
|
|
if (inum != nLocalTotal_)
|
|
throw ATC_Error(0,"size mismatch in neighbor list");
|
|
int *ilist = lammpsInterface_->neighbor_list_ilist();
|
|
int *numneigh = lammpsInterface_->neighbor_list_numneigh();
|
|
int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
|
|
int * mask = lammpsInterface_->atom_mask();
|
|
|
|
pairMap_.clear();
|
|
int pair_index = 0;
|
|
pair< int,int > pair_ij;
|
|
for (int i = 0; i < nLocalTotal_; i++) { // NOTE why nLocal__Total__?
|
|
int lammps_i = ilist[i];
|
|
// filter out atoms not in the internal group
|
|
if (mask[lammps_i] & groupbit_) {
|
|
//int tag_i = (lammpsInterface_->atom_tag())[lammps_i];
|
|
for (int j = 0; j < numneigh[lammps_i]; j++) {
|
|
int lammps_j = firstneigh[lammps_i][j];
|
|
//int tag_j = (lammpsInterface_->atom_tag())[lammps_j];
|
|
//if (lammps_i > lammps_j) continue; // skip a > b pairs
|
|
//if (tag_i > tag_j) continue; // skip a > b pairs
|
|
pair_ij.first = lammps_i; // alpha
|
|
pair_ij.second = lammps_j; // beta
|
|
pairMap_[pair_ij] = pair_index;
|
|
pair_index++;
|
|
}
|
|
}
|
|
}
|
|
nPairs_ = pair_index;
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_bond_matrix()
|
|
{
|
|
set_xPointer();
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << " ATC:: computing bond matrix " << flush;
|
|
}
|
|
// compute pair map
|
|
compute_pair_map();
|
|
|
|
// neighbor lists
|
|
int *numneigh = lammpsInterface_->neighbor_list_numneigh();
|
|
int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
|
|
|
|
Array<bool> latticePeriodicity(3);
|
|
latticePeriodicity(0) = (bool) periodicity[0];
|
|
latticePeriodicity(1) = (bool) periodicity[1];
|
|
latticePeriodicity(2) = (bool) periodicity[2];
|
|
|
|
|
|
double lam1,lam2;
|
|
double bond_value;
|
|
pair< int,int > pair_jk;
|
|
map< pair< int,int >,int >::iterator pairMapIterator;
|
|
int pair_index;
|
|
|
|
|
|
atomicBondTable_.reset(nNodes_,nPairs_);
|
|
|
|
// process differently for mesh vs translation-invariant kernels
|
|
if (kernelFunction_ ) {
|
|
int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
|
|
//int heartbeatFreq = (nNodesGlobal_ <= 10 ? 1 : (int) nNodesGlobal_ / 10);
|
|
// "normal" kernel functions
|
|
DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
|
|
double kernel_inv_vol = kernelFunction_->inv_vol();
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
//for (int i = 0; i < nNodesGlobal_; i++) {
|
|
if (lammpsInterface_->comm_rank() == 0 && i % heartbeatFreq == 0 ) {
|
|
cout << "." << flush;
|
|
}
|
|
int ii = i;
|
|
//int ii = feEngine_->get_feMesh()->map_global_to_unique(i);
|
|
// Hardy point
|
|
xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
|
|
//xI = (feEngine_->get_feMesh())->global_coordinates(i);
|
|
for (pairMapIterator = pairMap_.begin();
|
|
pairMapIterator != pairMap_.end(); pairMapIterator++){
|
|
int lammps_j = (pairMapIterator->first).first ;
|
|
int lammps_k = (pairMapIterator->first).second;
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
xaI = xa - xI;
|
|
periodicity_correction(xaI);
|
|
xb.copy(xPointer_[lammps_k],3);
|
|
xba = xb - xa;
|
|
xbI = xba + xaI;
|
|
// intersect segment with kernel support
|
|
kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
|
|
if (lam1 < lam2) {
|
|
bond_value
|
|
= kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
|
|
pair_index = pairMapIterator->second;
|
|
atomicBondTable_.add(ii,pair_index,bond_value);
|
|
} // if lam1 < lam2
|
|
} // pair map
|
|
}// end nodes loop
|
|
atomicBondTable_.compress();
|
|
}
|
|
else {
|
|
// mesh-based kernel functions
|
|
int heartbeatFreq = (int) pairMap_.size() / 10;
|
|
int i=0;
|
|
int nodes_per_element = feEngine_->get_feMesh()->get_nNodesPerElement();
|
|
Array<int> node_list(nodes_per_element);
|
|
DENS_VEC shp(nodes_per_element);
|
|
DENS_VEC xa(nsd_),xb(nsd_),xab(nsd_),xlambda(nsd_);
|
|
double **xatom = lammpsInterface_->xatom();
|
|
for (pairMapIterator = pairMap_.begin();
|
|
pairMapIterator != pairMap_.end(); pairMapIterator++){
|
|
if (lammpsInterface_->comm_rank() == 0 && i++ % heartbeatFreq == 0 ) {
|
|
cout << "." << flush;
|
|
}
|
|
int lammps_j = (pairMapIterator->first).first ;
|
|
int lammps_k = (pairMapIterator->first).second;
|
|
pair_index = pairMapIterator->second;
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
xb.copy(xPointer_[lammps_k],3);
|
|
lam1 = 0.0; lam2 = 1.0;
|
|
double del_lambda = 0.5*(lam2 - lam1);
|
|
double avg_lambda = 0.5*(lam2 + lam1);
|
|
xab = xa - xb;
|
|
for (int i = 0; i < line_ngauss; i++) {
|
|
double lambda = del_lambda*line_xg[i] +avg_lambda;
|
|
xlambda = lambda*xab + xb;
|
|
// NOTE doesn't work for the case that fe_region \subset md_region
|
|
int dummyEltID;
|
|
feEngine_->shape_functions(xlambda,shp,dummyEltID,
|
|
node_list,latticePeriodicity);
|
|
// accumulate to nodes whose support overlaps the integration point
|
|
for (int I = 0; I < nodes_per_element; I++) {
|
|
// Use factor of 0.5 to account for line integration
|
|
// domain of 0 to 1 instead of -1 to 1
|
|
int inode = node_list(I);
|
|
double val = invNodeVolumes_(inode,inode)*shp(I)*line_wg[i]*0.5;
|
|
atomicBondTable_.add(inode,pair_index,val);
|
|
}
|
|
}
|
|
}
|
|
atomicBondTable_.compress();
|
|
}// end if kernelFunction_
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << "done\n";
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_gradient_matrix()
|
|
{
|
|
feEngine_->compute_gradient_matrix(gradientTable_);
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// checks if atoms have changed neighbors
|
|
void ATC_TransferHardy::check_pair_map()
|
|
{
|
|
// (1) compute fresh temporary pairMap
|
|
// (2) check if tmpPairMap == pairMap_
|
|
|
|
// neighbor lists
|
|
int *numneigh = lammpsInterface_->neighbor_list_numneigh();
|
|
int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
|
|
|
|
pair< int,int > pair_ij;
|
|
map< pair< int,int >,int >::iterator pairMapIterator;
|
|
|
|
for (int i = 0; i < nLocal_; i++) {
|
|
int lammps_i = internalToAtom_(i);
|
|
for (int k = 0; k < numneigh[lammps_i]; k++) {
|
|
int j = firstneigh[lammps_i][k];
|
|
pair_ij.first = lammps_i;
|
|
pair_ij.second = j;
|
|
pairMapIterator = pairMap_.find(pair_ij);
|
|
// if pair is not found in pairMap_, recompute pair map and bond table
|
|
// NOTE need to check if pairMap_ has a stale pair <<<<
|
|
// construct new pair map and use "=="
|
|
if (pairMapIterator == pairMap_.end()) {
|
|
cout << "ATC::check_pair_map has found that pair map and bond table need to be recomputed" << endl;
|
|
compute_bond_matrix();
|
|
return;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_force_matrix()
|
|
{
|
|
// reset force table
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
atomicForceTable_.reset(nPairs_,9);
|
|
}
|
|
else if (atomToElementMapType_ == EULERIAN) {
|
|
atomicForceTable_.reset(nPairs_,6);
|
|
}
|
|
// calculate force table
|
|
double **xatom = lammpsInterface_->xatom();
|
|
map< pair< int,int >,int >::iterator pairMapIterator;
|
|
for (pairMapIterator = pairMap_.begin();
|
|
pairMapIterator != pairMap_.end(); pairMapIterator++){
|
|
int lammps_i = (pairMapIterator->first).first ;
|
|
int lammps_j = (pairMapIterator->first).second;
|
|
int pair_index = pairMapIterator->second;
|
|
double * xi = xatom[lammps_i];
|
|
double * xj = xatom[lammps_j];
|
|
double delx = xi[0] - xj[0];
|
|
double dely = xi[1] - xj[1];
|
|
double delz = xi[2] - xj[2];
|
|
double rsq = delx*delx + dely*dely + delz*delz;
|
|
double fforce = 0;
|
|
lammpsInterface_->pair_force(lammps_i,lammps_j,rsq,fforce);
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
double delX = xref_[lammps_i][0] - xref_[lammps_j][0];
|
|
double delY = xref_[lammps_i][1] - xref_[lammps_j][1];
|
|
double delZ = xref_[lammps_i][2] - xref_[lammps_j][2];
|
|
atomicForceTable_(pair_index,0)=-delx*fforce*delX;
|
|
atomicForceTable_(pair_index,1)=-delx*fforce*delY;
|
|
atomicForceTable_(pair_index,2)=-delx*fforce*delZ;
|
|
atomicForceTable_(pair_index,3)=-dely*fforce*delX;
|
|
atomicForceTable_(pair_index,4)=-dely*fforce*delY;
|
|
atomicForceTable_(pair_index,5)=-dely*fforce*delZ;
|
|
atomicForceTable_(pair_index,6)=-delz*fforce*delX;
|
|
atomicForceTable_(pair_index,7)=-delz*fforce*delY;
|
|
atomicForceTable_(pair_index,8)=-delz*fforce*delZ;
|
|
}
|
|
else if (atomToElementMapType_ == EULERIAN) {
|
|
atomicForceTable_(pair_index,0)=-delx*delx*fforce;
|
|
atomicForceTable_(pair_index,1)=-dely*dely*fforce;
|
|
atomicForceTable_(pair_index,2)=-delz*delz*fforce;
|
|
atomicForceTable_(pair_index,3)=-delx*dely*fforce;
|
|
atomicForceTable_(pair_index,4)=-delx*delz*fforce;
|
|
atomicForceTable_(pair_index,5)=-dely*delz*fforce;
|
|
}
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
void ATC_TransferHardy::compute_heat_matrix()
|
|
{
|
|
// reset heat flux and energy tables
|
|
atomicHeatTable_.reset(nPairs_,3);
|
|
// calculate force table
|
|
map< pair< int,int >,int >::iterator pairMapIterator;
|
|
double ** xatom = lammpsInterface_->xatom();
|
|
double fforce;
|
|
for (pairMapIterator = pairMap_.begin();
|
|
pairMapIterator != pairMap_.end(); pairMapIterator++){
|
|
int lammps_i = (pairMapIterator->first).first ;
|
|
int lammps_j = (pairMapIterator->first).second;
|
|
int pair_index = pairMapIterator->second;
|
|
double * xi = xatom[lammps_i];
|
|
double * xj = xatom[lammps_j];
|
|
double energy_i = 0.0;
|
|
double delx = xi[0] - xj[0];
|
|
double dely = xi[1] - xj[1];
|
|
double delz = xi[2] - xj[2];
|
|
double rsq = delx*delx + dely*dely + delz*delz;
|
|
double pair_pe=lammpsInterface_->pair_force(lammps_i,lammps_j,rsq,fforce);
|
|
fforce *= 0.5;
|
|
// This term is correct ONLY for pair potentials.
|
|
// fforce = - (1/x_{ij})*(d(phi_i)/d(x_{ij}) + d(phi_j)/d(x_{ij}))
|
|
// What we need for the heat flux calculation is the portion of this that belongs only to
|
|
// atom j: - (1/x_{ij})*d(phi_j)/d(x_{ij}).
|
|
int i = lammps_i; // NOTE is this right?
|
|
fforce *= (delx*uVariationVelocity_(i,0) +
|
|
dely*uVariationVelocity_(i,1) +
|
|
delz*uVariationVelocity_(i,2));
|
|
// NOTE REJ don't think this is right, where is energy_i from above?
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
double delX = xref_[lammps_i][0] - xref_[lammps_j][0];
|
|
double delY = xref_[lammps_i][1] - xref_[lammps_j][1];
|
|
double delZ = xref_[lammps_i][2] - xref_[lammps_j][2];
|
|
atomicHeatTable_(pair_index,0)=delX*fforce;
|
|
atomicHeatTable_(pair_index,1)=delY*fforce;
|
|
atomicHeatTable_(pair_index,2)=delZ*fforce;
|
|
}
|
|
else if (atomToElementMapType_ == EULERIAN) {
|
|
atomicHeatTable_(pair_index,0)=delx*fforce;
|
|
atomicHeatTable_(pair_index,1)=dely*fforce;
|
|
atomicHeatTable_(pair_index,2)=delz*fforce;
|
|
}
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// on-the-fly calculation of stress
|
|
void ATC_TransferHardy::compute_potential_stress(DENS_MAT& stress)
|
|
{
|
|
set_xPointer();
|
|
stress.zero();
|
|
// neighbor lists
|
|
int *numneigh = lammpsInterface_->neighbor_list_numneigh();
|
|
int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
|
|
double ** xatom = lammpsInterface_->xatom();
|
|
Array<bool> latticePeriodicity(3);
|
|
latticePeriodicity(0) = (bool) periodicity[0];
|
|
latticePeriodicity(1) = (bool) periodicity[1];
|
|
latticePeriodicity(2) = (bool) periodicity[2];
|
|
double lam1,lam2;
|
|
double bond_value;
|
|
// process differently for mesh vs translation-invariant kernels
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << "ATC:: computing potential stress: " << flush;
|
|
}
|
|
int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
|
|
if (kernelFunction_ ) {
|
|
// "normal" kernel functions
|
|
DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
|
|
double kernel_inv_vol = kernelFunction_->inv_vol();
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
//for (int i = 0; i < nNodesGlobal_; i++) {
|
|
if (lammpsInterface_->comm_rank() == 0 && i % heartbeatFreq == 0 ) {
|
|
cout << "." << flush;
|
|
}
|
|
// Hardy point
|
|
xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
|
|
//xI = (feEngine_->get_feMesh())->global_coordinates(i);
|
|
int inode = i;
|
|
//int inode = feEngine_->get_feMesh()->map_global_to_unique(i);
|
|
for (int j = 0; j < nLocal_; j++) {
|
|
// second (neighbor) atom location
|
|
int lammps_j = internalToAtom_(j);
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
double * xj = xatom[lammps_j];
|
|
// difference vector
|
|
xaI = xa - xI;
|
|
periodicity_correction(xaI);
|
|
for (int k = 0; k < numneigh[lammps_j]; ++k) {
|
|
int lammps_k = firstneigh[lammps_j][k];
|
|
// first atom location
|
|
xb.copy(xPointer_[lammps_k],3);
|
|
double * xk = xatom[lammps_k];
|
|
// difference vector
|
|
xba = xb - xa;
|
|
xbI = xba + xaI;
|
|
kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
|
|
// compute virial
|
|
if (lam1 < lam2) {
|
|
bond_value
|
|
= kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
|
|
double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
|
|
double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
|
|
double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
|
|
double rsq = delx*delx + dely*dely + delz*delz;
|
|
double fforce = 0;
|
|
lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
|
|
fforce *= 0.5; // dbl count
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
|
|
double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
|
|
double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
|
|
stress(inode,0) +=-delx*fforce*delX*bond_value;
|
|
stress(inode,1) +=-delx*fforce*delY*bond_value;
|
|
stress(inode,2) +=-delx*fforce*delZ*bond_value;
|
|
stress(inode,3) +=-dely*fforce*delX*bond_value;
|
|
stress(inode,4) +=-dely*fforce*delY*bond_value;
|
|
stress(inode,5) +=-dely*fforce*delZ*bond_value;
|
|
stress(inode,6) +=-delz*fforce*delX*bond_value;
|
|
stress(inode,7) +=-delz*fforce*delY*bond_value;
|
|
stress(inode,8) +=-delz*fforce*delZ*bond_value;
|
|
}
|
|
else { //if (atomToElementMapType_ == EULERIAN) {
|
|
stress(inode,0) +=-delx*delx*fforce*bond_value;
|
|
stress(inode,1) +=-dely*dely*fforce*bond_value;
|
|
stress(inode,2) +=-delz*delz*fforce*bond_value;
|
|
stress(inode,3) +=-delx*dely*fforce*bond_value;
|
|
stress(inode,4) +=-delx*delz*fforce*bond_value;
|
|
stress(inode,5) +=-dely*delz*fforce*bond_value;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
// mesh-based kernel functions
|
|
int nodes_per_element = feEngine_->get_feMesh()->get_nNodesPerElement();
|
|
Array<int> node_list(nodes_per_element);
|
|
DENS_VEC shp(nodes_per_element);
|
|
DENS_VEC xa(nsd_),xb(nsd_),xab(nsd_),xlambda(nsd_);
|
|
for (int j = 0; j < nLocal_; j++) {
|
|
if (lammpsInterface_->comm_rank() == 0 && j % heartbeatFreq == 0 ) {
|
|
cout << "." << flush;
|
|
}
|
|
// first atom location
|
|
int lammps_j = internalToAtom_(j);
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
for (int k = 0; k < numneigh[lammps_j]; ++k) {
|
|
int lammps_k = firstneigh[lammps_j][k];
|
|
//if (lammps_k < lammps_j) continue; // full neighbor list
|
|
// second (neighbor) atom location
|
|
xb.copy(xPointer_[lammps_k],3);
|
|
double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
|
|
double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
|
|
double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
|
|
double rsq = delx*delx + dely*dely + delz*delz;
|
|
double fforce = 0;
|
|
lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
|
|
fforce *= 0.5; // 1/2 sum_ab = sum_(ab)
|
|
double virial[9];
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
|
|
double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
|
|
double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
|
|
virial[0] =-delx*fforce*delX;
|
|
virial[1] =-delx*fforce*delY;
|
|
virial[2] =-delx*fforce*delZ;
|
|
virial[3] =-dely*fforce*delX;
|
|
virial[4] =-dely*fforce*delY;
|
|
virial[5] =-dely*fforce*delZ;
|
|
virial[6] =-delz*fforce*delX;
|
|
virial[7] =-delz*fforce*delY;
|
|
virial[8] =-delz*fforce*delZ;
|
|
}
|
|
else {//if (atomToElementMapType_ == EULERIAN) {
|
|
virial[0] =-delx*delx*fforce;
|
|
virial[1] =-dely*dely*fforce;
|
|
virial[2] =-delz*delz*fforce;
|
|
virial[3] =-delx*dely*fforce;
|
|
virial[4] =-delx*delz*fforce;
|
|
virial[5] =-dely*delz*fforce;
|
|
}
|
|
lam1 = 0.0; lam2 = 1.0;
|
|
double del_lambda = 0.5*(lam2 - lam1);
|
|
double avg_lambda = 0.5*(lam2 + lam1);
|
|
xab = xa - xb;
|
|
for (int i = 0; i < line_ngauss; i++) {
|
|
double lambda = del_lambda*line_xg[i] +avg_lambda;
|
|
xlambda = lambda*xab + xb;
|
|
// NOTE doesn't work for the case that fe_region \subset md_region
|
|
int dummyEltID;
|
|
feEngine_->shape_functions(xlambda,shp,dummyEltID,
|
|
node_list,latticePeriodicity);
|
|
// accumulate to nodes whose support overlaps the integration point
|
|
for (int I = 0; I < nodes_per_element; I++) {
|
|
// Use factor of 0.5 to account for line integration
|
|
// domain of 0 to 1 instead of -1 to 1
|
|
int inode = node_list(I);
|
|
double bond_value
|
|
= invNodeVolumes_(inode,inode)*shp(I)*line_wg[i]*0.5;
|
|
stress(inode,0) += virial[0]*bond_value;
|
|
stress(inode,1) += virial[1]*bond_value;
|
|
stress(inode,2) += virial[2]*bond_value;
|
|
stress(inode,3) += virial[3]*bond_value;
|
|
stress(inode,4) += virial[4]*bond_value;
|
|
stress(inode,5) += virial[5]*bond_value;
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
stress(inode,6) += virial[6]*bond_value;
|
|
stress(inode,7) += virial[7]*bond_value;
|
|
stress(inode,8) += virial[8]*bond_value;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}// end if kernelFunction_
|
|
if (lammpsInterface_->comm_rank() == 0) {
|
|
cout << "done\n" << flush;
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// on-the-fly calculation of the heat flux
|
|
void ATC_TransferHardy::compute_potential_heatflux(DENS_MAT& flux)
|
|
{
|
|
set_xPointer();
|
|
flux.zero();
|
|
// neighbor lists
|
|
int *numneigh = lammpsInterface_->neighbor_list_numneigh();
|
|
int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
|
|
double ** xatom = lammpsInterface_->xatom();
|
|
Array<bool> latticePeriodicity(3);
|
|
latticePeriodicity(0) = (bool) periodicity[0];
|
|
latticePeriodicity(1) = (bool) periodicity[1];
|
|
latticePeriodicity(2) = (bool) periodicity[2];
|
|
double lam1,lam2;
|
|
double bond_value;
|
|
// process differently for mesh vs translation-invariant kernels
|
|
if (kernelFunction_ ) {
|
|
// "normal" kernel functions
|
|
DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
|
|
double kernel_inv_vol = kernelFunction_->inv_vol();
|
|
for (int i = 0; i < nNodes_; i++) {
|
|
//for (int i = 0; i < nNodesGlobal_; i++) {
|
|
int inode = i;
|
|
//int inode = feEngine_->get_feMesh()->map_global_to_unique(i);
|
|
// Hardy point
|
|
xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
|
|
//xI = (feEngine_->get_feMesh())->global_coordinates(i);
|
|
for (int j = 0; j < nLocal_; j++) {
|
|
// second (neighbor) atom location
|
|
int lammps_j = internalToAtom_(j);
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
double * xj = xatom[lammps_j];
|
|
// difference vector
|
|
xaI = xa - xI;
|
|
periodicity_correction(xaI);
|
|
for (int k = 0; k < numneigh[lammps_j]; ++k) {
|
|
int lammps_k = firstneigh[lammps_j][k];
|
|
// first atom location
|
|
xb.copy(xPointer_[lammps_k],3);
|
|
double * xk = xatom[lammps_k];
|
|
// difference vector
|
|
xba = xb - xa;
|
|
xbI = xba + xaI;
|
|
kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
|
|
// compute virial
|
|
if (lam1 < lam2) {
|
|
bond_value
|
|
= kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
|
|
double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
|
|
double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
|
|
double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
|
|
double rsq = delx*delx + dely*dely + delz*delz;
|
|
double fforce = 0;
|
|
double pair_pe = lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
|
|
fforce *= 0.5; // dbl count
|
|
fforce *= (delx*uVariationVelocity_(j,0) +
|
|
dely*uVariationVelocity_(j,1) +
|
|
delz*uVariationVelocity_(j,2));
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
|
|
double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
|
|
double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
|
|
flux(inode,0) +=fforce*delX*bond_value;
|
|
flux(inode,1) +=fforce*delY*bond_value;
|
|
flux(inode,2) +=fforce*delZ*bond_value;
|
|
}
|
|
else { //if (atomToElementMapType_ == EULERIAN) {
|
|
flux(inode,0) +=fforce*delx*bond_value;
|
|
flux(inode,1) +=fforce*dely*bond_value;
|
|
flux(inode,2) +=fforce*delz*bond_value;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
// mesh-based kernel functions
|
|
int nodes_per_element = feEngine_->get_feMesh()->get_nNodesPerElement();
|
|
Array<int> node_list(nodes_per_element);
|
|
DENS_VEC shp(nodes_per_element);
|
|
DENS_VEC xa(nsd_),xb(nsd_),xab(nsd_),xlambda(nsd_);
|
|
for (int j = 0; j < nLocal_; j++) {
|
|
// first atom location
|
|
int lammps_j = internalToAtom_(j);
|
|
xa.copy(xPointer_[lammps_j],3);
|
|
for (int k = 0; k < numneigh[lammps_j]; ++k) {
|
|
int lammps_k = firstneigh[lammps_j][k];
|
|
//if (lammps_k < lammps_j) continue; // full neighbor list
|
|
// second (neighbor) atom location
|
|
xb.copy(xPointer_[lammps_k],3);
|
|
double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
|
|
double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
|
|
double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
|
|
double rsq = delx*delx + dely*dely + delz*delz;
|
|
double fforce = 0;
|
|
double pair_pe = lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
|
|
fforce *= 0.5; // 1/2 sum_ab = sum_(ab)
|
|
fforce *= (delx*uVariationVelocity_(j,0) +
|
|
dely*uVariationVelocity_(j,1) +
|
|
delz*uVariationVelocity_(j,2));
|
|
double flux_vec[3];
|
|
if (atomToElementMapType_ == LAGRANGIAN) {
|
|
double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
|
|
double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
|
|
double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
|
|
flux_vec[0] =fforce*delX;
|
|
flux_vec[1] =fforce*delY;
|
|
flux_vec[2] =fforce*delZ;
|
|
}
|
|
else {//if (atomToElementMapType_ == EULERIAN) {
|
|
flux_vec[0] =fforce*delx;
|
|
flux_vec[1] =fforce*dely;
|
|
flux_vec[2] =fforce*delz;
|
|
}
|
|
lam1 = 0.0; lam2 = 1.0;
|
|
double del_lambda = 0.5*(lam2 - lam1);
|
|
double avg_lambda = 0.5*(lam2 + lam1);
|
|
xab = xa - xb;
|
|
for (int i = 0; i < line_ngauss; i++) {
|
|
double lambda = del_lambda*line_xg[i] +avg_lambda;
|
|
xlambda = lambda*xab + xb;
|
|
// NOTE doesn't work for the case that fe_region \subset md_region
|
|
int dummyEltID;
|
|
feEngine_->shape_functions(xlambda,shp,dummyEltID,
|
|
node_list,latticePeriodicity);
|
|
// accumulate to nodes whose support overlaps the integration point
|
|
for (int I = 0; I < nodes_per_element; I++) {
|
|
// Use factor of 0.5 to account for line integration
|
|
// domain of 0 to 1 instead of -1 to 1
|
|
int inode = node_list(I);
|
|
double bond_value
|
|
= invNodeVolumes_(inode,inode)*shp(I)*line_wg[i]*0.5;
|
|
flux(inode,0) += flux_vec[0]*bond_value;
|
|
flux(inode,1) += flux_vec[1]*bond_value;
|
|
flux(inode,2) += flux_vec[2]*bond_value;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}// end if kernelFunction_
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// correction of node-pair distances due to periodicity
|
|
void ATC_TransferHardy::periodicity_correction(DENS_VEC& xaI)
|
|
{
|
|
for (int m = 0; m < nsd_; m++) {
|
|
xaI(m) -= periodicity[m]*box_length[m]*rnd(xaI(m)/box_length[m]);
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// set xPointer_ to xref or xatom depending on Lagrangian/Eulerian analysis
|
|
void ATC_TransferHardy::set_xPointer()
|
|
{
|
|
xPointer_ = xref_;
|
|
if (atomToElementMapType_ == EULERIAN) {
|
|
xPointer_ = lammpsInterface_->xatom();
|
|
}
|
|
}
|
|
|
|
//-------------------------------------------------------------------
|
|
// check consistency of fieldFlags_
|
|
void ATC_TransferHardy::check_fieldFlags_consistency()
|
|
{
|
|
if (fieldFlags_(HARDY_TRANSFORMED_STRESS)) {
|
|
fieldFlags_(HARDY_STRESS) = true;
|
|
fieldFlags_(HARDY_DISPLACEMENT) = true;
|
|
}
|
|
if (fieldFlags_(HARDY_ESHELBY_STRESS)) {
|
|
fieldFlags_(HARDY_STRESS) = true;
|
|
fieldFlags_(HARDY_ENERGY) = true;
|
|
fieldFlags_(HARDY_DISPLACEMENT) = true;
|
|
}
|
|
if (fieldFlags_(HARDY_CAUCHY_BORN_STRESS)) {
|
|
if (! (cauchyBornStress_) ) {
|
|
throw ATC_Error(0,"can't compute cauchy-born stress w/o cauchy born model");
|
|
}
|
|
fieldFlags_(HARDY_DISPLACEMENT) = true;
|
|
}
|
|
if (fieldFlags_(HARDY_VACANCY_CONCENTRATION)) {
|
|
fieldFlags_(HARDY_DISPLACEMENT) = true;
|
|
fieldFlags_(HARDY_NUMBER_DENSITY) = true;
|
|
}
|
|
if (fieldFlags_(HARDY_TEMPERATURE) || fieldFlags_(HARDY_HEAT_FLUX) ||
|
|
fieldFlags_(HARDY_ENERGY) || (fieldFlags_(HARDY_STRESS) &&
|
|
atomToElementMapType_ == EULERIAN) ) {
|
|
fieldFlags_(HARDY_VELOCITY) = true;
|
|
}
|
|
if (fieldFlags_(HARDY_VELOCITY)) {
|
|
fieldFlags_(HARDY_DENSITY) = true;
|
|
fieldFlags_(HARDY_MOMENTUM) = true;
|
|
}
|
|
if (fieldFlags_(HARDY_DISPLACEMENT)) {
|
|
fieldFlags_(HARDY_DENSITY) = true;
|
|
}
|
|
if (fieldFlags_(HARDY_TEMPERATURE) ||
|
|
fieldFlags_(HARDY_TYPE_CONCENTRATION)) {
|
|
fieldFlags_(HARDY_NUMBER_DENSITY) = true;
|
|
}
|
|
|
|
}
|
|
|
|
} // end namespace ATC
|