forked from lijiext/lammps
454 lines
17 KiB
C++
454 lines
17 KiB
C++
// ATC headers
|
|
#include "ATC_CouplingMomentum.h"
|
|
#include "ATC_Error.h"
|
|
#include "LammpsInterface.h"
|
|
#include "PrescribedDataManager.h"
|
|
#include "PerAtomQuantity.h"
|
|
#include "TransferOperator.h"
|
|
|
|
// Other Headers
|
|
#include <vector>
|
|
#include <map>
|
|
#include <set>
|
|
#include <utility>
|
|
#include <iostream>
|
|
|
|
using std::string;
|
|
|
|
namespace ATC {
|
|
|
|
//--------------------------------------------------------
|
|
//--------------------------------------------------------
|
|
// Class ATC_CouplingMomentum
|
|
//--------------------------------------------------------
|
|
//--------------------------------------------------------
|
|
|
|
//--------------------------------------------------------
|
|
// Constructor
|
|
//--------------------------------------------------------
|
|
ATC_CouplingMomentum::ATC_CouplingMomentum(string groupName,
|
|
double **& perAtomArray,
|
|
LAMMPS_NS::Fix * thisFix,
|
|
string matParamFile,
|
|
PhysicsType intrinsicModel,
|
|
ExtrinsicModelType extrinsicModel)
|
|
: ATC_Coupling(groupName,perAtomArray,thisFix),
|
|
refPE_(0)
|
|
{
|
|
// Allocate PhysicsModel
|
|
create_physics_model(intrinsicModel, matParamFile);
|
|
|
|
// create extrinsic physics model
|
|
if (extrinsicModel != NO_MODEL) {
|
|
extrinsicModelManager_.create_model(extrinsicModel,matParamFile);
|
|
}
|
|
|
|
// set up field data based on physicsModel
|
|
physicsModel_->num_fields(fieldSizes_,fieldMask_);
|
|
|
|
// Defaults
|
|
set_time();
|
|
bndyIntType_ = FE_INTERPOLATION;
|
|
trackCharge_ = false;
|
|
|
|
// use a kinetostat
|
|
atomicRegulator_ = new Kinetostat(this);
|
|
|
|
// set time integrator and change any defaults based on model type
|
|
if (intrinsicModel == ELASTIC) {
|
|
trackDisplacement_ = true;
|
|
fieldSizes_[DISPLACEMENT] = fieldSizes_[VELOCITY];
|
|
timeIntegrators_[VELOCITY] = new MomentumTimeIntegrator(this,TimeIntegrator::VERLET);
|
|
ghostManager_.set_boundary_dynamics(GhostManager::PRESCRIBED);
|
|
}
|
|
else if (intrinsicModel == SHEAR) {
|
|
atomToElementMapType_ = EULERIAN;
|
|
atomToElementMapFrequency_ = 1;
|
|
timeIntegrators_[VELOCITY] = new MomentumTimeIntegrator(this,TimeIntegrator::GEAR);
|
|
ghostManager_.set_boundary_dynamics(GhostManager::NO_BOUNDARY_DYNAMICS);
|
|
}
|
|
|
|
// output variable vector info:
|
|
// output[1] = total coarse scale kinetic energy
|
|
// output[2] = total coarse scale potential energy
|
|
// output[3] = total coarse scale energy
|
|
scalarFlag_ = 1;
|
|
vectorFlag_ = 1;
|
|
sizeVector_ = 5;
|
|
scalarVectorFreq_ = 1;
|
|
extVector_ = 1;
|
|
thermoEnergyFlag_ = 1;
|
|
if (extrinsicModel != NO_MODEL)
|
|
sizeVector_ += extrinsicModelManager_.size_vector(sizeVector_);
|
|
}
|
|
|
|
//--------------------------------------------------------
|
|
// Destructor
|
|
//--------------------------------------------------------
|
|
ATC_CouplingMomentum::~ATC_CouplingMomentum()
|
|
{
|
|
interscaleManager_.clear();
|
|
}
|
|
|
|
//--------------------------------------------------------
|
|
// initialize
|
|
// sets up all the necessary data
|
|
//--------------------------------------------------------
|
|
void ATC_CouplingMomentum::initialize()
|
|
{
|
|
// clear displacement entries if requested
|
|
if (!trackDisplacement_) {
|
|
fieldSizes_.erase(DISPLACEMENT);
|
|
for (int i = 0; i < NUM_FLUX; i++)
|
|
fieldMask_(DISPLACEMENT,i) = false;
|
|
}
|
|
|
|
// Base class initalizations
|
|
ATC_Coupling::initialize();
|
|
|
|
// reset integration field mask
|
|
intrinsicMask_.reset(NUM_FIELDS,NUM_FLUX);
|
|
intrinsicMask_ = false;
|
|
for (int i = 0; i < NUM_FLUX; i++)
|
|
intrinsicMask_(VELOCITY,i) = fieldMask_(VELOCITY,i);
|
|
|
|
|
|
refPE_=0;
|
|
refPE_=potential_energy();
|
|
}
|
|
|
|
//--------------------------------------------------------
|
|
// construct_transfers
|
|
// constructs needed transfer operators
|
|
//--------------------------------------------------------
|
|
void ATC_CouplingMomentum::construct_transfers()
|
|
{
|
|
ATC_Coupling::construct_transfers();
|
|
|
|
// momentum of each atom
|
|
AtomicMomentum * atomicMomentum = new AtomicMomentum(this);
|
|
interscaleManager_.add_per_atom_quantity(atomicMomentum,
|
|
"AtomicMomentum");
|
|
|
|
// nodal momentum for RHS
|
|
AtfShapeFunctionRestriction * nodalAtomicMomentum = new AtfShapeFunctionRestriction(this,
|
|
atomicMomentum,
|
|
shpFcn_);
|
|
interscaleManager_.add_dense_matrix(nodalAtomicMomentum,
|
|
"NodalAtomicMomentum");
|
|
|
|
// nodal forces
|
|
FundamentalAtomQuantity * atomicForce = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE);
|
|
AtfShapeFunctionRestriction * nodalAtomicForce = new AtfShapeFunctionRestriction(this,
|
|
atomicForce,
|
|
shpFcn_);
|
|
interscaleManager_.add_dense_matrix(nodalAtomicForce,
|
|
"NodalAtomicForce");
|
|
|
|
// nodal velocity derived only from atoms
|
|
AtfShapeFunctionMdProjection * nodalAtomicVelocity = new AtfShapeFunctionMdProjection(this,
|
|
nodalAtomicMomentum,
|
|
VELOCITY);
|
|
interscaleManager_.add_dense_matrix(nodalAtomicVelocity,
|
|
"NodalAtomicVelocity");
|
|
|
|
if (trackDisplacement_) {
|
|
// mass-weighted (center-of-mass) displacement of each atom
|
|
AtomicMassWeightedDisplacement * atomicMassWeightedDisplacement;
|
|
if (needXrefProcessorGhosts_ || groupbitGhost_) { // explicit construction on internal group
|
|
PerAtomQuantity<double> * atomReferencePositions = interscaleManager_.per_atom_quantity("AtomicInternalReferencePositions");
|
|
atomicMassWeightedDisplacement = new AtomicMassWeightedDisplacement(this,atomPositions_,
|
|
atomMasses_,
|
|
atomReferencePositions,
|
|
INTERNAL);
|
|
}
|
|
else
|
|
atomicMassWeightedDisplacement = new AtomicMassWeightedDisplacement(this);
|
|
interscaleManager_.add_per_atom_quantity(atomicMassWeightedDisplacement,
|
|
"AtomicMassWeightedDisplacement");
|
|
|
|
// nodal (RHS) mass-weighted displacement
|
|
AtfShapeFunctionRestriction * nodalAtomicMassWeightedDisplacement = new AtfShapeFunctionRestriction(this,
|
|
atomicMassWeightedDisplacement,
|
|
shpFcn_);
|
|
interscaleManager_.add_dense_matrix(nodalAtomicMassWeightedDisplacement,
|
|
"NodalAtomicMassWeightedDisplacement");
|
|
|
|
// nodal displacement derived only from atoms
|
|
AtfShapeFunctionMdProjection * nodalAtomicDisplacement = new AtfShapeFunctionMdProjection(this,
|
|
nodalAtomicMassWeightedDisplacement,
|
|
VELOCITY);
|
|
interscaleManager_.add_dense_matrix(nodalAtomicDisplacement,
|
|
"NodalAtomicDisplacement");
|
|
}
|
|
|
|
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
|
|
(_tiIt_->second)->construct_transfers();
|
|
}
|
|
atomicRegulator_->construct_transfers();
|
|
}
|
|
|
|
//---------------------------------------------------------
|
|
// init_filter
|
|
// sets up the time filtering operations in all objects
|
|
//---------------------------------------------------------
|
|
void ATC_CouplingMomentum::init_filter()
|
|
{
|
|
|
|
ATC_Coupling::init_filter();
|
|
|
|
if (timeFilterManager_.end_equilibrate() && equilibriumStart_) // set up correct initial lambda forces to enforce initial accerlation
|
|
if (atomicRegulator_->coupling_mode()==AtomicRegulator::FLUX || atomicRegulator_->coupling_mode()==AtomicRegulator::GHOST_FLUX)
|
|
// nothing needed in other cases since kinetostat force is balanced by boundary flux in FE equations
|
|
atomicRegulator_->reset_lambda_contribution(nodalAtomicFieldsRoc_[VELOCITY].quantity());
|
|
}
|
|
|
|
//--------------------------------------------------------
|
|
// modify
|
|
// parses inputs and modifies state of the filter
|
|
//--------------------------------------------------------
|
|
bool ATC_CouplingMomentum::modify(int narg, char **arg)
|
|
{
|
|
bool foundMatch = false;
|
|
int argIndex = 0;
|
|
|
|
// check to see if it is a transfer class command
|
|
// check derived class before base class
|
|
|
|
// pass-through to kinetostat
|
|
if (strcmp(arg[argIndex],"control")==0) {
|
|
argIndex++;
|
|
foundMatch = atomicRegulator_->modify(narg-argIndex,&arg[argIndex]);
|
|
}
|
|
|
|
// pass-through to timeIntegrator class
|
|
else if (strcmp(arg[argIndex],"time_integration")==0) {
|
|
argIndex++;
|
|
foundMatch = timeIntegrators_[VELOCITY]->modify(narg-argIndex,&arg[argIndex]);
|
|
}
|
|
|
|
// switch for if displacement is tracked or not
|
|
/*! \page man_track_displacement fix_modify AtC track_displacement
|
|
\section syntax
|
|
fix_modify AtC track_displacement <on/off> \n
|
|
\section examples
|
|
<TT> fix_modify atc track_displacement on </TT> \n
|
|
\section description
|
|
Determines whether displacement is tracked or not. For solids problems this is a useful quantity, but for fluids it is not relevant.
|
|
\section restrictions
|
|
Some constitutive models require the displacement field
|
|
\section default
|
|
on
|
|
*/
|
|
else if (strcmp(arg[argIndex],"track_displacement")==0) {
|
|
argIndex++;
|
|
if (strcmp(arg[argIndex],"on")==0) {
|
|
trackDisplacement_ = true;
|
|
foundMatch = true;
|
|
}
|
|
else if (strcmp(arg[argIndex],"off")==0) {
|
|
trackDisplacement_ = false;
|
|
foundMatch = true;
|
|
}
|
|
if (foundMatch) {
|
|
needReset_ = true;
|
|
}
|
|
}
|
|
|
|
else if (strcmp(arg[argIndex],"boundary_dynamics")==0) {
|
|
argIndex++;
|
|
foundMatch = ghostManager_.modify(narg-argIndex,&arg[argIndex]);
|
|
}
|
|
|
|
// no match, call base class parser
|
|
if (!foundMatch) {
|
|
foundMatch = ATC_Coupling::modify(narg, arg);
|
|
}
|
|
|
|
return foundMatch;
|
|
|
|
}
|
|
|
|
//--------------------------------------------------------
|
|
// min_pre_force
|
|
// add to interatomic forces for minimize
|
|
//--------------------------------------------------------
|
|
void ATC_CouplingMomentum::min_pre_force()
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------
|
|
// min_post_force
|
|
// add to interatomic forces for minimize
|
|
// this determines the search direction
|
|
//--------------------------------------------------------
|
|
void ATC_CouplingMomentum::min_post_force()
|
|
{
|
|
// reset positions and shape functions
|
|
ATC_Method::min_post_force();
|
|
|
|
// Set sources
|
|
|
|
prescribedDataMgr_->set_sources(time(),sources_);
|
|
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
|
|
extrinsicModelManager_.pre_final_integrate();
|
|
|
|
|
|
|
|
|
|
if (outputNow_) {
|
|
update_time(1.0);
|
|
update_step();
|
|
output();
|
|
outputNow_ = false;
|
|
}
|
|
|
|
|
|
localStep_ += 1;
|
|
}
|
|
|
|
//--------------------------------------------------------
|
|
// output
|
|
// does post-processing steps and outputs data
|
|
//--------------------------------------------------------
|
|
void ATC_CouplingMomentum::output()
|
|
{
|
|
if (output_now()) {
|
|
feEngine_->departition_mesh();
|
|
OUTPUT_LIST outputData;
|
|
|
|
// base class output
|
|
ATC_Method::output();
|
|
|
|
// push atc fields time integrator modifies into output arrays
|
|
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
|
|
(_tiIt_->second)->post_process();
|
|
}
|
|
|
|
// auxilliary data
|
|
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
|
|
(_tiIt_->second)->output(outputData);
|
|
}
|
|
atomicRegulator_->output(outputData);
|
|
extrinsicModelManager_.output(outputData);
|
|
|
|
DENS_MAT & velocity(nodalAtomicFields_[VELOCITY].set_quantity());
|
|
DENS_MAT & rhs(rhs_[VELOCITY].set_quantity());
|
|
if (lammpsInterface_->rank_zero()) {
|
|
// mesh data
|
|
outputData["NodalAtomicVelocity"] = &velocity;
|
|
outputData["FE_Force"] = &rhs;
|
|
if (trackDisplacement_) {
|
|
outputData["NodalAtomicDisplacement"] = & nodalAtomicFields_[DISPLACEMENT].set_quantity();
|
|
}
|
|
|
|
feEngine_->write_data(output_index(), fields_, & outputData);
|
|
}
|
|
// force optional variables to reset to keep in sync
|
|
if (trackDisplacement_) {
|
|
nodalAtomicFields_[DISPLACEMENT].force_reset();
|
|
}
|
|
|
|
feEngine_->partition_mesh();
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------
|
|
// compute_scalar : added energy
|
|
// this is used in the line search
|
|
//--------------------------------------------------------------------
|
|
double ATC_CouplingMomentum::compute_scalar(void)
|
|
{
|
|
double energy = extrinsicModelManager_.compute_scalar();
|
|
return energy;
|
|
}
|
|
|
|
//--------------------------------------------------------------------
|
|
// kinetic energy
|
|
//--------------------------------------------------------------------
|
|
double ATC_CouplingMomentum::kinetic_energy(const IntegrationDomainType domain) // const
|
|
{
|
|
const MATRIX & M = massMats_[VELOCITY].quantity();
|
|
const DENS_MAT & velocity(fields_[VELOCITY].quantity());
|
|
double kineticEnergy = 0;
|
|
for (int j = 0; j < nsd_; j++) {
|
|
CLON_VEC v = column(velocity,j);
|
|
kineticEnergy += v.dot(M*v);
|
|
}
|
|
if (domain == FE_DOMAIN) {
|
|
|
|
Array<FieldName> massMask(1);
|
|
massMask(0) = VELOCITY;
|
|
feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel_,atomMaterialGroups_,
|
|
atomVolume_->quantity(),shpFcn_->quantity(),
|
|
Ma_);
|
|
const MATRIX & Ma = Ma_[VELOCITY].quantity();
|
|
for (int j = 0; j < nsd_; j++) {
|
|
CLON_VEC v = column(velocity,j);
|
|
kineticEnergy -= v.dot(Ma*v);
|
|
}
|
|
}
|
|
double mvv2e = lammpsInterface_->mvv2e();
|
|
kineticEnergy *= 0.5*mvv2e; // convert to LAMMPS units
|
|
|
|
return kineticEnergy;
|
|
}
|
|
//--------------------------------------------------------------------
|
|
// potential/strain energy
|
|
//--------------------------------------------------------------------
|
|
double ATC_CouplingMomentum::potential_energy(const IntegrationDomainType domain) const
|
|
{
|
|
Array<FieldName> mask(1);
|
|
mask(0) = VELOCITY;
|
|
FIELD_MATS energy;
|
|
feEngine_->compute_energy(mask,
|
|
fields_,
|
|
physicsModel_,
|
|
elementToMaterialMap_,
|
|
energy,
|
|
&(elementMask_->quantity()),
|
|
domain);
|
|
double potentialEnergy = energy[VELOCITY].col_sum();
|
|
double mvv2e = lammpsInterface_->mvv2e();
|
|
potentialEnergy *= mvv2e; // convert to LAMMPS units
|
|
return potentialEnergy-refPE_;
|
|
}
|
|
//--------------------------------------------------------------------
|
|
// compute_vector
|
|
//--------------------------------------------------------------------
|
|
// this is for direct output to lammps thermo
|
|
double ATC_CouplingMomentum::compute_vector(int n)
|
|
{
|
|
// output[1] = total coarse scale kinetic energy
|
|
// output[2] = total coarse scale potential energy
|
|
// output[3] = total coarse scale energy
|
|
// output[4] = fe-only coarse scale kinetic energy
|
|
// output[5] = fe-only coarse scale potential energy
|
|
|
|
|
|
|
|
if (n == 0) {
|
|
return kinetic_energy();
|
|
}
|
|
else if (n == 1) {
|
|
return potential_energy();
|
|
}
|
|
else if (n == 2) {
|
|
return kinetic_energy()+potential_energy();
|
|
}
|
|
else if (n == 3) {
|
|
return kinetic_energy(FE_DOMAIN);
|
|
}
|
|
else if (n == 4) {
|
|
return potential_energy(FE_DOMAIN);
|
|
}
|
|
else if (n > 4) {
|
|
double extrinsicValue = extrinsicModelManager_.compute_vector(n);
|
|
return extrinsicValue;
|
|
}
|
|
return 0.;
|
|
|
|
}
|
|
|
|
};
|