Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F120132688
ATC_CouplingMomentum.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jul 2, 04:10
Size
16 KB
Mime Type
text/x-c
Expires
Fri, Jul 4, 04:10 (2 d)
Engine
blob
Format
Raw Data
Handle
26730241
Attached To
rLAMMPS lammps
ATC_CouplingMomentum.cpp
View Options
// 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.;
}
};
Event Timeline
Log In to Comment