Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73835112
Kinetostat.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 24, 18:27
Size
86 KB
Mime Type
text/x-c
Expires
Fri, Jul 26, 18:27 (2 d)
Engine
blob
Format
Raw Data
Handle
19278025
Attached To
rLAMMPS lammps
Kinetostat.cpp
View Options
#include "Kinetostat.h"
#include "ATC_Error.h"
#include "ATC_Coupling.h"
#include "LammpsInterface.h"
#include "PerAtomQuantityLibrary.h"
#include "PrescribedDataManager.h"
#include "ElasticTimeIntegrator.h"
#include "TransferOperator.h"
namespace ATC {
//--------------------------------------------------------
//--------------------------------------------------------
// Class Kinetostat
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
Kinetostat::Kinetostat(ATC_Coupling * atc,
const string & regulatorPrefix) :
AtomicRegulator(atc,regulatorPrefix)
{
// do nothing
}
//--------------------------------------------------------
// modify:
// parses and adjusts kinetostat state based on
// user input, in the style of LAMMPS user input
//--------------------------------------------------------
bool Kinetostat::modify(int narg, char **arg)
{
bool foundMatch = false;
int argIndex = 0;
if (strcmp(arg[argIndex],"momentum")==0) {
argIndex++;
// fluxstat type
/*! \page man_control_momentum fix_modify AtC control momentum
\section syntax
fix_modify AtC control momentum none \n
fix_modify AtC control momentum rescale <frequency>\n
- frequency (int) = time step frequency for applying displacement and velocity rescaling \n
fix_modify AtC control momentum glc_displacement \n
fix_modify AtC control momentum glc_velocity \n
fix_modify AtC control momentum hoover \n
fix_modify AtC control momentum flux [faceset face_set_id, interpolate]
- face_set_id (string) = id of boundary face set, if not specified
(or not possible when the atomic domain does not line up with
mesh boundaries) defaults to an atomic-quadrature approximate
evaulation\n
\section examples
fix_modify AtC control momentum glc_velocity \n
fix_modify AtC control momentum stress_flux faceset bndy_faces \n
\section description
\section restrictions
only for be used with specific transfers :
elastic \n
rescale not valid with time filtering activated
\section related
\section default
none
*/
boundaryIntegrationType_ = NO_QUADRATURE;
howOften_ = 1;
if (strcmp(arg[argIndex],"none")==0) { // restore defaults
regulatorTarget_ = NONE;
couplingMode_ = UNCOUPLED;
foundMatch = true;
}
else if (strcmp(arg[argIndex],"glc_displacement")==0) {
regulatorTarget_ = FIELD;
couplingMode_ = FIXED;
foundMatch = true;
}
else if (strcmp(arg[argIndex],"glc_velocity")==0) {
regulatorTarget_ = DERIVATIVE;
couplingMode_ = FIXED;
foundMatch = true;
}
else if (strcmp(arg[argIndex],"hoover")==0) {
regulatorTarget_ = DYNAMICS;
couplingMode_ = FIXED;
foundMatch = true;
}
else if (strcmp(arg[argIndex],"flux")==0) {
regulatorTarget_ = DYNAMICS;
couplingMode_ = FLUX;
argIndex++;
boundaryIntegrationType_ = atc_->parse_boundary_integration(narg-argIndex,&arg[argIndex],boundaryFaceSet_);
foundMatch = true;
}
else if (strcmp(arg[argIndex],"ghost_flux")==0) {
regulatorTarget_ = DYNAMICS;
couplingMode_ = GHOST_FLUX;
foundMatch = true;
}
}
if (!foundMatch)
foundMatch = AtomicRegulator::modify(narg,arg);
if (foundMatch)
needReset_ = true;
return foundMatch;
}
//--------------------------------------------------------
// reset_lambda_contribution
// resets the kinetostat generated force to a
// prescribed value
//--------------------------------------------------------
void Kinetostat::reset_lambda_contribution(const DENS_MAT & target)
{
DENS_MAN * lambdaForceFiltered = regulator_data("LambdaForceFiltered",nsd_);
lambdaForceFiltered->set_quantity() = target;
}
//--------------------------------------------------------
// initialize:
// sets up methods before a run
// dependence, but in general there is also a
// time integrator dependence. In general the
// precedence order is:
// time filter -> time integrator -> kinetostat
// In the future this may need to be added if
// different types of time integrators can be
// specified.
//--------------------------------------------------------
void Kinetostat::construct_methods()
{
// get data associated with stages 1 & 2 of ATC_Method::initialize
AtomicRegulator::construct_methods();
if (atc_->reset_methods()) {
// eliminate existing methods
delete_method();
DENS_MAT nodalGhostForceFiltered;
TimeIntegrator::TimeIntegrationType myIntegrationType = (atc_->time_integrator(VELOCITY))->time_integration_type();
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (timeFilterManager->end_equilibrate() && regulatorTarget_==AtomicRegulator::DYNAMICS) {
StressFlux * myMethod;
myMethod = dynamic_cast<StressFlux *>(regulatorMethod_);
nodalGhostForceFiltered = (myMethod->filtered_ghost_force()).quantity();
}
// update time filter
if (timeFilterManager->need_reset()) {
timeFilter_ = timeFilterManager->construct(TimeFilterManager::IMPLICIT_UPDATE);
}
if (timeFilterManager->filter_dynamics()) {
switch (regulatorTarget_) {
case NONE: {
regulatorMethod_ = new RegulatorMethod(this);
break;
}
case FIELD: {
regulatorMethod_ = new DisplacementGlcFiltered(this);
break;
}
case DERIVATIVE: {
regulatorMethod_ = new VelocityGlcFiltered(this);
break;
}
case DYNAMICS: {
throw ATC_Error("Kinetostat::initialize - force based kinetostats not yet implemented with time filtering");
regulatorMethod_ = new StressFluxFiltered(this);
if (timeFilterManager->end_equilibrate()) {
StressFlux * myMethod;
myMethod = dynamic_cast<StressFlux *>(regulatorMethod_);
myMethod->reset_filtered_ghost_force(nodalGhostForceFiltered);
}
break;
}
default:
throw ATC_Error("Unknown kinetostat type in Kinetostat::initialize");
}
}
else {
switch (regulatorTarget_) {
case NONE: {
regulatorMethod_ = new RegulatorMethod(this);
break;
}
case FIELD: {
regulatorMethod_ = new DisplacementGlc(this);
break;
}
case DERIVATIVE: {
regulatorMethod_ = new VelocityGlc(this);
break;
}
case DYNAMICS: {
if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
if (couplingMode_ == GHOST_FLUX) {
regulatorMethod_ = new KinetostatFluxGhost(this);
}
else if (couplingMode_ == FIXED) {
regulatorMethod_ = new KinetostatFixed(this);
}
else if (couplingMode_ == FLUX) {
regulatorMethod_ = new KinetostatFlux(this);
}
break;
}
if (myIntegrationType == TimeIntegrator::GEAR) {
if (couplingMode_ == FIXED) {
regulatorMethod_ = new KinetostatFixed(this);
}
else if (couplingMode_ == FLUX) {
regulatorMethod_ = new KinetostatFlux(this);
}
break;
}
else {
if (couplingMode_ == GHOST_FLUX) {
regulatorMethod_ = new StressFluxGhost(this);
}
else {
regulatorMethod_ = new StressFlux(this);
}
break;
}
}
default:
throw ATC_Error("Unknown kinetostat type in Kinetostat::initialize");
}
AtomicRegulator::reset_method();
}
}
else {
set_all_data_to_used();
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class KinetostatShapeFunction
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
KinetostatShapeFunction::KinetostatShapeFunction(Kinetostat *kinetostat,
const string & regulatorPrefix) :
RegulatorShapeFunction(kinetostat,regulatorPrefix),
kinetostat_(kinetostat),
timeFilter_(atomicRegulator_->time_filter()),
nodalAtomicLambdaForce_(NULL),
lambdaForceFiltered_(NULL),
atomKinetostatForce_(NULL),
atomVelocities_(NULL),
atomMasses_(NULL)
{
// data associated with stage 3 in ATC_Method::initialize
lambda_ = kinetostat->regulator_data(regulatorPrefix_+"LambdaMomentum",nsd_);
lambdaForceFiltered_ = kinetostat_->regulator_data("LambdaForceFiltered",nsd_);
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void KinetostatShapeFunction::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// needed fundamental quantities
atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
// base class transfers
RegulatorShapeFunction::construct_transfers();
// lambda interpolated to the atomic coordinates
atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
lambda_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_per_atom_quantity(atomLambdas_,
regulatorPrefix_+"AtomLambdaMomentum");
}
//--------------------------------------------------------
// set_weights
// sets diagonal weighting matrix used in
// solve_for_lambda
//--------------------------------------------------------
void KinetostatShapeFunction::set_weights()
{
if (this->use_local_shape_functions()) {
ConstantQuantityMapped<double> * myWeights = new ConstantQuantityMapped<double>(atc_,1.,lambdaAtomMap_);
weights_ = myWeights;
(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
"AtomOnesMapped");
}
else {
weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicOnes");
if (!weights_) {
weights_ = new ConstantQuantity<double>(atc_,1.);
(atc_->interscale_manager()).add_per_atom_quantity(weights_,
"AtomicOnes");
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class GlcKinetostat
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
GlcKinetostat::GlcKinetostat(Kinetostat *kinetostat) :
KinetostatShapeFunction(kinetostat),
mdMassMatrix_(atc_->set_mass_mat_md(VELOCITY)),
atomPositions_(NULL)
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void GlcKinetostat::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// needed fundamental quantities
atomPositions_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION);
// base class transfers
KinetostatShapeFunction::construct_transfers();
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void GlcKinetostat::initialize()
{
KinetostatShapeFunction::initialize();
// set up list of nodes using Hoover coupling
// (a) nodes with prescribed values
PrescribedDataManager * prescribedDataMgr(atc_->prescribed_data_manager());
for (int i = 0; i < nNodes_; ++i)
for (int j = 0; j < nsd_; ++j)
if (prescribedDataMgr->is_fixed(i,VELOCITY,j))
hooverNodes_.insert(pair<int,int>(i,j));
// (b) AtC coupling nodes
if (atomicRegulator_->coupling_mode()==AtomicRegulator::FIXED) {
InterscaleManager & interscaleManager(atc_->interscale_manager());
const INT_ARRAY & nodeType((interscaleManager.dense_matrix_int("NodalGeometryType"))->quantity());
if (atomicRegulator_->use_localized_lambda()) {
for (int i = 0; i < nNodes_; ++i) {
if (nodeType(i,0)==BOUNDARY) {
for (int j = 0; j < nsd_; ++j) {
hooverNodes_.insert(pair<int,int>(i,j));
}
}
}
}
else {
for (int i = 0; i < nNodes_; ++i) {
if (nodeType(i,0)==BOUNDARY || nodeType(i,0)==MD_ONLY) {
for (int j = 0; j < nsd_; ++j) {
hooverNodes_.insert(pair<int,int>(i,j));
}
}
}
}
}
}
//--------------------------------------------------------
// apply_lambda_to_atoms
// uses existing lambda to modify given
// atomic quantity
//--------------------------------------------------------
void GlcKinetostat::apply_to_atoms(PerAtomQuantity<double> * quantity,
const DENS_MAT & lambdaAtom,
double dt)
{
*quantity -= lambdaAtom;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class DisplacementGlc
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
DisplacementGlc::DisplacementGlc(Kinetostat * kinetostat) :
GlcKinetostat(kinetostat),
nodalAtomicMassWeightedDisplacement_(NULL),
nodalDisplacements_(atc_->field(DISPLACEMENT))
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void DisplacementGlc::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// set up shape function matrix
if (this->use_local_shape_functions()) {
lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
regulatorPrefix_+"LambdaAtomMap");
shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
lambdaAtomMap_,
nodeToOverlapMap_);
}
else {
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
}
interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixMomentum");
// set linear solver strategy
if (atomicRegulator_->use_lumped_lambda_solve()) {
linearSolverType_ = AtomicRegulator::RSL_SOLVE;
}
else {
linearSolverType_ = AtomicRegulator::CG_SOLVE;
}
// base class transfers
GlcKinetostat::construct_transfers();
// atomic force induced by kinetostat
atomKinetostatForce_ = new AtomicKinetostatForceDisplacement(atc_);
interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
regulatorPrefix_+"AtomKinetostatForce");
// restricted force due to kinetostat
nodalAtomicLambdaForce_ = new AtfShapeFunctionRestriction(atc_,
atomKinetostatForce_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_dense_matrix(nodalAtomicLambdaForce_,
regulatorPrefix_+"NodalAtomicLambdaForce");
// nodal displacement restricted from atoms
nodalAtomicMassWeightedDisplacement_ = interscaleManager.dense_matrix("NodalAtomicMassWeightedDisplacement");
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void DisplacementGlc::initialize()
{
GlcKinetostat::initialize();
// sets up time filter for cases where variables temporally filtered
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (!timeFilterManager->end_equilibrate()) {
*lambdaForceFiltered_ = 0.;
timeFilter_->initialize(lambdaForceFiltered_->quantity());
}
}
//--------------------------------------------------------
// apply:
// apply the kinetostat to the atoms
//--------------------------------------------------------
void DisplacementGlc::apply_post_predictor(double dt)
{
compute_kinetostat(dt);
}
//--------------------------------------------------------
// compute_kinetostat
// manages the solution and application of the
// kinetostat equations and variables
//--------------------------------------------------------
void DisplacementGlc::compute_kinetostat(double dt)
{
// initial filtering update
apply_pre_filtering(dt);
// set up rhs
DENS_MAT rhs(nNodes_,nsd_);
set_kinetostat_rhs(rhs,dt);
// solve linear system for lambda
solve_for_lambda(rhs,lambda_->set_quantity());
// compute nodal atomic power
compute_nodal_lambda_force(dt);
// apply kinetostat to atoms
apply_to_atoms(atomPositions_,atomLambdas_->quantity());
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the right-hand side of the
// kinetostat equations
//--------------------------------------------------------
void DisplacementGlc::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// form rhs : sum_a (hatN_Ia * x_ai) - (Upsilon)_Ii
rhs = nodalAtomicMassWeightedDisplacement_->quantity();
rhs -= ((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalDisplacements_.quantity());
}
//--------------------------------------------------------
// compute_nodal_lambda_force
// compute the effective FE force applied
// by the kinetostat
//--------------------------------------------------------
void DisplacementGlc::compute_nodal_lambda_force(double dt)
{
const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
timeFilter_->apply_post_step1(lambdaForceFiltered_->set_quantity(),
myNodalAtomicLambdaForce,dt);
// update FE displacements for localized thermostats
apply_localization_correction(myNodalAtomicLambdaForce,
nodalDisplacements_.set_quantity(),
dt*dt);
}
//--------------------------------------------------------
// apply_pre_filtering
// applies first step of filtering to
// relevant variables
//--------------------------------------------------------
void DisplacementGlc::apply_pre_filtering(double dt)
{
// apply time filtered lambda force
DENS_MAT lambdaZero(nNodes_,nsd_);
timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),(-1./dt/dt)*lambdaZero,dt);
}
//--------------------------------------------------------
// set_weights
// sets diagonal weighting matrix used in
// solve_for_lambda
//--------------------------------------------------------
void DisplacementGlc::set_weights()
{
if (lambdaAtomMap_) {
MappedAtomQuantity * myWeights = new MappedAtomQuantity(atc_,atomMasses_,lambdaAtomMap_);
weights_ = myWeights;
(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
"AtomMassesMapped");
}
else {
weights_ = atomMasses_;
}
}
//--------------------------------------------------------
// apply_localization_correction
// corrects for localized kinetostats only
// solving kinetostat equations on a subset
// of the MD region
//--------------------------------------------------------
void DisplacementGlc::apply_localization_correction(const DENS_MAT & source,
DENS_MAT & nodalField,
double weight)
{
DENS_MAT nodalLambdaRoc(nNodes_,nsd_);
atc_->apply_inverse_mass_matrix(source,
nodalLambdaRoc,
VELOCITY);
set<pair<int,int> >::const_iterator iter;
for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
nodalLambdaRoc(iter->first,iter->second) = 0.;
}
nodalField += weight*nodalLambdaRoc;
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void DisplacementGlc::output(OUTPUT_LIST & outputData)
{
_nodalAtomicLambdaForceOut_ = nodalAtomicLambdaForce_->quantity();
DENS_MAT & lambda(lambda_->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = λ
outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_nodalAtomicLambdaForceOut_);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class DisplacementGlcFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
DisplacementGlcFiltered::DisplacementGlcFiltered(Kinetostat * kinetostat) :
DisplacementGlc(kinetostat),
nodalAtomicDisplacements_(atc_->nodal_atomic_field(DISPLACEMENT))
{
// do nothing
}
//--------------------------------------------------------
// apply_pre_filtering
// applies first step of filtering to
// relevant variables
//--------------------------------------------------------
void DisplacementGlcFiltered::apply_pre_filtering(double dt)
{
// apply time filtered lambda to atomic fields
DisplacementGlc::apply_pre_filtering(dt);
DENS_MAT nodalAcceleration(nNodes_,nsd_);
atc_->apply_inverse_md_mass_matrix(lambdaForceFiltered_->set_quantity(),
nodalAcceleration,
VELOCITY);
nodalAtomicDisplacements_ += dt*dt*nodalAcceleration;
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the right-hand side of the
// kinetostat equations
//--------------------------------------------------------
void DisplacementGlcFiltered::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// form rhs : sum_a (hatN_Ia * x_ai) - (Upsilon)_Ii
double coef = 1./(timeFilter_->unfiltered_coefficient_pre_s1(dt));
rhs = coef*((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalAtomicDisplacements_.quantity() - nodalDisplacements_.quantity());
}
//--------------------------------------------------------
// compute_nodal_lambda_force
// compute the effective FE force applied
// by the kinetostat
//--------------------------------------------------------
void DisplacementGlcFiltered::compute_nodal_lambda_force(double dt)
{
const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
DENS_MAT & myLambdaForceFiltered(lambdaForceFiltered_->set_quantity());
timeFilter_->apply_post_step1(myLambdaForceFiltered,
myNodalAtomicLambdaForce,dt);
// update filtered atomic displacements
DENS_MAT nodalLambdaRoc(myNodalAtomicLambdaForce.nRows(),myNodalAtomicLambdaForce.nCols());
atc_->apply_inverse_md_mass_matrix(myNodalAtomicLambdaForce,
nodalLambdaRoc,
VELOCITY);
timeFilter_->apply_post_step1(nodalAtomicDisplacements_.set_quantity(),dt*dt*nodalLambdaRoc,dt);
// update FE displacements for localized thermostats
apply_localization_correction(myLambdaForceFiltered,
nodalDisplacements_.set_quantity(),
dt*dt);
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void DisplacementGlcFiltered::output(OUTPUT_LIST & outputData)
{
DENS_MAT & lambda(lambda_->set_quantity());
DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = λ
outputData[regulatorPrefix_+"NodalLambdaForce"] = &lambdaForceFiltered;
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class VelocityGlc
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
VelocityGlc::VelocityGlc(Kinetostat * kinetostat) :
GlcKinetostat(kinetostat),
nodalAtomicMomentum_(NULL),
nodalVelocities_(atc_->field(VELOCITY))
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void VelocityGlc::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// set up shape function matrix
if (this->use_local_shape_functions()) {
lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
regulatorPrefix_+"LambdaAtomMap");
shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
lambdaAtomMap_,
nodeToOverlapMap_);
}
else {
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
}
interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixMomentum");
// set linear solver strategy
if (atomicRegulator_->use_lumped_lambda_solve()) {
linearSolverType_ = AtomicRegulator::RSL_SOLVE;
}
else {
linearSolverType_ = AtomicRegulator::CG_SOLVE;
}
// base class transfers
GlcKinetostat::construct_transfers();
// atomic force induced by kinetostat
atomKinetostatForce_ = new AtomicKinetostatForceVelocity(atc_);
interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
regulatorPrefix_+"AtomKinetostatForce");
// restricted force due to kinetostat
nodalAtomicLambdaForce_ = new AtfShapeFunctionRestriction(atc_,
atomKinetostatForce_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_dense_matrix(nodalAtomicLambdaForce_,
regulatorPrefix_+"NodalAtomicLambdaForce");
// nodal momentum restricted from atoms
nodalAtomicMomentum_ = interscaleManager.dense_matrix("NodalAtomicMomentum");
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void VelocityGlc::initialize()
{
GlcKinetostat::initialize();
// sets up time filter for cases where variables temporally filtered
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (!timeFilterManager->end_equilibrate()) {
lambdaForceFiltered_->set_quantity() = 0.;
timeFilter_->initialize(lambdaForceFiltered_->quantity());
}
}
//--------------------------------------------------------
// apply_mid_corrector:
// apply the kinetostat during the middle of the
// predictor phase
//--------------------------------------------------------
void VelocityGlc::apply_mid_predictor(double dt)
{
double dtLambda = 0.5*dt;
compute_kinetostat(dtLambda);
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the kinetostat after the corrector phase
//--------------------------------------------------------
void VelocityGlc::apply_post_corrector(double dt)
{
double dtLambda = 0.5*dt;
compute_kinetostat(dtLambda);
}
//--------------------------------------------------------
// apply_pre_filtering
// applies first step of filtering to
// relevant variables
//--------------------------------------------------------
void VelocityGlc::apply_pre_filtering(double dt)
{
// apply time filtered lambda to atomic fields
DENS_MAT lambdaZero(nNodes_,nsd_);
timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),(-1./dt)*lambdaZero,dt);
}
//--------------------------------------------------------
// compute_kinetostat
// manages the solution and application of the
// kinetostat equations and variables
//--------------------------------------------------------
void VelocityGlc::compute_kinetostat(double dt)
{
// initial filtering update
apply_pre_filtering(dt);
// set up rhs
DENS_MAT rhs(nNodes_,nsd_);
set_kinetostat_rhs(rhs,dt);
// solve linear system for lambda
solve_for_lambda(rhs,lambda_->set_quantity());
// compute nodal atomic power
compute_nodal_lambda_force(dt);
// apply kinetostat to atoms
apply_to_atoms(atomVelocities_,atomLambdas_->quantity());
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the right-hand side of the
// kinetostat equations
//--------------------------------------------------------
void VelocityGlc::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// form rhs : sum_a (hatN_Ia * x_ai) - (\dot{Upsilon})_Ii
rhs = nodalAtomicMomentum_->quantity();
rhs -= ((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalVelocities_.quantity());
}
//--------------------------------------------------------
// compute_nodal_lambda_force
// compute the effective FE force applied
// by the kinetostat
//--------------------------------------------------------
void VelocityGlc::compute_nodal_lambda_force(double dt)
{
const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),
myNodalAtomicLambdaForce,dt);
// update FE displacements for localized thermostats
apply_localization_correction(myNodalAtomicLambdaForce,
nodalVelocities_.set_quantity(),
dt);
}
//--------------------------------------------------------
// set_weights
// sets diagonal weighting matrix used in
// solve_for_lambda
//--------------------------------------------------------
void VelocityGlc::set_weights()
{
if (lambdaAtomMap_) {
MappedAtomQuantity * myWeights = new MappedAtomQuantity(atc_,atomMasses_,lambdaAtomMap_);
weights_ = myWeights;
(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
"AtomMassesMapped");
}
else {
weights_ = atomMasses_;
}
}
//--------------------------------------------------------
// apply_localization_correction
// corrects for localized kinetostats only
// solving kinetostat equations on a subset
// of the MD region
//--------------------------------------------------------
void VelocityGlc::apply_localization_correction(const DENS_MAT & source,
DENS_MAT & nodalField,
double weight)
{
DENS_MAT nodalLambdaRoc(nNodes_,nsd_);
atc_->apply_inverse_mass_matrix(source,
nodalLambdaRoc,
VELOCITY);
set<pair<int,int> >::const_iterator iter;
for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
nodalLambdaRoc(iter->first,iter->second) = 0.;
}
nodalField += weight*nodalLambdaRoc;
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void VelocityGlc::output(OUTPUT_LIST & outputData)
{
_nodalAtomicLambdaForceOut_ = nodalAtomicLambdaForce_->quantity();
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = &(lambda_->set_quantity());
outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_nodalAtomicLambdaForceOut_);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class VelocityGlcFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
VelocityGlcFiltered::VelocityGlcFiltered(Kinetostat *kinetostat)
: VelocityGlc(kinetostat),
nodalAtomicVelocities_(atc_->nodal_atomic_field(VELOCITY))
{
// do nothing
}
//--------------------------------------------------------
// apply_pre_filtering
// applies first step of filtering to
// relevant variables
//--------------------------------------------------------
void VelocityGlcFiltered::apply_pre_filtering(double dt)
{
// apply time filtered lambda to atomic fields
VelocityGlc::apply_pre_filtering(dt);
DENS_MAT nodalAcceleration(nNodes_,nsd_);
atc_->apply_inverse_md_mass_matrix(lambdaForceFiltered_->quantity(),
nodalAcceleration,
VELOCITY);
nodalAtomicVelocities_ += dt*nodalAcceleration;
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the right-hand side of the
// kinetostat equations
//--------------------------------------------------------
void VelocityGlcFiltered::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// form rhs : sum_a (hatN_Ia * x_ai) - (Upsilon)_Ii
double coef = 1./(timeFilter_->unfiltered_coefficient_pre_s1(dt));
rhs = coef*((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalAtomicVelocities_.quantity() - nodalVelocities_.quantity());
}
//--------------------------------------------------------
// compute_nodal_lambda_force
// compute the effective FE force applied
// by the kinetostat
//--------------------------------------------------------
void VelocityGlcFiltered::compute_nodal_lambda_force(double dt)
{
const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
DENS_MAT & myLambdaForceFiltered(lambdaForceFiltered_->set_quantity());
timeFilter_->apply_post_step1(myLambdaForceFiltered,myNodalAtomicLambdaForce,dt);
// update filtered atomic displacements
DENS_MAT nodalLambdaRoc(myNodalAtomicLambdaForce.nRows(),myNodalAtomicLambdaForce.nCols());
atc_->apply_inverse_md_mass_matrix(myNodalAtomicLambdaForce,
nodalLambdaRoc,
VELOCITY);
timeFilter_->apply_post_step1(nodalAtomicVelocities_.set_quantity(),dt*nodalLambdaRoc,dt);
// update FE displacements for localized thermostats
apply_localization_correction(myLambdaForceFiltered,
nodalVelocities_.set_quantity(),
dt);
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void VelocityGlcFiltered::output(OUTPUT_LIST & outputData)
{
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = &(lambda_->set_quantity());
outputData[regulatorPrefix_+"NodalLambdaForce"] = &(lambdaForceFiltered_->set_quantity());
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class StressFlux
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
StressFlux::StressFlux(Kinetostat * kinetostat) :
GlcKinetostat(kinetostat),
nodalForce_(atc_->field_rhs(VELOCITY)),
nodalAtomicForce_(NULL),
nodalGhostForce_(NULL),
momentumSource_(atc_->atomic_source(VELOCITY))
{
// flag for performing boundary flux calculation
fieldMask_(VELOCITY,FLUX) = true;
}
StressFlux::~StressFlux()
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void StressFlux::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// set up shape function matrix
if (this->use_local_shape_functions()) {
lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
regulatorPrefix_+"LambdaAtomMap");
shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
lambdaAtomMap_,
nodeToOverlapMap_);
}
else {
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
}
interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixMomentum");
// set linear solver strategy
if (atomicRegulator_->use_lumped_lambda_solve()) {
linearSolverType_ = AtomicRegulator::RSL_SOLVE;
}
else {
linearSolverType_ = AtomicRegulator::CG_SOLVE;
}
// base class transfers
GlcKinetostat::construct_transfers();
// force at nodes due to atoms
nodalAtomicForce_ = interscaleManager.dense_matrix("NodalAtomicForce");
// atomic force induced by kinetostat
atomKinetostatForce_ = new AtomicKinetostatForceStress(atc_,atomLambdas_);
interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
regulatorPrefix_+"AtomKinetostatForce");
// restricted force due to kinetostat
nodalAtomicLambdaForce_ = new AtfShapeFunctionRestriction(atc_,
atomKinetostatForce_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_dense_matrix(nodalAtomicLambdaForce_,
regulatorPrefix_+"NodalAtomicLambdaForce");
// sets up space for ghost force related variables
if (atc_->groupbit_ghost()) {
GhostCouplingMatrix * shapeFunctionGhost = new GhostCouplingMatrix(atc_,interscaleManager.per_atom_sparse_matrix("InterpolantGhost"),
regulatedNodes_,nodeToOverlapMap_);
interscaleManager.add_sparse_matrix(shapeFunctionGhost,
regulatorPrefix_+"GhostCouplingMatrix");
FundamentalAtomQuantity * atomGhostForce = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE,
GHOST);
nodalGhostForce_ = new AtfShapeFunctionRestriction(atc_,atomGhostForce,
shapeFunctionGhost);
interscaleManager.add_dense_matrix(nodalGhostForce_,
regulatorPrefix_+"NodalGhostForce");
nodalGhostForceFiltered_.reset(nNodes_,nsd_);
}
}
//--------------------------------------------------------
// compute_boundary_flux:
// computes the boundary flux to be consistent with
// the controller
//--------------------------------------------------------
void StressFlux::compute_boundary_flux(FIELDS & fields)
{
GlcKinetostat::compute_boundary_flux(fields);
}
//--------------------------------------------------------
// apply_pre_predictor:
// apply the kinetostat to the atoms in the
// mid-predictor integration phase
//--------------------------------------------------------
void StressFlux::apply_mid_predictor(double dt)
{
double dtLambda = 0.5*dt;
// apply lambda force to atoms
apply_to_atoms(atomVelocities_,atomKinetostatForce_->quantity(),dtLambda);
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the kinetostat to the atoms in the
// post-corrector integration phase
//--------------------------------------------------------
void StressFlux::apply_post_corrector(double dt)
{
double dtLambda = 0.5*dt;
// apply lambda force to atoms
apply_to_atoms(atomVelocities_,atomKinetostatForce_->quantity(),dtLambda);
}
//--------------------------------------------------------
// compute_kinetostat
// manages the solution and application of the
// kinetostat equations and variables
//--------------------------------------------------------
void StressFlux::compute_kinetostat(double dt)
{
// initial filtering update
apply_pre_filtering(dt);
// set up rhs
DENS_MAT rhs(nNodes_,nsd_);
set_kinetostat_rhs(rhs,dt);
// solve linear system for lambda
solve_for_lambda(rhs,lambda_->set_quantity());
// compute nodal atomic power
compute_nodal_lambda_force(dt);
}
//--------------------------------------------------------
// apply_pre_filtering
// applies first step of filtering to
// relevant variables
//--------------------------------------------------------
void StressFlux::apply_pre_filtering(double dt)
{
// apply time filtered lambda force
DENS_MAT lambdaZero(nNodes_,nsd_);
timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),lambdaZero,dt);
if (nodalGhostForce_) {
timeFilter_->apply_pre_step1(nodalGhostForceFiltered_.set_quantity(),
nodalGhostForce_->quantity(),dt);
}
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the RHS of the kinetostat equations
// for the coupling parameter lambda
//--------------------------------------------------------
void StressFlux::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// (a) for flux based :
// form rhs : \int N_I r dV - \sum_g N_Ig^* f_g
// sources are set in ATC transfer
rhs.reset(nNodes_,nsd_);
rhs = momentumSource_.quantity();
if (nodalGhostForce_) {
rhs -= nodalGhostForce_->quantity();
}
// (b) for ess. bcs
// form rhs : {sum_a (N_Ia * f_ia) - M_md * (ddupsilon/dt)_I}
DENS_MAT rhsPrescribed = -1.*nodalForce_.quantity();
atc_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
rhsPrescribed = (mdMassMatrix_.quantity())*rhsPrescribed;
rhsPrescribed += nodalAtomicForce_->quantity();
set<pair<int,int> >::const_iterator iter;
for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
rhs(iter->first,iter->second) = rhsPrescribed(iter->first,iter->second);
}
}
//--------------------------------------------------------
// compute_nodal_lambda_force
// computes the force induced on the FE
// by applying lambdaForce on the atoms
//--------------------------------------------------------
void StressFlux::compute_nodal_lambda_force(double dt)
{
DENS_MAT myNodalAtomicLambdaForce = nodalAtomicLambdaForce_->quantity();
set<pair<int,int> >::const_iterator iter;
for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
myNodalAtomicLambdaForce(iter->first,iter->second) = 0.;
}
timeFilter_->apply_post_step1(lambdaForceFiltered_->set_quantity(),
myNodalAtomicLambdaForce,dt);
}
//--------------------------------------------------------
// add_to_rhs:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the kinetostat
//--------------------------------------------------------
void StressFlux::add_to_rhs(FIELDS & rhs)
{
// compute the kinetostat force
compute_kinetostat(atc_->dt());
rhs[VELOCITY] += nodalAtomicLambdaForce_->quantity() + boundaryFlux_[VELOCITY].quantity();
}
//--------------------------------------------------------
// apply_lambda_to_atoms
// uses existing lambda to modify given
// atomic quantity
//--------------------------------------------------------
void StressFlux::apply_to_atoms(PerAtomQuantity<double> * atomVelocities,
const DENS_MAT & lambdaForce,
double dt)
{
_deltaVelocity_ = lambdaForce;
_deltaVelocity_ /= atomMasses_->quantity();
_deltaVelocity_ *= dt;
*atomVelocities += _deltaVelocity_;
}
//--------------------------------------------------------
// reset_filtered_ghost_force:
// resets the kinetostat generated ghost force to a
// prescribed value
//--------------------------------------------------------
void StressFlux::reset_filtered_ghost_force(DENS_MAT & target)
{
nodalGhostForceFiltered_ = target;
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void StressFlux::output(OUTPUT_LIST & outputData)
{
_nodalAtomicLambdaForceOut_ = nodalAtomicLambdaForce_->quantity();
DENS_MAT & lambda(lambda_->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = λ
outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_nodalAtomicLambdaForceOut_);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class StressFluxGhost
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
StressFluxGhost::StressFluxGhost(Kinetostat * kinetostat) :
StressFlux(kinetostat)
{
// flag for performing boundary flux calculation
fieldMask_(VELOCITY,FLUX) = false;
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void StressFluxGhost::construct_transfers()
{
StressFlux::construct_transfers();
if (!nodalGhostForce_) {
throw ATC_Error("StressFluxGhost::StressFluxGhost - ghost atoms must be specified");
}
}
//--------------------------------------------------------
// compute_boundary_flux:
// computes the boundary flux to be consistent with
// the controller
//--------------------------------------------------------
void StressFluxGhost::compute_boundary_flux(FIELDS & fields)
{
// This is only used in computation of atomic sources
boundaryFlux_[VELOCITY] = 0.;
}
//--------------------------------------------------------
// add_to_rhs:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the kinetostat
//--------------------------------------------------------
void StressFluxGhost::add_to_rhs(FIELDS & rhs)
{
// compute the kinetostat force
compute_kinetostat(atc_->dt());
// uses ghost force as the boundary flux to add to the RHS
rhs[VELOCITY] += nodalAtomicLambdaForce_->quantity() + nodalGhostForce_->quantity();
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the RHS of the kinetostat equations
// for the coupling parameter lambda
//--------------------------------------------------------
void StressFluxGhost::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// (a) for flux based :
// form rhs : \int N_I r dV - \sum_g N_Ig^* f_g
// sources are set in ATC transfer
rhs.reset(nNodes_,nsd_);
rhs = momentumSource_.quantity();
// (b) for ess. bcs
// form rhs : {sum_a (N_Ia * f_ia) - M_md * (ddupsilon/dt)_I}
DENS_MAT rhsPrescribed = -1.*nodalForce_.quantity();
atc_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
rhsPrescribed = (mdMassMatrix_.quantity())*rhsPrescribed;
rhsPrescribed += nodalAtomicForce_->quantity();
set<pair<int,int> >::const_iterator iter;
for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
rhs(iter->first,iter->second) = rhsPrescribed(iter->first,iter->second);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class StressFluxFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
StressFluxFiltered::StressFluxFiltered(Kinetostat * kinetostat) :
StressFlux(kinetostat),
nodalAtomicVelocity_(atc_->nodal_atomic_field(VELOCITY))
{
// do nothing
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the RHS of the kinetostat equations
// for the coupling parameter lambda
//--------------------------------------------------------
void StressFluxFiltered::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// set basic terms
// (a) for flux based :
// form rhs : \int N_I r dV - \sum_g N_Ig^* f_g
// sources are set in ATC transfer
rhs.reset(nNodes_,nsd_);
rhs = momentumSource_.quantity() - nodalGhostForceFiltered_.quantity();
// (b) for ess. bcs
// form rhs : {sum_a (N_Ia * f_ia) - M_md * (ddupsilon/dt)_I}
DENS_MAT rhsPrescribed = -1.*nodalForce_.quantity();
atc_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
rhsPrescribed = (mdMassMatrix_.quantity())*rhsPrescribed;
rhsPrescribed += nodalAtomicForce_->quantity();
set<pair<int,int> >::const_iterator iter;
for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
rhs(iter->first,iter->second) = rhsPrescribed(iter->first,iter->second);
}
// adjust for application of current lambda force
rhs += lambdaForceFiltered_->quantity();
// correct for time filtering
rhs *= 1./(timeFilter_->unfiltered_coefficient_pre_s1(dt));
}
//--------------------------------------------------------
// apply_lambda_to_atoms
// uses existing lambda to modify given
// atomic quantity
//--------------------------------------------------------
void StressFluxFiltered::apply_to_atoms(PerAtomQuantity<double> * atomVelocities,
const DENS_MAT & lambdaForce,
double dt)
{
StressFlux::apply_to_atoms(atomVelocities,lambdaForce,dt);
// add in corrections to filtered nodal atomice velocity
DENS_MAT velocityRoc(nNodes_,nsd_);
atc_->apply_inverse_md_mass_matrix(lambdaForceFiltered_->quantity(),
velocityRoc,
VELOCITY);
nodalAtomicVelocity_ += dt*velocityRoc;
}
//--------------------------------------------------------
// add_to_rhs:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the kinetostat
//--------------------------------------------------------
void StressFluxFiltered::add_to_rhs(FIELDS & rhs)
{
// compute kinetostat forces
compute_kinetostat(atc_->dt());
rhs[VELOCITY] += lambdaForceFiltered_->quantity() + boundaryFlux_[VELOCITY].quantity();
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void StressFluxFiltered::output(OUTPUT_LIST & outputData)
{
DENS_MAT & lambda(lambda_->set_quantity());
DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = λ
outputData[regulatorPrefix_+"NodalLambdaForce"] = &lambdaForceFiltered;
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class KinetostatGlcFs
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
KinetostatGlcFs::KinetostatGlcFs(Kinetostat * kinetostat,
const string & regulatorPrefix) :
KinetostatShapeFunction(kinetostat,regulatorPrefix),
velocity_(atc_->field(VELOCITY))
{
// constuct/obtain data corresponding to stage 3 of ATC_Method::initialize
nodalAtomicLambdaForce_ = kinetostat_->regulator_data(regulatorPrefix_+"NodalAtomicLambdaForce",nsd_);
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void KinetostatGlcFs::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// base class transfers
KinetostatShapeFunction::construct_transfers();
// get data from manager
nodalAtomicMomentum_ = interscaleManager.dense_matrix("NodalAtomicMomentum");
// atomic force induced by kinetostat
PerAtomQuantity<double> * atomLambdas = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaMomentum");
atomKinetostatForce_ = new AtomicKinetostatForceStress(atc_,atomLambdas);
interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
regulatorPrefix_+"AtomKinetostatForce");
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void KinetostatGlcFs::initialize()
{
KinetostatShapeFunction::initialize();
// set up workspaces
_deltaMomentum_.reset(nNodes_,nsd_);
_lambdaForceOutput_.reset(nNodes_,nsd_);
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (!timeFilterManager->end_equilibrate()) {
// we should reset lambda and lambdaForce to zero in this case
// implies an initial condition of 0 for the filtered nodal lambda power
// initial conditions will always be needed when using time filtering
// however, the fractional step scheme must assume the instantaneous
// nodal lambda power is 0 initially because all quantities are in delta form
*lambda_ = 0.; // ensures initial lambda force is zero
*nodalAtomicLambdaForce_ = 0.; // momentum change due to kinetostat
*lambdaForceFiltered_ = 0.; // filtered momentum change due to kinetostats
}
else {
// we can grab lambda power variables using time integrator and atc transfer in cases for equilibration
}
// sets up time filter for cases where variables temporally filtered
if (timeFilterManager->need_reset()) {
// the form of this integrator implies no time filters that require history data can be used
timeFilter_->initialize(nodalAtomicLambdaForce_->quantity());
}
compute_rhs_map();
}
//--------------------------------------------------------
// compute_rhs_map
// creates mapping from all nodes to those to which
// the kinetostat applies
//--------------------------------------------------------
void KinetostatGlcFs::compute_rhs_map()
{
rhsMap_.resize(overlapToNodeMap_->nRows(),1);
DENS_MAT rhsMapGlobal(nNodes_,1);
const set<int> & applicationNodes(applicationNodes_->quantity());
for (int i = 0; i < nNodes_; i++) {
if (applicationNodes.find(i) != applicationNodes.end()) {
rhsMapGlobal(i,0) = 1.;
}
else {
rhsMapGlobal(i,0) = 0.;
}
}
map_unique_to_overlap(rhsMapGlobal,rhsMap_);
}
//--------------------------------------------------------
// apply_pre_predictor:
// apply the kinetostat to the atoms in the
// pre-predictor integration phase
//--------------------------------------------------------
void KinetostatGlcFs::apply_pre_predictor(double dt)
{
DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
DENS_MAT & nodalAtomicLambdaForce(nodalAtomicLambdaForce_->set_quantity());
// update filtered forces
timeFilter_->apply_pre_step1(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
// apply lambda force to atoms and compute instantaneous lambda force
this->apply_to_atoms(atomVelocities_,nodalAtomicMomentum_,
atomKinetostatForce_->quantity(),
nodalAtomicLambdaForce,0.5*dt);
// update nodal variables for first half of timestep
this->add_to_momentum(nodalAtomicLambdaForce,_deltaMomentum_,0.5*dt);
atc_->apply_inverse_mass_matrix(_deltaMomentum_,VELOCITY);
velocity_ += _deltaMomentum_;
// start update of filtered lambda force
nodalAtomicLambdaForce = 0.;
timeFilter_->apply_post_step1(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the kinetostat to the atoms in the
// post-corrector integration phase
//--------------------------------------------------------
void KinetostatGlcFs::apply_post_corrector(double dt)
{
// compute the kinetostat equation and update lambda
this->compute_lambda(dt);
DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
DENS_MAT & nodalAtomicLambdaForce(nodalAtomicLambdaForce_->set_quantity());
// update filtered force
timeFilter_->apply_pre_step1(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
// apply lambda force to atoms and compute instantaneous lambda force
this->apply_to_atoms(atomVelocities_,nodalAtomicMomentum_,
atomKinetostatForce_->quantity(),
nodalAtomicLambdaForce,0.5*dt);
// update nodal variables for first half of timestep
this->add_to_momentum(nodalAtomicLambdaForce,_deltaMomentum_,0.5*dt);
nodalAtomicLambdaForce *= 2./dt;
atc_->apply_inverse_mass_matrix(_deltaMomentum_,VELOCITY);
velocity_ += _deltaMomentum_;
// start update of filtered lambda force
timeFilter_->apply_post_step2(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
}
//--------------------------------------------------------
// compute_kinetostat
// manages the solution and application of the
// kinetostat equations and variables
//--------------------------------------------------------
void KinetostatGlcFs::compute_lambda(double dt)
{
// set up rhs for lambda equation
this->set_kinetostat_rhs(rhs_,0.5*dt);
// solve linear system for lambda
DENS_MAT & lambda(lambda_->set_quantity());
solve_for_lambda(rhs_,lambda);
}
//--------------------------------------------------------
// apply_lambda_to_atoms
// uses existing lambda to modify given
// atomic quantity
//--------------------------------------------------------
void KinetostatGlcFs::apply_to_atoms(PerAtomQuantity<double> * atomVelocity,
const DENS_MAN * nodalAtomicMomentum,
const DENS_MAT & lambdaForce,
DENS_MAT & nodalAtomicLambdaForce,
double dt)
{
// compute initial contributions to lambda force
nodalAtomicLambdaForce = nodalAtomicMomentum->quantity();
nodalAtomicLambdaForce *= -1.;
// apply lambda force to atoms
_velocityDelta_ = lambdaForce;
_velocityDelta_ /= atomMasses_->quantity();
_velocityDelta_ *= dt;
(*atomVelocity) += _velocityDelta_;
// finalize lambda force
nodalAtomicLambdaForce += nodalAtomicMomentum->quantity();
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void KinetostatGlcFs::output(OUTPUT_LIST & outputData)
{
_lambdaForceOutput_ = nodalAtomicLambdaForce_->quantity();
// approximate value for lambda force
double dt = LammpsInterface::instance()->dt();
_lambdaForceOutput_ *= (2./dt);
DENS_MAT & lambda(lambda_->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = λ
outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_lambdaForceOutput_);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class KinetostatFlux
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
KinetostatFlux::KinetostatFlux(Kinetostat * kinetostat,
const string & regulatorPrefix) :
KinetostatGlcFs(kinetostat,regulatorPrefix),
momentumSource_(atc_->atomic_source(VELOCITY)),
nodalGhostForce_(NULL),
nodalGhostForceFiltered_(NULL)
{
// flag for performing boundary flux calculation
fieldMask_(VELOCITY,FLUX) = true;
// constuct/obtain data corresponding to stage 3 of ATC_Method::initialize
nodalGhostForceFiltered_ = kinetostat_->regulator_data(regulatorPrefix_+"NodalGhostForceFiltered",nsd_);
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void KinetostatFlux::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// set up linear solver
// set up data for linear solver
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixMomentum");
if (elementMask_) {
lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
regulatorPrefix_+"LambdaAtomMap");
}
if (atomicRegulator_->use_localized_lambda()) {
linearSolverType_ = AtomicRegulator::RSL_SOLVE;
}
else {
linearSolverType_ = AtomicRegulator::CG_SOLVE;
}
// base class transfers
KinetostatGlcFs::construct_transfers();
// sets up space for ghost force related variables
if (atc_->groupbit_ghost()) {
MatrixDependencyManager<DenseMatrix, int> * nodeToOverlapMap =
interscaleManager.dense_matrix_int(regulatorPrefix_+"NodeToOverlapMap");
GhostCouplingMatrix * shapeFunctionGhost = new GhostCouplingMatrix(atc_,interscaleManager.per_atom_sparse_matrix("InterpolantGhost"),
regulatedNodes_,
nodeToOverlapMap);
interscaleManager.add_sparse_matrix(shapeFunctionGhost,
regulatorPrefix_+"GhostCouplingMatrix");
FundamentalAtomQuantity * atomGhostForce = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE,
GHOST);
nodalGhostForce_ = new AtfShapeFunctionRestriction(atc_,atomGhostForce,
shapeFunctionGhost);
interscaleManager.add_dense_matrix(nodalGhostForce_,
regulatorPrefix_+"NodalGhostForce");
}
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void KinetostatFlux::initialize()
{
KinetostatGlcFs::initialize();
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (!timeFilterManager->end_equilibrate()) {
// we should reset lambda and lambdaForce to zero in this case
// implies an initial condition of 0 for the filtered nodal lambda power
// initial conditions will always be needed when using time filtering
// however, the fractional step scheme must assume the instantaneous
// nodal lambda power is 0 initially because all quantities are in delta form
*nodalGhostForceFiltered_ = 0.; // filtered force from ghost atoms
}
else {
// we can grab lambda power variables using time integrator and atc transfer in cases for equilibration
}
}
//--------------------------------------------------------
// construct_regulated_nodes:
// constructs the set of nodes being regulated
//--------------------------------------------------------
void KinetostatFlux::construct_regulated_nodes()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// matrix requires all entries even if localized for correct lumping
regulatedNodes_ = new RegulatedNodes(atc_);
interscaleManager.add_set_int(regulatedNodes_,
regulatorPrefix_+"KinetostatRegulatedNodes");
// if localized monitor nodes with applied fluxes
if (atomicRegulator_->use_localized_lambda()) {
if ((kinetostat_->coupling_mode() == Kinetostat::FLUX) && (atomicRegulator_->boundary_integration_type() != NO_QUADRATURE)) {
// include boundary nodes
applicationNodes_ = new FluxBoundaryNodes(atc_);
boundaryNodes_ = new BoundaryNodes(atc_);
interscaleManager.add_set_int(boundaryNodes_,
regulatorPrefix_+"KinetostatBoundaryNodes");
}
else {
// fluxed nodes only
applicationNodes_ = new FluxNodes(atc_);
}
interscaleManager.add_set_int(applicationNodes_,
regulatorPrefix_+"KinetostatApplicationNodes");
}
else {
applicationNodes_ = regulatedNodes_;
}
// special set of boundary elements for boundary flux quadrature
if ((atomicRegulator_->boundary_integration_type() == FE_INTERPOLATION)
&& (atomicRegulator_->use_localized_lambda())) {
elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
interscaleManager.add_dense_matrix_bool(elementMask_,
regulatorPrefix_+"BoundaryElementMask");
}
}
//--------------------------------------------------------
// apply_pre_predictor:
// apply the kinetostat to the atoms in the
// pre-predictor integration phase
//--------------------------------------------------------
void KinetostatFlux::apply_pre_predictor(double dt)
{
// update filtered forces
if (nodalGhostForce_) {
timeFilter_->apply_pre_step1(nodalGhostForceFiltered_->set_quantity(),
nodalGhostForce_->quantity(),dt);
}
KinetostatGlcFs::apply_pre_predictor(dt);
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the kinetostat to the atoms in the
// post-corrector integration phase
//--------------------------------------------------------
void KinetostatFlux::apply_post_corrector(double dt)
{
// update filtered ghost force
if (nodalGhostForce_) {
timeFilter_->apply_post_step1(nodalGhostForceFiltered_->set_quantity(),
nodalGhostForce_->quantity(),dt);
}
// compute the kinetostat equation and update lambda
KinetostatGlcFs::apply_post_corrector(dt);
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the RHS of the kinetostat equations
// for the coupling parameter lambda
//--------------------------------------------------------
void KinetostatFlux::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// (a) for flux based :
// form rhs : \int N_I r dV - \sum_g N_Ig^* f_g
// sources are set in ATC transfer
rhs.reset(nNodes_,nsd_);
const DENS_MAT & momentumSource(momentumSource_.quantity());
const set<int> & applicationNodes(applicationNodes_->quantity());
set<int>::const_iterator iNode;
for (iNode = applicationNodes.begin(); iNode != applicationNodes.end(); iNode++) {
for (int j = 0; j < nsd_; j++) {
rhs(*iNode,j) = momentumSource(*iNode,j);
}
}
// add ghost forces, if needed
if (nodalGhostForce_) {
const DENS_MAT & nodalGhostForce(nodalGhostForce_->quantity());
for (iNode = applicationNodes.begin(); iNode != applicationNodes.end(); iNode++) {
for (int j = 0; j < nsd_; j++) {
rhs(*iNode,j) -= nodalGhostForce(*iNode,j);
}
}
}
}
//--------------------------------------------------------
// add_to_momentum:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the kinetostat
//--------------------------------------------------------
void KinetostatFlux::add_to_momentum(const DENS_MAT & nodalLambdaForce,
DENS_MAT & deltaMomentum,
double dt)
{
deltaMomentum.resize(nNodes_,nsd_);
const DENS_MAT & boundaryFlux(boundaryFlux_[VELOCITY].quantity());
for (int i = 0; i < nNodes_; i++) {
for (int j = 0; j < nsd_; j++) {
deltaMomentum(i,j) = nodalLambdaForce(i,j) + dt*boundaryFlux(i,j);
}
}
}
//--------------------------------------------------------
// reset_filtered_ghost_force:
// resets the kinetostat generated ghost force to a
// prescribed value
//--------------------------------------------------------
void KinetostatFlux::reset_filtered_ghost_force(DENS_MAT & target)
{
(*nodalGhostForceFiltered_) = target;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class KinetostatFluxGhost
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
KinetostatFluxGhost::KinetostatFluxGhost(Kinetostat * kinetostat,
const string & regulatorPrefix) :
KinetostatFlux(kinetostat,regulatorPrefix)
{
// flag for performing boundary flux calculation
fieldMask_(VELOCITY,FLUX) = false;
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void KinetostatFluxGhost::construct_transfers()
{
KinetostatFlux::construct_transfers();
if (!nodalGhostForce_) {
throw ATC_Error("StressFluxGhost::StressFluxGhost - ghost atoms must be specified");
}
}
//--------------------------------------------------------
// compute_boundary_flux:
// computes the boundary flux to be consistent with
// the controller
//--------------------------------------------------------
void KinetostatFluxGhost::compute_boundary_flux(FIELDS & fields)
{
// This is only used in computation of atomic sources
boundaryFlux_[VELOCITY] = 0.;
}
//--------------------------------------------------------
// add_to_momentum:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the kinetostat
//--------------------------------------------------------
void KinetostatFluxGhost::add_to_momentum(const DENS_MAT & nodalLambdaForce,
DENS_MAT & deltaMomentum,
double dt)
{
deltaMomentum.resize(nNodes_,nsd_);
const DENS_MAT & boundaryFlux(nodalGhostForce_->quantity());
for (int i = 0; i < nNodes_; i++) {
for (int j = 0; j < nsd_; j++) {
deltaMomentum(i,j) = nodalLambdaForce(i,j) + dt*boundaryFlux(i,j);
}
}
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the RHS of the kinetostat equations
// for the coupling parameter lambda
//--------------------------------------------------------
void KinetostatFluxGhost::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// (a) for flux based :
// form rhs : \int N_I r dV - \sum_g N_Ig^* f_g
// sources are set in ATC transfer
rhs.reset(nNodes_,nsd_);
const DENS_MAT & momentumSource(momentumSource_.quantity());
const set<int> & applicationNodes(applicationNodes_->quantity());
set<int>::const_iterator iNode;
for (iNode = applicationNodes.begin(); iNode != applicationNodes.end(); iNode++) {
for (int j = 0; j < nsd_; j++) {
rhs(*iNode,j) = momentumSource(*iNode,j);
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class KinetostatFixed
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and kinetostat data
//--------------------------------------------------------
KinetostatFixed::KinetostatFixed(Kinetostat * kinetostat,
const string & regulatorPrefix) :
KinetostatGlcFs(kinetostat,regulatorPrefix),
mdMassMatrix_(atc_->set_mass_mat_md(VELOCITY)),
isFirstTimestep_(true)
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void KinetostatFixed::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// determine if map is needed and set up if so
if (this->use_local_shape_functions()) {
lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
regulatorPrefix_+"LambdaAtomMap");
shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
lambdaAtomMap_,
nodeToOverlapMap_);
}
else {
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
}
interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixMomentum");
linearSolverType_ = AtomicRegulator::CG_SOLVE;
// base class transfers
KinetostatGlcFs::construct_transfers();
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void KinetostatFixed::initialize()
{
KinetostatGlcFs::initialize();
// reset data to zero
deltaFeMomentum_.reset(nNodes_,nsd_);
deltaNodalAtomicMomentum_.reset(nNodes_,nsd_);
}
//--------------------------------------------------------
// halve_force:
// flag to halve the lambda force for improved
// accuracy
//--------------------------------------------------------
bool KinetostatFixed::halve_force()
{
if (isFirstTimestep_ || ((atc_->atom_to_element_map_type() == EULERIAN)
&& (atc_->atom_to_element_map_frequency() > 1)
&& (atc_->step() % atc_->atom_to_element_map_frequency() == 1))) {
return true;
}
return false;
}
//--------------------------------------------------------
// construct_regulated_nodes:
// constructs the set of nodes being regulated
//--------------------------------------------------------
void KinetostatFixed::construct_regulated_nodes()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
if (!atomicRegulator_->use_localized_lambda()) {
regulatedNodes_ = new RegulatedNodes(atc_);
}
else if (kinetostat_->coupling_mode() == Kinetostat::FLUX) {
regulatedNodes_ = new FixedNodes(atc_);
}
else if (kinetostat_->coupling_mode() == Kinetostat::FIXED) {
// include boundary nodes
regulatedNodes_ = new FixedBoundaryNodes(atc_);
}
else {
throw ATC_Error("ThermostatFixed::construct_regulated_nodes - couldn't determine set of regulated nodes");
}
interscaleManager.add_set_int(regulatedNodes_,
regulatorPrefix_+"RegulatedNodes");
applicationNodes_ = regulatedNodes_;
// special set of boundary elements for defining regulated atoms
if (atomicRegulator_->use_localized_lambda()) {
elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
interscaleManager.add_dense_matrix_bool(elementMask_,
regulatorPrefix_+"BoundaryElementMask");
}
}
//--------------------------------------------------------
// initialize_delta_nodal_atomic_momentum:
// initializes storage for the variable tracking
// the change in the nodal atomic momentum
// that has occured over the past timestep
//--------------------------------------------------------
void KinetostatFixed::initialize_delta_nodal_atomic_momentum(double dt)
{
// initialize delta energy
const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity());
initialNodalAtomicMomentum_ = myNodalAtomicMomentum;
initialNodalAtomicMomentum_ *= -1.; // initially stored as negative for efficiency
timeFilter_->apply_pre_step1(nodalAtomicMomentumFiltered_.set_quantity(),
myNodalAtomicMomentum,dt);
}
//--------------------------------------------------------
// compute_delta_nodal_atomic_momentum:
// computes the change in the nodal atomic momentum
// that has occured over the past timestep
//--------------------------------------------------------
void KinetostatFixed::compute_delta_nodal_atomic_momentum(double dt)
{
// set delta energy based on predicted atomic velocities
const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity());
timeFilter_->apply_post_step1(nodalAtomicMomentumFiltered_.set_quantity(),
myNodalAtomicMomentum,dt);
deltaNodalAtomicMomentum_ = initialNodalAtomicMomentum_;
deltaNodalAtomicMomentum_ += myNodalAtomicMomentum;
}
//--------------------------------------------------------
// apply_pre_predictor:
// apply the kinetostat to the atoms in the
// pre-predictor integration phase
//--------------------------------------------------------
void KinetostatFixed::apply_pre_predictor(double dt)
{
// initialize values to be track change in finite element energy over the timestep
initialize_delta_nodal_atomic_momentum(dt);
initialFeMomentum_ = -1.*((mdMassMatrix_.quantity())*(velocity_.quantity())); // initially stored as negative for efficiency
KinetostatGlcFs::apply_pre_predictor(dt);
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the kinetostat to the atoms in the
// post-corrector integration phase
//--------------------------------------------------------
void KinetostatFixed::apply_post_corrector(double dt)
{
KinetostatGlcFs::apply_post_corrector(dt);
// update filtered momentum with lambda force
DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->set_quantity());
timeFilter_->apply_post_step2(nodalAtomicMomentumFiltered_.set_quantity(),
myNodalAtomicLambdaForce,dt);
if (halve_force()) {
// Halve lambda force due to fixed temperature constraints
// 1) makes up for poor initial condition
// 2) accounts for possibly large value of lambda when atomic shape function values change
// from eulerian mapping after more than 1 timestep
// avoids unstable oscillations arising from
// thermostat having to correct for error introduced in lambda changing the
// shape function matrices
*lambda_ *= 0.5;
}
isFirstTimestep_ = false;
}
//--------------------------------------------------------
// compute_kinetostat
// manages the solution and application of the
// kinetostat equations and variables
//--------------------------------------------------------
void KinetostatFixed::compute_lambda(double dt)
{
// compute predicted changes in nodal atomic momentum
compute_delta_nodal_atomic_momentum(dt);
// change in finite element momentum
deltaFeMomentum_ = initialFeMomentum_;
deltaFeMomentum_ += (mdMassMatrix_.quantity())*(velocity_.quantity());
// set up rhs for lambda equation
KinetostatGlcFs::compute_lambda(dt);
}
//--------------------------------------------------------
// set_kinetostat_rhs
// sets up the RHS of the kinetostat equations
// for the coupling parameter lambda
//--------------------------------------------------------
void KinetostatFixed::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
{
// for essential bcs (fixed nodes) :
// form rhs : (delUpsV - delUps)/dt
const set<int> & regulatedNodes(regulatedNodes_->quantity());
double factor = (1./dt);
for (int i = 0; i < nNodes_; i++) {
if (regulatedNodes.find(i) != regulatedNodes.end()) {
for (int j = 0; j < nsd_; j++) {
rhs(i,j) = factor*(deltaNodalAtomicMomentum_(i,j) - deltaFeMomentum_(i,j));
}
}
else {
for (int j = 0; j < nsd_; j++) {
rhs(i,j) = 0.;
}
}
}
}
//--------------------------------------------------------
// add_to_momentum:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the kinetostat
//--------------------------------------------------------
void KinetostatFixed::add_to_momentum(const DENS_MAT & nodalLambdaForce,
DENS_MAT & deltaMomentum,
double dt)
{
deltaMomentum.resize(nNodes_,nsd_);
const set<int> & regulatedNodes(regulatedNodes_->quantity());
for (int i = 0; i < nNodes_; i++) {
if (regulatedNodes.find(i) != regulatedNodes.end()) {
for (int j = 0; j < nsd_; j++) {
deltaMomentum(i,j) = 0.;
}
}
else {
for (int j = 0; j < nsd_; j++) {
deltaMomentum(i,j) = nodalLambdaForce(i,j);
}
}
}
}
};
Event Timeline
Log In to Comment