Homec4science
Diffusion lammps (lammm-master-deprecated) Browse

lammps/lib/atc/Thermostat.cpp
494ee3b26c32lammm-master-deprecated

Last Change
Raw File

Thermostat.cpp
View Options

#include "Thermostat.h"
#include "ATC_Coupling.h"
#include "ATC_Error.h"
#include "PrescribedDataManager.h"
#include "ThermalTimeIntegrator.h"
#include "TransferOperator.h"
using namespace std;
namespace ATC {
//--------------------------------------------------------
//--------------------------------------------------------
// Class Thermostat
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
Thermostat::Thermostat(ATC_Coupling * atc,
const string & regulatorPrefix) :
AtomicRegulator(atc,regulatorPrefix),
lambdaMaxIterations_(myLambdaMaxIterations)
{
// nothing to do
}
//--------------------------------------------------------
// modify:
// parses and adjusts thermostat state based on
// user input, in the style of LAMMPS user input
//--------------------------------------------------------
bool Thermostat::modify(int narg, char **arg)
{
bool foundMatch = false;
int argIndex = 0;
if (strcmp(arg[argIndex],"thermal")==0) {
argIndex++;
// thermostat type
/*! \page man_control_thermal fix_modify AtC control thermal
\section syntax
fix_modify AtC control thermal <control_type> <optional_args>
- control_type (string) = none | rescale | hoover | flux\n
fix_modify AtC control thermal rescale <frequency> \n
- frequency (int) = time step frequency for applying velocity rescaling \n
fix_modify AtC control thermal hoover \n
fix_modify AtC control thermal flux <boundary_integration_type(optional)> <face_set_id(optional)>\n
- boundary_integration_type (string) = faceset | interpolate\n
- face_set_id (string), optional = 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, does not work with interpolate\n
\section examples
<TT> fix_modify AtC control thermal none </TT> \n
<TT> fix_modify AtC control thermal rescale 10 </TT> \n
<TT> fix_modify AtC control thermal hoover </TT> \n
<TT> fix_modify AtC control thermal flux </TT> \n
<TT> fix_modify AtC control thermal flux faceset bndy_faces </TT> \n
\section description
Sets the energy exchange mechansim from the finite elements to the atoms, managed through a control algorithm. Rescale computes a scale factor for each atom to match the finite element temperature. Hoover is a Gaussian least-constraint isokinetic thermostat enforces that the nodal restricted atomic temperature matches the finite element temperature. Flux is a similar mode, but rather adds energy to the atoms based on conservation of energy. Hoover and flux allows the prescription of sources or fixed temperatures on the atoms.
\section restrictions
only for be used with specific transfers :
thermal (rescale, hoover, flux), two_temperature (flux) \n
rescale not valid with time filtering activated
\section related
\section default
none\n
rescale frequency is 1\n
flux boundary_integration_type is interpolate
*/
if (strcmp(arg[argIndex],"none")==0) { // restore defaults
regulatorTarget_ = NONE;
couplingMode_ = UNCOUPLED;
howOften_ = 1;
boundaryIntegrationType_ = NO_QUADRATURE;
foundMatch = true;
}
else if (strcmp(arg[argIndex],"rescale")==0) {
argIndex++;
howOften_ = atoi(arg[argIndex]);
if (howOften_ < 1) {
throw ATC_Error("Bad rescaling thermostat frequency");
}
else {
regulatorTarget_ = FIELD;
couplingMode_ = UNCOUPLED;
boundaryIntegrationType_ = NO_QUADRATURE;
foundMatch = true;
}
}
else if (strcmp(arg[argIndex],"hoover")==0) {
regulatorTarget_ = DYNAMICS;
couplingMode_ = FIXED;
howOften_ = 1;
boundaryIntegrationType_ = NO_QUADRATURE;
foundMatch = true;
}
else if (strcmp(arg[argIndex],"flux")==0) {
regulatorTarget_ = DYNAMICS;
couplingMode_ = FLUX;
howOften_ = 1;
argIndex++;
boundaryIntegrationType_ = atc_->parse_boundary_integration(narg-argIndex,&arg[argIndex],boundaryFaceSet_);
foundMatch = true;
}
// set parameters for numerical matrix solutions unique to this thermostat
/*! \page man_control_thermal_correction_max_iterations fix_modify AtC control thermal correction_max_iterations
\section syntax
fix_modify AtC control thermal correction_max_iterations <max_iterations>
- max_iterations (int) = maximum number of iterations that will be used by iterative matrix solvers\n
\section examples
<TT> fix_modify AtC control thermal correction_max_iterations 10 </TT> \n
\section description
Sets the maximum number of iterations to compute the 2nd order in time correction term for lambda with the fractional step method. The method uses the same tolerance as the controller's matrix solver.
\section restrictions
only for use with thermal physics using the fractional step method.
\section related
\section default
correction_max_iterations is 20
*/
else if (strcmp(arg[argIndex],"correction_max_iterations")==0) {
argIndex++;
lambdaMaxIterations_ = atoi(arg[argIndex]);
if (lambdaMaxIterations_ < 1) {
throw ATC_Error("Bad correction maximum iteration count");
}
foundMatch = true;
}
}
if (!foundMatch)
foundMatch = AtomicRegulator::modify(narg,arg);
if (foundMatch)
needReset_ = true;
return foundMatch;
}
//--------------------------------------------------------
// reset_lambda_contribution:
// resets the thermostat generated power to a
// prescribed value
//--------------------------------------------------------
void Thermostat::reset_lambda_contribution(const DENS_MAT & target)
{
DENS_MAN * lambdaPowerFiltered = regulator_data("LambdaPowerFiltered",1);
*lambdaPowerFiltered = target;
}
//--------------------------------------------------------
// construct_methods:
// instantiations desired regulator method(s)
// dependence, but in general there is also a
// time integrator dependence. In general the
// precedence order is:
// time filter -> time integrator -> thermostat
// In the future this may need to be added if
// different types of time integrators can be
// specified.
//--------------------------------------------------------
void Thermostat::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();
// update time filter
TimeIntegrator::TimeIntegrationType myIntegrationType = (atc_->time_integrator(TEMPERATURE))->time_integration_type();
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (timeFilterManager->need_reset() ) {
if (myIntegrationType == TimeIntegrator::GEAR) {
timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT);
}
else if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT_IMPLICIT);
}
}
if (timeFilterManager->filter_dynamics()) {
switch (regulatorTarget_) {
case NONE: {
regulatorMethod_ = new RegulatorMethod(this);
break;
}
case FIELD: { // error check, rescale and filtering not supported together
throw ATC_Error("Cannot use rescaling thermostat with time filtering");
break;
}
case DYNAMICS: {
switch (couplingMode_) {
case FIXED: {
if (use_lumped_lambda_solve()) {
throw ATC_Error("Thermostat:construct_methods - lumped lambda solve cannot be used with Hoover thermostats");
}
if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
if (md_flux_nodes(TEMPERATURE)) {
if (!md_fixed_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
// there are fluxes but no fixed or coupled nodes
regulatorMethod_ = new ThermostatIntegratorFluxFiltered(this,lambdaMaxIterations_);
}
else {
// there are both fixed and flux nodes
regulatorMethod_ = new ThermostatFluxFixedFiltered(this,lambdaMaxIterations_);
}
}
else {
// there are only fixed nodes
regulatorMethod_ = new ThermostatIntegratorFixedFiltered(this,lambdaMaxIterations_);
}
}
else {
regulatorMethod_ = new ThermostatHooverVerletFiltered(this);
}
break;
}
case FLUX: {
if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
if (use_lumped_lambda_solve()) {
throw ATC_Error("Thermostat:construct_methods - lumped lambda solve has been depricated for fractional step thermostats");
}
if (md_fixed_nodes(TEMPERATURE)) {
if (!md_flux_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
// there are fixed nodes but no fluxes
regulatorMethod_ = new ThermostatIntegratorFixedFiltered(this,lambdaMaxIterations_);
}
else {
// there are both fixed and flux nodes
regulatorMethod_ = new ThermostatFluxFixedFiltered(this,lambdaMaxIterations_);
}
}
else {
// there are only flux nodes
regulatorMethod_ = new ThermostatIntegratorFluxFiltered(this,lambdaMaxIterations_);
}
}
else {
if (use_localized_lambda()) {
if (!((atc_->prescribed_data_manager())->no_fluxes(TEMPERATURE)) &&
atc_->boundary_integration_type() != NO_QUADRATURE) {
throw ATC_Error("Cannot use flux coupling with localized lambda");
}
}
regulatorMethod_ = new ThermostatPowerVerletFiltered(this);
}
break;
}
default:
throw ATC_Error("Unknown coupling mode in Thermostat::initialize");
}
break;
}
default:
throw ATC_Error("Unknown thermostat type in Thermostat::initialize");
}
}
else {
switch (regulatorTarget_) {
case NONE: {
regulatorMethod_ = new RegulatorMethod(this);
break;
}
case FIELD: {
if (atc_->temperature_def()==KINETIC)
regulatorMethod_ = new ThermostatRescale(this);
else if (atc_->temperature_def()==TOTAL)
regulatorMethod_ = new ThermostatRescaleMixedKePe(this);
else
throw ATC_Error("Unknown temperature definition");
break;
}
case DYNAMICS: {
switch (couplingMode_) {
case FIXED: {
if (use_lumped_lambda_solve()) {
throw ATC_Error("Thermostat:construct_methods - lumped lambda solve cannot be used with Hoover thermostats");
}
if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
if (md_flux_nodes(TEMPERATURE)) {
if (!md_fixed_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
// there are fluxes but no fixed or coupled nodes
regulatorMethod_ = new ThermostatIntegratorFlux(this,lambdaMaxIterations_);
}
else {
// there are both fixed and flux nodes
regulatorMethod_ = new ThermostatFluxFixed(this,lambdaMaxIterations_);
}
}
else {
// there are only fixed nodes
regulatorMethod_ = new ThermostatIntegratorFixed(this,lambdaMaxIterations_);
}
}
else {
regulatorMethod_ = new ThermostatHooverVerlet(this);
}
break;
}
case FLUX: {
if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
if (use_lumped_lambda_solve()) {
throw ATC_Error("Thermostat:construct_methods - lumped lambda solve has been depricated for fractional step thermostats");
}
if (md_fixed_nodes(TEMPERATURE)) {
if (!md_flux_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
// there are fixed nodes but no fluxes
regulatorMethod_ = new ThermostatIntegratorFixed(this,lambdaMaxIterations_);
}
else {
// there are both fixed and flux nodes
regulatorMethod_ = new ThermostatFluxFixed(this,lambdaMaxIterations_);
}
}
else {
// there are only flux nodes
regulatorMethod_ = new ThermostatIntegratorFlux(this,lambdaMaxIterations_);
}
}
else {
if (use_localized_lambda()) {
if (!((atc_->prescribed_data_manager())->no_fluxes(TEMPERATURE)) &&
atc_->boundary_integration_type() != NO_QUADRATURE) {
throw ATC_Error("Cannot use flux coupling with localized lambda");
}
}
regulatorMethod_ = new ThermostatPowerVerlet(this);
}
break;
}
default:
throw ATC_Error("Unknown coupling mode in Thermostat::initialize");
}
break;
}
default:
throw ATC_Error("Unknown thermostat target in Thermostat::initialize");
}
}
AtomicRegulator::reset_method();
}
else {
set_all_data_to_used();
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatShapeFunction
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatShapeFunction::ThermostatShapeFunction(AtomicRegulator * thermostat,
const string & regulatorPrefix) :
RegulatorShapeFunction(thermostat,regulatorPrefix),
mdMassMatrix_(atc_->set_mass_mat_md(TEMPERATURE)),
atomVelocities_(NULL)
{
fieldMask_(TEMPERATURE,FLUX) = true;
lambda_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1); // data associated with stage 3 in ATC_Method::initialize
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatShapeFunction::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
RegulatorShapeFunction::construct_transfers();
// get atom velocity data from manager
atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
// construct lambda evaluated at atom locations
atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
lambda_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
}
//---------------------------------------------------------
// set_weights:
// set the diagonal weighting matrix to be the atomic
// temperatures
//---------------------------------------------------------
void ThermostatShapeFunction::set_weights()
{
if (this->use_local_shape_functions()) {
VelocitySquaredMapped * myWeights = new VelocitySquaredMapped(atc_,lambdaAtomMap_);
weights_ = myWeights;
(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
regulatorPrefix_+"AtomVelocitySquaredMapped");
}
else {
VelocitySquared * myWeights = new VelocitySquared(atc_);
weights_ = myWeights;
(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
regulatorPrefix_+"AtomVelocitySquared");
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatRescale
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatRescale::ThermostatRescale(AtomicRegulator * thermostat) :
ThermostatShapeFunction(thermostat),
nodalTemperature_(atc_->field(TEMPERATURE)),
atomVelocityRescalings_(NULL)
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatRescale::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// set up data for linear solver
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
(atc_->interscale_manager()).add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixEnergy");
linearSolverType_ = AtomicRegulator::CG_SOLVE;
// base class transfers
ThermostatShapeFunction::construct_transfers();
// velocity rescaling factor
atomVelocityRescalings_ = new AtomicVelocityRescaleFactor(atc_,atomLambdas_);
interscaleManager.add_per_atom_quantity(atomVelocityRescalings_,
regulatorPrefix_+"AtomVelocityRescaling");
}
//---------------------------------------------------------
// set_weights:
// set the diagonal weighting matrix to be the atomic
// temperatures
//---------------------------------------------------------
void ThermostatRescale::set_weights()
{
weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicEnergyForTemperature");
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the thermostat in the post corrector phase
//--------------------------------------------------------
void ThermostatRescale::apply_post_corrector(double dt)
{
compute_thermostat(dt);
// application of rescaling lambda due
apply_to_atoms(atomVelocities_);
}
//--------------------------------------------------------
// compute_thermostat
// manages the solution of the
// thermostat equations and variables
//--------------------------------------------------------
void ThermostatRescale::compute_thermostat(double dt)
{
// compute right-hand side
this->set_rhs(_rhs_);
// solve equations
solve_for_lambda(_rhs_,lambda_->set_quantity());
}
//--------------------------------------------------------
// set_rhs:
// constructs the RHS vector with the target
// temperature
//--------------------------------------------------------
void ThermostatRescale::set_rhs(DENS_MAT & rhs)
{
rhs = mdMassMatrix_.quantity()*nodalTemperature_.quantity();
}
//--------------------------------------------------------
// apply_lambda_to_atoms:
// applies the velocity rescale with an existing lambda
// note oldAtomicQuantity and dt are not used
//--------------------------------------------------------
void ThermostatRescale::apply_to_atoms(PerAtomQuantity<double> * atomVelocities)
{
*atomVelocities *= atomVelocityRescalings_->quantity();
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void ThermostatRescale::output(OUTPUT_LIST & outputData)
{
DENS_MAT & lambda(lambda_->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData["LambdaEnergy"] = &lambda;
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatRescaleMixedKePe
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatRescaleMixedKePe::ThermostatRescaleMixedKePe(AtomicRegulator * thermostat) :
ThermostatRescale(thermostat),
nodalAtomicFluctuatingPotentialEnergy_(NULL)
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatRescaleMixedKePe::construct_transfers()
{
ThermostatRescale::construct_transfers();
InterscaleManager & interscaleManager(atc_->interscale_manager());
// get fluctuating PE at nodes
nodalAtomicFluctuatingPotentialEnergy_ =
interscaleManager.dense_matrix("NodalAtomicFluctuatingPotentialEnergy");
}
//---------------------------------------------------------
// set_weights:
// set the diagonal weighting matrix to be the atomic
// temperatures
//---------------------------------------------------------
void ThermostatRescaleMixedKePe::set_weights()
{
weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicTwiceKineticEnergy");
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void ThermostatRescaleMixedKePe::initialize()
{
ThermostatRescale::initialize();
InterscaleManager & interscaleManager(atc_->interscale_manager());
// multipliers for KE and PE
AtomicEnergyForTemperature * atomEnergyForTemperature =
static_cast<AtomicEnergyForTemperature * >(interscaleManager.per_atom_quantity("AtomicEnergyForTemperature"));
keMultiplier_ = atomEnergyForTemperature->kinetic_energy_multiplier();
peMultiplier_ = 2. - keMultiplier_;
keMultiplier_ /= 2.; // account for use of 2 X KE in matrix equation
}
//--------------------------------------------------------
// set_rhs:
// accounts for potential energy contribution to
// definition of atomic temperature
//--------------------------------------------------------
void ThermostatRescaleMixedKePe::set_rhs(DENS_MAT & rhs)
{
ThermostatRescale::set_rhs(rhs);
rhs -= peMultiplier_*(nodalAtomicFluctuatingPotentialEnergy_->quantity());
rhs /= keMultiplier_;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatFsSolver
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatFsSolver::ThermostatFsSolver(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
RegulatorShapeFunction(thermostat,regulatorPrefix),
lambdaMaxIterations_(lambdaMaxIterations),
rhsLambdaSquared_(NULL),
dtFactor_(1.)
{
fieldMask_(TEMPERATURE,FLUX) = true;
lambda_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1); // data associated with stage 3 in ATC_Method::initialize
}
//--------------------------------------------------------
// initialize
// creates mapping from all nodes to those to which
// the thermostat applies
//--------------------------------------------------------
void ThermostatFsSolver::initialize()
{
RegulatorShapeFunction::initialize();
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_);
}
//---------------------------------------------------------
// set_weights:
// set the diagonal weighting matrix to be the atomic
// temperatures
//---------------------------------------------------------
void ThermostatFsSolver::set_weights()
{
if (this->use_local_shape_functions()) {
VelocitySquaredMapped * myWeights = new VelocitySquaredMapped(atc_,lambdaAtomMap_);
weights_ = myWeights;
(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
regulatorPrefix_+"AtomVelocitySquaredMapped");
}
else {
VelocitySquared * myWeights = new VelocitySquared(atc_);
weights_ = myWeights;
(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
regulatorPrefix_+"AtomVelocitySquared");
}
}
//--------------------------------------------------------
// compute_lambda:
// solves linear system for lambda, if the
// bool is true it iterators to a non-linear solution
//--------------------------------------------------------
void ThermostatFsSolver::compute_lambda(const DENS_MAT & rhs,
bool iterateSolution)
{
// solve linear system for lambda guess
DENS_MAT & lambda(lambda_->set_quantity());
solve_for_lambda(rhs,lambda);
// iterate to solution
if (iterateSolution) {
iterate_lambda(rhs);
}
}
//--------------------------------------------------------
// iterate_lambda:
// iteratively solves the equations for lambda
// for the higher order dt corrections, assuming
// an initial guess for lambda
//--------------------------------------------------------
void ThermostatFsSolver::iterate_lambda(const MATRIX & rhs)
{
int nNodeOverlap = overlapToNodeMap_->nRows();
DENS_VEC _lambdaOverlap_(nNodeOverlap);
DENS_MAT & lambda(lambda_->set_quantity());
map_unique_to_overlap(lambda,_lambdaOverlap_);
double factor = 0.5*dtFactor_*atc_->dt();
_lambdaOld_.resize(nNodes_,1);
_rhsOverlap_.resize(nNodeOverlap,1);
map_unique_to_overlap(rhs,_rhsOverlap_);
_rhsTotal_.resize(nNodeOverlap);
// solve assuming we get initial guess for lambda
double error(-1.);
for (int i = 0; i < lambdaMaxIterations_; ++i) {
_lambdaOld_ = lambda;
// solve the system with the new rhs
const DENS_MAT & rhsLambdaSquared(rhsLambdaSquared_->quantity());
for (int i = 0; i < nNodeOverlap; i++) {
if (rhsMap_(i,0) == 1.) {
_rhsTotal_(i) = _rhsOverlap_(i,0) + factor*rhsLambdaSquared(i,0);
}
else {
_rhsTotal_(i) = 0.;
}
}
matrixSolver_->execute(_rhsTotal_,_lambdaOverlap_);
// check convergence
map_overlap_to_unique(_lambdaOverlap_,lambda);
lambda_->force_reset();
DENS_MAT difference = lambda-_lambdaOld_;
error = difference.col_norm()/_lambdaOld_.col_norm();
if (error < tolerance_)
break;
}
if (error >= tolerance_) {
stringstream message;
message << "WARNING: Iterative solve for lambda failed to converge after " << lambdaMaxIterations_ << " iterations, final tolerance was " << error << "\n";
ATC::LammpsInterface::instance()->print_msg(message.str());
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatGlcFs
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatGlcFs::ThermostatGlcFs(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
RegulatorMethod(thermostat,regulatorPrefix),
lambdaSolver_(NULL),
mdMassMatrix_(atc_->set_mass_mat_md(TEMPERATURE)),
atomVelocities_(NULL),
temperature_(atc_->field(TEMPERATURE)),
timeFilter_(atomicRegulator_->time_filter()),
nodalAtomicLambdaPower_(NULL),
lambdaPowerFiltered_(NULL),
atomLambdas_(NULL),
atomThermostatForces_(NULL),
atomMasses_(NULL),
isFirstTimestep_(true),
nodalAtomicEnergy_(NULL),
atomPredictedVelocities_(NULL),
nodalAtomicPredictedEnergy_(NULL),
firstHalfAtomForces_(NULL)
{
// construct/obtain data corresponding to stage 3 of ATC_Method::initialize
nodalAtomicLambdaPower_ = thermostat->regulator_data(regulatorPrefix_+"NodalAtomicLambdaPower",1);
lambdaPowerFiltered_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaPowerFiltered",1);
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatGlcFs::construct_transfers()
{
lambdaSolver_->construct_transfers();
InterscaleManager & interscaleManager(atc_->interscale_manager());
// get atom velocity data from manager
atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
// construct lambda evaluated at atom locations
atomLambdas_ = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy");
if (!atomLambdas_) {
DENS_MAN * lambda = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1);
atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
lambda,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
}
// get data from manager
atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
nodalAtomicEnergy_ = interscaleManager.dense_matrix("NodalAtomicEnergy");
// thermostat forces based on lambda and the atomic velocities
atomThermostatForces_ = new AtomicThermostatForce(atc_,atomLambdas_);
interscaleManager.add_per_atom_quantity(atomThermostatForces_,
regulatorPrefix_+"AtomThermostatForce");
// predicted temperature quantities: atom velocities, atom energies, and restricted atom energies
atomPredictedVelocities_ = new AtcAtomQuantity<double>(atc_,atc_->nsd());
// MAKE THINGS WORK WITH ONLY ONE PREDICTED VELOCITY, CHECK IT EXISTS
interscaleManager.add_per_atom_quantity(atomPredictedVelocities_,
regulatorPrefix_+"AtomicPredictedVelocities");
AtomicEnergyForTemperature * atomPredictedEnergyForTemperature = new TwiceKineticEnergy(atc_,
atomPredictedVelocities_);
interscaleManager.add_per_atom_quantity(atomPredictedEnergyForTemperature,
regulatorPrefix_+"AtomicPredictedTwiceKineticEnergy");
nodalAtomicPredictedEnergy_ = new AtfShapeFunctionRestriction(atc_,
atomPredictedEnergyForTemperature,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_dense_matrix(nodalAtomicPredictedEnergy_,
regulatorPrefix_+"NodalAtomicPredictedEnergy");
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void ThermostatGlcFs::initialize()
{
RegulatorMethod::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
lambdaSolver_->initialize(); // ensures initial lambda force is zero
lambdaSolver_->set_lambda_to_value(0.);
*nodalAtomicLambdaPower_ = 0.; // energy change due to thermostats
*lambdaPowerFiltered_ = 0.; // filtered energy change due to thermostats
}
else {
lambdaSolver_->initialize();
// 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(nodalAtomicLambdaPower_->quantity());
}
atomThermostatForces_->quantity(); // initialize
atomThermostatForces_->fix_quantity();
firstHalfAtomForces_ = atomThermostatForces_; // initialize
#ifdef OBSOLETE
compute_rhs_map();
#endif
}
//--------------------------------------------------------
// reset_atom_materials:
// resets the localized atom to material map
//--------------------------------------------------------
void ThermostatGlcFs::reset_atom_materials(const Array<int> & elementToMaterialMap,
const MatrixDependencyManager<DenseMatrix, int> * atomElement)
{
lambdaSolver_->reset_atom_materials(elementToMaterialMap,
atomElement);
}
//--------------------------------------------------------
// apply_to_atoms:
// determines what if any contributions to the
// atomic moition is needed for
// consistency with the thermostat
// and computes the instantaneous induced power
//--------------------------------------------------------
void ThermostatGlcFs::apply_to_atoms(PerAtomQuantity<double> * atomicVelocity,
const DENS_MAN * nodalAtomicEnergy,
const DENS_MAT & lambdaForce,
DENS_MAT & nodalAtomicLambdaPower,
double dt)
{
// compute initial contributions to lambda power
nodalAtomicLambdaPower = nodalAtomicEnergy->quantity();
nodalAtomicLambdaPower *= -1.;
// apply lambda force to atoms
_velocityDelta_ = lambdaForce;
_velocityDelta_ /= atomMasses_->quantity();
_velocityDelta_ *= dt;
(*atomicVelocity) += _velocityDelta_;
// finalize lambda power
nodalAtomicLambdaPower += nodalAtomicEnergy->quantity();
}
//--------------------------------------------------------
// full_prediction:
// flag to perform a full prediction calcalation
// for lambda rather than using the old value
//--------------------------------------------------------
bool ThermostatGlcFs::full_prediction()
{
if (isFirstTimestep_ || ((atc_->atom_to_element_map_type() == EULERIAN)
&& (atc_->atom_to_element_map_frequency() > 1)
&& (atc_->step() % atc_->atom_to_element_map_frequency() == 0 ))) {
return true;
}
return false;
}
//--------------------------------------------------------
// apply_predictor:
// apply the thermostat to the atoms in the first step
// of the Verlet algorithm
//--------------------------------------------------------
void ThermostatGlcFs::apply_pre_predictor(double dt)
{
DENS_MAT & myLambdaPowerFiltered(lambdaPowerFiltered_->set_quantity());
DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());
// update filtered power
timeFilter_->apply_pre_step1(myLambdaPowerFiltered,myNodalAtomicLambdaPower,dt); // equivalent update to measure power as change in energy due to thermostat
// apply lambda force to atoms and compute instantaneous lambda power for first half of time step
this->apply_to_atoms(atomVelocities_,nodalAtomicEnergy_,
firstHalfAtomForces_->quantity(),
myNodalAtomicLambdaPower,0.5*dt);
// update nodal variables for first half of time step
this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy1_,0.5*dt);
// start update of filtered lambda power
myNodalAtomicLambdaPower = 0.; // temporary power for first part of update
timeFilter_->apply_post_step1(myLambdaPowerFiltered,myNodalAtomicLambdaPower,dt);
}
//--------------------------------------------------------
// apply_pre_corrector:
// apply the thermostat to the atoms in the first part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatGlcFs::apply_pre_corrector(double dt)
{
(*atomPredictedVelocities_) = atomVelocities_->quantity();
// do full prediction if we just redid the shape functions
if (full_prediction()) {
this->compute_lambda(dt);
atomThermostatForces_->unfix_quantity(); // allow update of atomic force applied by lambda
}
// apply lambda force to atoms and compute instantaneous lambda power to predict second half of time step
DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());
apply_to_atoms(atomPredictedVelocities_,
nodalAtomicPredictedEnergy_,
firstHalfAtomForces_->quantity(),
myNodalAtomicLambdaPower,0.5*dt);
if (full_prediction())
atomThermostatForces_->fix_quantity();
// update predicted nodal variables for second half of time step
this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy2_,0.5*dt);
// following manipulations performed this way for efficiency
deltaEnergy1_ += deltaEnergy2_;
atc_->apply_inverse_mass_matrix(deltaEnergy1_,TEMPERATURE);
temperature_ += deltaEnergy1_;
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the thermostat to the atoms in the second part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatGlcFs::apply_post_corrector(double dt)
{
// remove predicted power effects
DENS_MAT & myTemperature(temperature_.set_quantity());
atc_->apply_inverse_mass_matrix(deltaEnergy2_,TEMPERATURE);
myTemperature -= deltaEnergy2_;
// set up equation and update lambda
this->compute_lambda(dt);
// apply lambda force to atoms and compute instantaneous lambda power for second half of time step
DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());
// allow computation of force applied by lambda using current velocities
atomThermostatForces_->unfix_quantity();
atomThermostatForces_->quantity();
atomThermostatForces_->fix_quantity();
apply_to_atoms(atomVelocities_,nodalAtomicEnergy_,
atomThermostatForces_->quantity(),
myNodalAtomicLambdaPower,0.5*dt);
// finalize filtered lambda power by adding latest contribution
timeFilter_->apply_post_step2(lambdaPowerFiltered_->set_quantity(),
myNodalAtomicLambdaPower,dt);
// update nodal variables for second half of time step
this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy2_,0.5*dt);
atc_->apply_inverse_mass_matrix(deltaEnergy2_,TEMPERATURE);
myTemperature += deltaEnergy2_;
isFirstTimestep_ = false;
}
//--------------------------------------------------------
// compute_lambda:
// sets up and solves linear system for lambda, if the
// bool is true it iterators to a non-linear solution
//--------------------------------------------------------
void ThermostatGlcFs::compute_lambda(double dt,
bool iterateSolution)
{
// set up rhs for lambda equation
this->set_thermostat_rhs(rhs_,0.5*dt);
// solve system
lambdaSolver_->compute_lambda(rhs_,iterateSolution);
#ifdef OBSOLETE
// solve linear system for lambda guess
DENS_MAT & lambda(lambda_->set_quantity());
solve_for_lambda(rhs_,lambda);
// iterate to solution
if (iterateSolution) {
iterate_lambda(rhs_);
}
#endif
}
//--------------------------------------------------------
// compute_boundary_flux
// default computation of boundary flux based on
// finite
//--------------------------------------------------------
void ThermostatGlcFs::compute_boundary_flux(FIELDS & fields)
{
lambdaSolver_->compute_boundary_flux(fields);
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void ThermostatGlcFs::output(OUTPUT_LIST & outputData)
{
_lambdaPowerOutput_ = nodalAtomicLambdaPower_->quantity();
// approximate value for lambda power
double dt = LammpsInterface::instance()->dt();
_lambdaPowerOutput_ *= (2./dt);
DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"LambdaEnergy"] = &lambda;
outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatSolverFlux
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatSolverFlux::ThermostatSolverFlux(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
ThermostatFsSolver(thermostat,lambdaMaxIterations,regulatorPrefix)
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatSolverFlux::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// set up data for linear solver
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixEnergy");
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
ThermostatFsSolver::construct_transfers();
// add transfers for computation of extra RHS term accounting of O(lambda^2)
// lambda squared followed by fractional step RHS contribution
atomLambdas_ = (interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy"));
if (!atomLambdas_) {
atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
lambda_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
}
LambdaSquared * lambdaSquared = new LambdaSquared(atc_,
interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
weights_,
atomLambdas_);
interscaleManager.add_per_atom_quantity(lambdaSquared,
regulatorPrefix_+"LambdaSquaredMapped");
rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_);
interscaleManager.add_dense_matrix(rhsLambdaSquared_,
regulatorPrefix_+"RhsLambdaSquared");
}
//--------------------------------------------------------
// construct_regulated_nodes:
// constructs the set of nodes being regulated
//--------------------------------------------------------
void ThermostatSolverFlux::construct_regulated_nodes()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// matrix requires all entries even if localized for correct lumping
regulatedNodes_ = interscaleManager.set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
if (!regulatedNodes_) {
regulatedNodes_ = new RegulatedNodes(atc_);
interscaleManager.add_set_int(regulatedNodes_,
regulatorPrefix_+"ThermostatRegulatedNodes");
}
// if localized monitor nodes with applied fluxes
if (atomicRegulator_->use_localized_lambda()) {
if ((atomicRegulator_->coupling_mode() == Thermostat::FLUX) && (atomicRegulator_->boundary_integration_type() != NO_QUADRATURE)) {
// include boundary nodes
applicationNodes_ = new FluxBoundaryNodes(atc_);
boundaryNodes_ = new BoundaryNodes(atc_);
interscaleManager.add_set_int(boundaryNodes_,
regulatorPrefix_+"ThermostatBoundaryNodes");
}
else {
// fluxed nodes only
applicationNodes_ = new FluxNodes(atc_);
}
interscaleManager.add_set_int(applicationNodes_,
regulatorPrefix_+"ThermostatApplicationNodes");
}
else {
applicationNodes_ = regulatedNodes_;
}
// special set of boundary elements for boundary flux quadrature
if ((atomicRegulator_->boundary_integration_type() == FE_INTERPOLATION)
&& (atomicRegulator_->use_localized_lambda())) {
elementMask_ = interscaleManager.dense_matrix_bool(regulatorPrefix_+"BoundaryElementMask");
if (!elementMask_) {
elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
interscaleManager.add_dense_matrix_bool(elementMask_,
regulatorPrefix_+"BoundaryElementMask");
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatIntegratorFlux
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatIntegratorFlux::ThermostatIntegratorFlux(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
ThermostatGlcFs(thermostat,lambdaMaxIterations,regulatorPrefix),
heatSource_(atc_->atomic_source(TEMPERATURE))
{
lambdaSolver_ = new ThermostatSolverFlux(thermostat,
lambdaMaxIterations,
regulatorPrefix);
}
//--------------------------------------------------------
// add_to_temperature
// add in contributions from lambda power and boundary
// flux to the FE temperature
//--------------------------------------------------------
void ThermostatIntegratorFlux::add_to_energy(const DENS_MAT & nodalLambdaPower,
DENS_MAT & deltaEnergy,
double dt)
{
deltaEnergy.resize(nNodes_,1);
const DENS_MAT & myBoundaryFlux(boundaryFlux_[TEMPERATURE].quantity());
for (int i = 0; i < nNodes_; i++) {
deltaEnergy(i,0) = nodalLambdaPower(i,0) + dt*myBoundaryFlux(i,0);
}
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void ThermostatIntegratorFlux::initialize()
{
ThermostatGlcFs::initialize();
// timestep factor
lambdaSolver_->set_timestep_factor(1.);
}
//--------------------------------------------------------
// set_thermostat_rhs:
// sets up the right-hand side including boundary
// fluxes (coupling & prescribed), heat sources, and
// fixed (uncoupled) nodes
//--------------------------------------------------------
void ThermostatIntegratorFlux::set_thermostat_rhs(DENS_MAT & rhs,
double dt)
{
// only tested with flux != 0 + ess bc = 0
// (a) for flux based :
// form rhs : 2/3kB * W_I^-1 * \int N_I r dV
// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
// fluxes are set in ATC transfer
const DENS_MAT & heatSource(heatSource_.quantity());
#if true
const set<int> & applicationNodes((lambdaSolver_->application_nodes())->quantity());
rhs.resize(nNodes_,1);
for (int i = 0; i < nNodes_; i++) {
if (applicationNodes.find(i) != applicationNodes.end()) {
rhs(i,0) = heatSource(i,0);
}
else {
rhs(i,0) = 0.;
}
}
#else
rhs.resize(nNodes_,1);
for (int i = 0; i < nNodes_; i++) {
rhs(i,0) = heatSource(i,0);
}
#endif
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatSolverFixed
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatSolverFixed::ThermostatSolverFixed(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
ThermostatFsSolver(thermostat,lambdaMaxIterations,regulatorPrefix)
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatSolverFixed::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_+"LambdaCouplingMatrixEnergy");
linearSolverType_ = AtomicRegulator::CG_SOLVE;
// base class transfers, e.g. weights
ThermostatFsSolver::construct_transfers();
// add transfers for computation of extra RHS term accounting of O(lambda^2)
// lambda squared followed by fractional step RHS contribution
atomLambdas_ = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy");
if (!atomLambdas_) {
atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
lambda_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
}
if (lambdaAtomMap_) {
LambdaSquaredMapped * lambdaSquared = new LambdaSquaredMapped(atc_,
lambdaAtomMap_,
interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
weights_,
atomLambdas_);
interscaleManager.add_per_atom_quantity(lambdaSquared,
regulatorPrefix_+"LambdaSquared");
rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_);
}
else {
LambdaSquared * lambdaSquared = new LambdaSquared(atc_,
interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
weights_,
atomLambdas_);
interscaleManager.add_per_atom_quantity(lambdaSquared,
regulatorPrefix_+"LambdaSquaredMapped");
rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_);
}
interscaleManager.add_dense_matrix(rhsLambdaSquared_,
regulatorPrefix_+"RhsLambdaSquared");
}
//--------------------------------------------------------
// construct_regulated_nodes:
// constructs the set of nodes being regulated
//--------------------------------------------------------
void ThermostatSolverFixed::construct_regulated_nodes()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
regulatedNodes_ = interscaleManager.set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
if (!regulatedNodes_) {
if (!atomicRegulator_->use_localized_lambda()) {
regulatedNodes_ = new RegulatedNodes(atc_);
}
else if (atomicRegulator_->coupling_mode() == AtomicRegulator::FLUX) {
regulatedNodes_ = new FixedNodes(atc_);
}
else if (atomicRegulator_->coupling_mode() == AtomicRegulator::FIXED) {
// include boundary nodes
regulatedNodes_ = new FixedBoundaryNodes(atc_);
}
else {
throw ATC_Error("ThermostatSolverFixed::construct_regulated_nodes - couldn't determine set of regulated nodes");
}
interscaleManager.add_set_int(regulatedNodes_,
regulatorPrefix_+"ThermostatRegulatedNodes");
}
applicationNodes_ = regulatedNodes_;
// special set of boundary elements for defining regulated atoms
if (atomicRegulator_->use_localized_lambda()) {
elementMask_ = interscaleManager.dense_matrix_bool(regulatorPrefix_+"BoundaryElementMask");
if (!elementMask_) {
elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
interscaleManager.add_dense_matrix_bool(elementMask_,
regulatorPrefix_+"BoundaryElementMask");
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatIntegratorFixed
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatIntegratorFixed::ThermostatIntegratorFixed(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
ThermostatGlcFs(thermostat,lambdaMaxIterations,regulatorPrefix),
atomThermostatForcesPredVel_(NULL),
filterCoefficient_(1.)
{
lambdaSolver_ = new ThermostatSolverFixed(thermostat,
lambdaMaxIterations,
regulatorPrefix);
}
//--------------------------------------------------------
// construct_regulated_nodes:
// constructs the set of nodes being regulated
//--------------------------------------------------------
void ThermostatIntegratorFixed::construct_transfers()
{
ThermostatGlcFs::construct_transfers();
InterscaleManager & interscaleManager(atc_->interscale_manager());
// predicted forces for halving update
atomThermostatForcesPredVel_ = new AtomicThermostatForce(atc_,atomLambdas_,atomPredictedVelocities_);
interscaleManager.add_per_atom_quantity(atomThermostatForcesPredVel_,
regulatorPrefix_+"AtomThermostatForcePredictedVelocity");
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void ThermostatIntegratorFixed::initialize()
{
ThermostatGlcFs::initialize();
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set KE multiplier
AtomicEnergyForTemperature * atomEnergyForTemperature =
static_cast<AtomicEnergyForTemperature * >(interscaleManager.per_atom_quantity("AtomicEnergyForTemperature"));
keMultiplier_ = atomEnergyForTemperature->kinetic_energy_multiplier();
// reset data to zero
deltaFeEnergy_.reset(nNodes_,1);
deltaNodalAtomicEnergy_.reset(nNodes_,1);
// initialize filtered energy
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (!timeFilterManager->end_equilibrate()) {
nodalAtomicEnergyFiltered_ = nodalAtomicEnergy_->quantity();
}
// timestep factor
lambdaSolver_->set_timestep_factor(0.5);
}
//--------------------------------------------------------
// halve_force:
// flag to halve the lambda force for improved
// accuracy
//--------------------------------------------------------
bool ThermostatIntegratorFixed::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;
}
//--------------------------------------------------------
// initialize_delta_nodal_atomic_energy:
// initializes storage for the variable tracking
// the change in the nodal atomic energy
// that has occured over the past timestep
//--------------------------------------------------------
void ThermostatIntegratorFixed::initialize_delta_nodal_atomic_energy(double dt)
{
// initialize delta energy
const DENS_MAT & myNodalAtomicEnergy(nodalAtomicEnergy_->quantity());
initialNodalAtomicEnergy_ = myNodalAtomicEnergy;
initialNodalAtomicEnergy_ *= -1.; // initially stored as negative for efficiency
timeFilter_->apply_pre_step1(nodalAtomicEnergyFiltered_.set_quantity(),
myNodalAtomicEnergy,dt);
}
//--------------------------------------------------------
// compute_delta_nodal_atomic_energy:
// computes the change in the nodal atomic energy
// that has occured over the past timestep
//--------------------------------------------------------
void ThermostatIntegratorFixed::compute_delta_nodal_atomic_energy(double dt)
{
// set delta energy based on predicted atomic velocities
const DENS_MAT & myNodalAtomicEnergy(nodalAtomicEnergy_->quantity());
timeFilter_->apply_post_step1(nodalAtomicEnergyFiltered_.set_quantity(),
myNodalAtomicEnergy,dt);
deltaNodalAtomicEnergy_ = initialNodalAtomicEnergy_;
deltaNodalAtomicEnergy_ += myNodalAtomicEnergy;
}
//--------------------------------------------------------
// compute_lambda:
// sets up and solves linear system for lambda
//--------------------------------------------------------
void ThermostatIntegratorFixed::compute_lambda(double dt,
bool iterateSolution)
{
// compute predicted changes in nodal atomic energy
compute_delta_nodal_atomic_energy(dt);
// change in finite element energy
deltaFeEnergy_ = initialFeEnergy_;
deltaFeEnergy_ += (mdMassMatrix_.quantity())*(temperature_.quantity());
ThermostatGlcFs::compute_lambda(dt,iterateSolution);
}
//--------------------------------------------------------
// apply_predictor:
// apply the thermostat to the atoms in the first step
// of the Verlet algorithm
//--------------------------------------------------------
void ThermostatIntegratorFixed::apply_pre_predictor(double dt)
{
// initialize values to be track change in finite element energy over the timestep
initialize_delta_nodal_atomic_energy(dt);
initialFeEnergy_ = -1.*((mdMassMatrix_.quantity())*(temperature_.quantity())); // initially stored as negative for efficiency
ThermostatGlcFs::apply_pre_predictor(dt);
}
//--------------------------------------------------------
// apply_pre_corrector:
// apply the thermostat to the atoms in the first part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatIntegratorFixed::apply_pre_corrector(double dt)
{
// do full prediction if we just redid the shape functions
if (full_prediction()) {
firstHalfAtomForces_ = atomThermostatForces_; // reset in case this time step needed special treatment
_tempNodalAtomicEnergyFiltered_ = nodalAtomicEnergyFiltered_.quantity();
}
ThermostatGlcFs::apply_pre_corrector(dt);
if (full_prediction()) {
// reset temporary variables
nodalAtomicEnergyFiltered_ = _tempNodalAtomicEnergyFiltered_;
}
if (halve_force()) {
// save old velocities if we are doing halving calculation of lambda force
// copy velocities over into temporary storage
(*atomPredictedVelocities_) = atomVelocities_->quantity();
}
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the thermostat to the atoms in the second part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatIntegratorFixed::apply_post_corrector(double dt)
{
bool halveForce = halve_force();
ThermostatGlcFs::apply_post_corrector(dt);
// update filtered energy with lambda power
DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());
timeFilter_->apply_post_step2(nodalAtomicEnergyFiltered_.set_quantity(),
myNodalAtomicLambdaPower,dt);
if (halveForce) {
// 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
lambdaSolver_->scale_lambda(0.5);
firstHalfAtomForces_ = atomThermostatForcesPredVel_;
atomThermostatForcesPredVel_->unfix_quantity();
}
else {
firstHalfAtomForces_ = atomThermostatForces_;
}
}
//--------------------------------------------------------
// add_to_temperature
// add in contributions from lambda power and boundary
// flux to the FE temperature
//--------------------------------------------------------
void ThermostatIntegratorFixed::add_to_energy(const DENS_MAT & nodalLambdaPower,
DENS_MAT & deltaEnergy,
double dt)
{
deltaEnergy.resize(nNodes_,1);
SetDependencyManager<int> * myRegulatedNodes =
(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
const set<int> & regulatedNodes(myRegulatedNodes->quantity());
for (int i = 0; i < nNodes_; i++) {
if (regulatedNodes.find(i) != regulatedNodes.end()) {
deltaEnergy(i,0) = 0.;
}
else {
deltaEnergy(i,0) = nodalLambdaPower(i,0);
}
}
}
//--------------------------------------------------------
// set_thermostat_rhs:
// sets up the right-hand side including boundary
// fluxes (coupling & prescribed), heat sources, and
// fixed (uncoupled) nodes
//--------------------------------------------------------
void ThermostatIntegratorFixed::set_thermostat_rhs(DENS_MAT & rhs,
double dt)
{
// for essential bcs (fixed nodes) :
// form rhs : (delThetaV - delTheta)/dt
SetDependencyManager<int> * myRegulatedNodes =
(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
const set<int> & regulatedNodes(myRegulatedNodes->quantity());
double factor = (1./dt)/keMultiplier_;
rhs.resize(nNodes_,1);
for (int i = 0; i < nNodes_; i++) {
if (regulatedNodes.find(i) != regulatedNodes.end()) {
rhs(i,0) = factor*(deltaNodalAtomicEnergy_(i,0) - deltaFeEnergy_(i,0));
}
else {
rhs(i,0) = 0.;
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatIntegratorFluxFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatIntegratorFluxFiltered::ThermostatIntegratorFluxFiltered(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
ThermostatIntegratorFlux(thermostat,lambdaMaxIterations,regulatorPrefix)
{
// do nothing
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void ThermostatIntegratorFluxFiltered::initialize()
{
ThermostatIntegratorFlux::initialize();
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (!timeFilterManager->end_equilibrate()) {
// always must start as zero because of filtering scheme
heatSourceOld_.reset(nNodes_,1);
instantHeatSource_.reset(nNodes_,1);
timeStepSource_.reset(nNodes_,1);
}
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the thermostat to the atoms in the second part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatIntegratorFluxFiltered::apply_post_corrector(double dt)
{
// compute lambda
ThermostatIntegratorFlux::apply_post_corrector(dt);
// store data needed for filter inversion of heat flux for thermostat rhs
instantHeatSource_ = rhs_;
heatSourceOld_ = heatSource_.quantity();
}
//--------------------------------------------------------
// add_to_temperature
// add in contributions from lambda power and boundary
// flux to the FE temperature
//--------------------------------------------------------
void ThermostatIntegratorFluxFiltered::add_to_energy(const DENS_MAT & nodalLambdaPower,
DENS_MAT & deltaEnergy,
double dt)
{
deltaEnergy.reset(nNodes_,1);
double coef = timeFilter_->unfiltered_coefficient_post_s1(2.*dt);
const DENS_MAT & myBoundaryFlux(boundaryFlux_[TEMPERATURE].quantity());
for (int i = 0; i < nNodes_; i++) {
deltaEnergy(i,0) = coef*nodalLambdaPower(i,0) + dt*myBoundaryFlux(i,0);
}
}
//--------------------------------------------------------
// set_thermostat_rhs:
// sets up the right-hand side including boundary
// fluxes (coupling & prescribed), heat sources, and
// fixed (uncoupled) nodes
//--------------------------------------------------------
void ThermostatIntegratorFluxFiltered::set_thermostat_rhs(DENS_MAT & rhs,
double dt)
{
// only tested with flux != 0 + ess bc = 0
// (a) for flux based :
// form rhs : 2/3kB * W_I^-1 * \int N_I r dV
// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
// fluxes are set in ATC transfer
// invert heatSource_ to get unfiltered source
// relevant coefficients from time filter
double coefF1 = timeFilter_->filtered_coefficient_pre_s1(2.*dt);
double coefF2 = timeFilter_->filtered_coefficient_post_s1(2.*dt);
double coefU1 = timeFilter_->unfiltered_coefficient_pre_s1(2.*dt);
double coefU2 = timeFilter_->unfiltered_coefficient_post_s1(2.*dt);
const DENS_MAT & heatSource(heatSource_.quantity());
SetDependencyManager<int> * myApplicationNodes =
(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatApplicationNodes");
const set<int> & applicationNodes(myApplicationNodes->quantity());
rhs.resize(nNodes_,1);
for (int i = 0; i < nNodes_; i++) {
if (applicationNodes.find(i) != applicationNodes.end()) {
rhs(i,0) = heatSource(i,0) - coefF1*coefF2*heatSourceOld_(i,0) - coefU1*coefF2*instantHeatSource_(i,0);
rhs(i,0) /= coefU2;
}
else {
rhs(i,0) = 0.;
}
}
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void ThermostatIntegratorFluxFiltered::output(OUTPUT_LIST & outputData)
{
_lambdaPowerOutput_ = lambdaPowerFiltered_->quantity();
// approximate value for lambda power
double dt = LammpsInterface::instance()->dt();
_lambdaPowerOutput_ *= (2./dt);
DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = &lambda;
outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatIntegratorFixedFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatIntegratorFixedFiltered::ThermostatIntegratorFixedFiltered(AtomicRegulator * thermostat,
int lambdaMaxIterations,
const string & regulatorPrefix) :
ThermostatIntegratorFixed(thermostat,lambdaMaxIterations,regulatorPrefix)
{
// do nothing
}
//--------------------------------------------------------
// initialize_delta_nodal_atomic_energy:
// initializes storage for the variable tracking
// the change in the nodal atomic energy
// that has occured over the past timestep
//--------------------------------------------------------
void ThermostatIntegratorFixedFiltered::initialize_delta_nodal_atomic_energy(double dt)
{
// initialize delta energy
DENS_MAT & myNodalAtomicEnergyFiltered(nodalAtomicEnergyFiltered_.set_quantity());
initialNodalAtomicEnergy_ = myNodalAtomicEnergyFiltered;
initialNodalAtomicEnergy_ *= -1.; // initially stored as negative for efficiency
timeFilter_->apply_pre_step1(myNodalAtomicEnergyFiltered,
nodalAtomicEnergy_->quantity(),dt);
}
//--------------------------------------------------------
// compute_delta_nodal_atomic_energy:
// computes the change in the nodal atomic energy
// that has occured over the past timestep
//--------------------------------------------------------
void ThermostatIntegratorFixedFiltered::compute_delta_nodal_atomic_energy(double dt)
{
// set delta energy based on predicted atomic velocities
DENS_MAT & myNodalAtomicEnergyFiltered(nodalAtomicEnergyFiltered_.set_quantity());
timeFilter_->apply_post_step1(myNodalAtomicEnergyFiltered,
nodalAtomicEnergy_->quantity(),dt);
deltaNodalAtomicEnergy_ = initialNodalAtomicEnergy_;
deltaNodalAtomicEnergy_ += myNodalAtomicEnergyFiltered;
}
//--------------------------------------------------------
// add_to_temperature
// add in contributions from lambda power and boundary
// flux to the FE temperature
//--------------------------------------------------------
void ThermostatIntegratorFixedFiltered::add_to_energy(const DENS_MAT & nodalLambdaPower,
DENS_MAT & deltaEnergy,
double dt)
{
deltaEnergy.resize(nNodes_,1);
SetDependencyManager<int> * myRegulatedNodes =
(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
const set<int> & regulatedNodes(myRegulatedNodes->quantity());
double coef = timeFilter_->unfiltered_coefficient_post_s1(2.*dt);
for (int i = 0; i < nNodes_; i++) {
if (regulatedNodes.find(i) != regulatedNodes.end()) {
deltaEnergy(i,0) = 0.;
}
else {
deltaEnergy(i,0) = coef*nodalLambdaPower(i,0);
}
}
}
//--------------------------------------------------------
// set_thermostat_rhs:
// sets up the right-hand side for fixed
// (coupling & prescribed) temperature values
//--------------------------------------------------------
void ThermostatIntegratorFixedFiltered::set_thermostat_rhs(DENS_MAT & rhs,
double dt)
{
// (b) for essential bcs (fixed nodes):
// form rhs : (delThetaV - delTheta)/dt
SetDependencyManager<int> * myRegulatedNodes =
(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
const set<int> & regulatedNodes(myRegulatedNodes->quantity());
double factor = (1./dt)/keMultiplier_;
factor /= timeFilter_->unfiltered_coefficient_post_s1(2.*dt);
rhs.resize(nNodes_,1);
for (int i = 0; i < nNodes_; i++) {
if (regulatedNodes.find(i) != regulatedNodes.end()) {
rhs(i,0) = factor*(deltaNodalAtomicEnergy_(i,0) - deltaFeEnergy_(i,0));
}
else {
rhs(i,0) = 0.;
}
}
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void ThermostatIntegratorFixedFiltered::output(OUTPUT_LIST & outputData)
{
_lambdaPowerOutput_ = lambdaPowerFiltered_->quantity();
// approximate value for lambda power
double dt = LammpsInterface::instance()->dt();
_lambdaPowerOutput_ *= (2./dt);
DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData[regulatorPrefix_+"Lambda"] = &lambda;
outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatFluxFixed
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatFluxFixed::ThermostatFluxFixed(AtomicRegulator * thermostat,
int lambdaMaxIterations,
bool constructThermostats) :
RegulatorMethod(thermostat),
thermostatFlux_(NULL),
thermostatFixed_(NULL),
thermostatBcs_(NULL)
{
if (constructThermostats) {
thermostatFlux_ = new ThermostatIntegratorFlux(thermostat,lambdaMaxIterations,regulatorPrefix_+"Flux");
thermostatFixed_ = new ThermostatIntegratorFixed(thermostat,lambdaMaxIterations,regulatorPrefix_+"Fixed");
// need to choose BC type based on coupling mode
if (thermostat->coupling_mode() == AtomicRegulator::FLUX) {
thermostatBcs_ = thermostatFlux_;
}
else if (thermostat->coupling_mode() == AtomicRegulator::FIXED) {
thermostatBcs_ = thermostatFixed_;
}
else {
throw ATC_Error("ThermostatFluxFixed:create_thermostats - invalid thermostat type provided");
}
}
}
//--------------------------------------------------------
// Destructor
//--------------------------------------------------------
ThermostatFluxFixed::~ThermostatFluxFixed()
{
if (thermostatFlux_) delete thermostatFlux_;
if (thermostatFixed_) delete thermostatFixed_;
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatFluxFixed::construct_transfers()
{
thermostatFlux_->construct_transfers();
thermostatFixed_->construct_transfers();
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void ThermostatFluxFixed::initialize()
{
thermostatFixed_->initialize();
thermostatFlux_->initialize();
}
//--------------------------------------------------------
// apply_predictor:
// apply the thermostat to the atoms in the first step
// of the Verlet algorithm
//--------------------------------------------------------
void ThermostatFluxFixed::apply_pre_predictor(double dt)
{
thermostatFixed_->apply_pre_predictor(dt);
thermostatFlux_->apply_pre_predictor(dt);
}
//--------------------------------------------------------
// apply_pre_corrector:
// apply the thermostat to the atoms in the first part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatFluxFixed::apply_pre_corrector(double dt)
{
thermostatFlux_->apply_pre_corrector(dt);
if (thermostatFixed_->full_prediction()) {
atc_->set_fixed_nodes();
}
thermostatFixed_->apply_pre_corrector(dt);
}
//--------------------------------------------------------
// apply_post_corrector:
// apply the thermostat to the atoms in the second part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatFluxFixed::apply_post_corrector(double dt)
{
thermostatFlux_->apply_post_corrector(dt);
atc_->set_fixed_nodes();
thermostatFixed_->apply_post_corrector(dt);
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void ThermostatFluxFixed::output(OUTPUT_LIST & outputData)
{
thermostatFlux_->output(outputData);
thermostatFixed_->output(outputData);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatFluxFixedFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatFluxFixedFiltered::ThermostatFluxFixedFiltered(AtomicRegulator * thermostat,
int lambdaMaxIterations) :
ThermostatFluxFixed(thermostat,lambdaMaxIterations,false)
{
thermostatFlux_ = new ThermostatIntegratorFluxFiltered(thermostat,lambdaMaxIterations,regulatorPrefix_+"Flux");
thermostatFixed_ = new ThermostatIntegratorFixedFiltered(thermostat,lambdaMaxIterations,regulatorPrefix_+"Fixed");
// need to choose BC type based on coupling mode
if (thermostat->coupling_mode() == AtomicRegulator::FLUX) {
thermostatBcs_ = thermostatFlux_;
}
else if (thermostat->coupling_mode() == AtomicRegulator::FIXED) {
thermostatBcs_ = thermostatFixed_;
}
else {
throw ATC_Error("ThermostatFluxFixed:create_thermostats - invalid thermostat type provided");
}
}
//--------------------------------------------------------
// Class ThermostatGlc
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ThermostatGlc::ThermostatGlc(AtomicRegulator * thermostat) :
ThermostatShapeFunction(thermostat),
timeFilter_(atomicRegulator_->time_filter()),
lambdaPowerFiltered_(NULL),
atomThermostatForces_(NULL),
prescribedDataMgr_(atc_->prescribed_data_manager()),
atomMasses_(NULL)
{
// consistent with stage 3 of ATC_Method::initialize
lambdaPowerFiltered_= atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaPowerFiltered",1);
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatGlc::construct_transfers()
{
ThermostatShapeFunction::construct_transfers();
InterscaleManager & interscaleManager(atc_->interscale_manager());
// get data from manager
atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
// thermostat forces based on lambda and the atomic velocities
AtomicThermostatForce * atomThermostatForces = new AtomicThermostatForce(atc_);
interscaleManager.add_per_atom_quantity(atomThermostatForces,
regulatorPrefix_+"AtomThermostatForce");
atomThermostatForces_ = atomThermostatForces;
}
//--------------------------------------------------------
// apply_to_atoms:
// determines what if any contributions to the
// atomic moition is needed for
// consistency with the thermostat
//--------------------------------------------------------
void ThermostatGlc::apply_to_atoms(PerAtomQuantity<double> * atomVelocities,
const DENS_MAT & lambdaForce,
double dt)
{
_velocityDelta_ = lambdaForce;
_velocityDelta_ /= atomMasses_->quantity();
_velocityDelta_ *= dt;
(*atomVelocities) += _velocityDelta_;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatPowerVerlet
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatPowerVerlet::ThermostatPowerVerlet(AtomicRegulator * thermostat) :
ThermostatGlc(thermostat),
nodalTemperatureRoc_(atc_->field_roc(TEMPERATURE)),
heatSource_(atc_->atomic_source(TEMPERATURE)),
nodalAtomicPower_(NULL),
nodalAtomicLambdaPower_(NULL)
{
// do nothing
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatPowerVerlet::construct_transfers()
{
InterscaleManager & interscaleManager(atc_->interscale_manager());
// set up node mappings
create_node_maps();
// determine if mapping is needed and set up if so
if (atomicRegulator_->use_localized_lambda()) {
lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
regulatorPrefix_+"LambdaAtomMap");
}
// set up linear solver
if (atomicRegulator_->use_lumped_lambda_solve()) {
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
linearSolverType_ = AtomicRegulator::RSL_SOLVE;
}
else {
if (lambdaAtomMap_) {
shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
lambdaAtomMap_,
nodeToOverlapMap_);
}
else {
shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
}
linearSolverType_ = AtomicRegulator::CG_SOLVE;
}
interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
regulatorPrefix_+"LambdaCouplingMatrixEnergy");
// base class transfers, e.g. weights
ThermostatGlc::construct_transfers();
// get managed data
nodalAtomicPower_ = interscaleManager.dense_matrix("NodalAtomicPower");
// power induced by lambda
DotTwiceKineticEnergy * atomicLambdaPower =
new DotTwiceKineticEnergy(atc_,atomThermostatForces_);
interscaleManager.add_per_atom_quantity(atomicLambdaPower,
regulatorPrefix_+"AtomicLambdaPower");
// restriction to nodes of power induced by lambda
nodalAtomicLambdaPower_ = new AtfShapeFunctionRestriction(atc_,
atomicLambdaPower,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_dense_matrix(nodalAtomicLambdaPower_,
regulatorPrefix_+"NodalAtomicLambdaPower");
}
//--------------------------------------------------------
// initialize
// initializes all method data
//--------------------------------------------------------
void ThermostatPowerVerlet::initialize()
{
ThermostatGlc::initialize();
// sets up time filter for cases where variables temporally filtered
TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
if (!timeFilterManager->end_equilibrate()) {
_nodalAtomicLambdaPowerOut_ = 0.;
*lambdaPowerFiltered_ = 0.;
timeFilter_->initialize(lambdaPowerFiltered_->quantity());
}
}
//--------------------------------------------------------
// apply_predictor:
// apply the thermostat to the atoms in the first step
// of the Verlet algorithm
//--------------------------------------------------------
void ThermostatPowerVerlet::apply_pre_predictor(double dt)
{
atomThermostatForces_->unfix_quantity();
compute_thermostat(0.5*dt);
// apply lambda force to atoms
const DENS_MAT & thermostatForces(atomThermostatForces_->quantity());
atomThermostatForces_->fix_quantity();
apply_to_atoms(atomVelocities_,thermostatForces,0.5*dt);
}
//--------------------------------------------------------
// apply_pre_corrector:
// apply the thermostat to the atoms in the first part
// of the corrector step of the Verlet algorithm
//--------------------------------------------------------
void ThermostatPowerVerlet::apply_pre_corrector(double dt)
{
atomThermostatForces_->unfix_quantity();
compute_thermostat(0.5*dt);
// apply lambda force to atoms
const DENS_MAT & thermostatForces(atomThermostatForces_->quantity());
atomThermostatForces_->fix_quantity();
apply_to_atoms(atomVelocities_,thermostatForces,0.5*dt);
}
//--------------------------------------------------------
// add_to_rhs:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the thermostat
//--------------------------------------------------------
void ThermostatPowerVerlet::add_to_rhs(FIELDS & rhs)
{
rhs[TEMPERATURE] += nodalAtomicLambdaPower_->quantity() + boundaryFlux_[TEMPERATURE].quantity();
}
//--------------------------------------------------------
// set_thermostat_rhs:
// sets up the right-hand side including boundary
// fluxes (coupling & prescribed), heat sources, and
// fixed (uncoupled) nodes
//--------------------------------------------------------
void ThermostatPowerVerlet::set_thermostat_rhs(DENS_MAT & rhs_nodes)
{
// (a) for flux based :
// form rhs : \int N_I r dV
// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
// fluxes are set in ATC transfer
rhs_nodes = heatSource_.quantity();
// (b) for ess. bcs
// form rhs : {sum_a (2 * N_Ia * v_ia * f_ia) - (dtheta/dt)_I}
// replace rhs for prescribed nodes
const DENS_MAT & myNodalAtomicPower(nodalAtomicPower_->quantity());
const DIAG_MAT & myMdMassMatrix(mdMassMatrix_.quantity());
const DENS_MAT & myNodalTemperatureRoc(nodalTemperatureRoc_.quantity());
for (int i = 0; i < nNodes_; i++) {
if (prescribedDataMgr_->is_fixed(i,TEMPERATURE,0)) {
rhs_nodes(i,0) = 0.5*(myNodalAtomicPower(i,0) - myMdMassMatrix(i,i)*myNodalTemperatureRoc(i,0));
}
}
}
//--------------------------------------------------------
// compute_thermostat:
// sets up and solves the thermostat equations since
// they are the same at different parts of the time
// step
//--------------------------------------------------------
void ThermostatPowerVerlet::compute_thermostat(double dt)
{
// set up rhs
set_thermostat_rhs(_rhs_);
// solve linear system for lambda
DENS_MAT & myLambda(lambda_->set_quantity());
solve_for_lambda(_rhs_,myLambda);
nodalAtomicLambdaPower_->unfix_quantity(); // enable computation of force applied by lambda
timeFilter_->apply_pre_step1(lambdaPowerFiltered_->set_quantity(),
nodalAtomicLambdaPower_->quantity(),dt);
nodalAtomicLambdaPower_->fix_quantity();
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void ThermostatPowerVerlet::output(OUTPUT_LIST & outputData)
{
_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity();
DENS_MAT & lambda(lambda_->set_quantity());
if ((atc_->lammps_interface())->rank_zero()) {
outputData["lambda"] = &lambda;
outputData["nodalLambdaPower"] = &(_nodalAtomicLambdaPowerOut_);
}
}
//--------------------------------------------------------
// finish:
// final tasks after a run
//--------------------------------------------------------
void ThermostatPowerVerlet::finish()
{
_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatHooverVerlet
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatHooverVerlet::ThermostatHooverVerlet(AtomicRegulator * thermostat) :
ThermostatPowerVerlet(thermostat),
lambdaHoover_(NULL),
nodalAtomicHooverLambdaPower_(NULL)
{
// set up data consistent with stage 3 of ATC_Method::initialize
lambdaHoover_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaHoover",1);
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatHooverVerlet::construct_transfers()
{
ThermostatPowerVerlet::construct_transfers();
InterscaleManager & interscaleManager(atc_->interscale_manager());
FtaShapeFunctionProlongation * atomHooverLambdas = new FtaShapeFunctionProlongation(atc_,
lambdaHoover_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_per_atom_quantity(atomHooverLambdas,
regulatorPrefix_+"AtomHooverLambda");
AtomicThermostatForce * atomHooverThermostatForces = new AtomicThermostatForce(atc_,atomHooverLambdas);
interscaleManager.add_per_atom_quantity(atomHooverThermostatForces,
regulatorPrefix_+"AtomHooverThermostatForce");
SummedAtomicQuantity<double> * atomTotalThermostatForces =
new SummedAtomicQuantity<double>(atc_,atomThermostatForces_,atomHooverThermostatForces);
interscaleManager.add_per_atom_quantity(atomTotalThermostatForces,
regulatorPrefix_+"AtomTotalThermostatForce");
atomThermostatForces_ = atomTotalThermostatForces;
// transfers dependent on time integration method
DotTwiceKineticEnergy * atomicHooverLambdaPower =
new DotTwiceKineticEnergy(atc_,atomHooverThermostatForces);
interscaleManager.add_per_atom_quantity(atomicHooverLambdaPower,
regulatorPrefix_+"AtomicHooverLambdaPower");
nodalAtomicHooverLambdaPower_ = new AtfShapeFunctionRestriction(atc_,
atomicHooverLambdaPower,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_dense_matrix(nodalAtomicHooverLambdaPower_,
regulatorPrefix_+"NodalAtomicHooverLambdaPower");
}
//--------------------------------------------------------
// add_to_rhs:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the thermostat
//--------------------------------------------------------
void ThermostatHooverVerlet::add_to_rhs(FIELDS & rhs)
{
rhs[TEMPERATURE] += _nodalAtomicLambdaPowerOut_;
}
//--------------------------------------------------------
// compute_thermostat:
// sets up and solves the thermostat equations since
// they are the same at different parts of the time
// step
//--------------------------------------------------------
void ThermostatHooverVerlet::compute_thermostat(double dt)
{
// apply prescribed/extrinsic sources and fixed nodes
ThermostatPowerVerlet::compute_thermostat(0.5*dt);
_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); // save power from lambda in power-based thermostat
// set up Hoover rhs
set_hoover_rhs(_rhs_);
// solve linear system for lambda
DENS_MAT & myLambda(lambdaHoover_->set_quantity());
solve_for_lambda(_rhs_,myLambda);
// compute force applied by lambda
// compute nodal atomic power from Hoover coupling
// only add in contribution to uncoupled nodes
if (atomicRegulator_->use_localized_lambda())
add_to_lambda_power(atomThermostatForces_->quantity(),0.5*dt);
}
//--------------------------------------------------------
// set_hoover_rhs:
// sets up the right-hand side for fixed value,
// i.e. Hoover coupling
//--------------------------------------------------------
void ThermostatHooverVerlet::set_hoover_rhs(DENS_MAT & rhs)
{
// form rhs : sum_a ( N_Ia * v_ia * f_ia) - 0.5*M_MD*(dtheta/dt)_I
rhs = nodalAtomicPower_->quantity();
rhs -= mdMassMatrix_.quantity()*nodalTemperatureRoc_.quantity();
rhs /= 2.;
}
//--------------------------------------------------------
// add_to_nodal_lambda_power:
// determines the power exerted by the Hoover
// thermostat at each FE node
//--------------------------------------------------------
void ThermostatHooverVerlet::add_to_lambda_power(const DENS_MAT & myLambdaForce,
double dt)
{
_myNodalLambdaPower_ = nodalAtomicHooverLambdaPower_->quantity();
const INT_ARRAY & nodeToOverlapMap(nodeToOverlapMap_->quantity());
for (int i = 0; i < nNodes_; ++i) {
if (nodeToOverlapMap(i,0)==-1)
_nodalAtomicLambdaPowerOut_(i,0) += _myNodalLambdaPower_(i,0);
else
_myNodalLambdaPower_(i,0) = 0.;
}
timeFilter_->apply_post_step1(lambdaPowerFiltered_->set_quantity(),_myNodalLambdaPower_,dt);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatPowerVerletFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatPowerVerletFiltered::ThermostatPowerVerletFiltered(AtomicRegulator * thermostat) :
ThermostatPowerVerlet(thermostat),
nodalTemperature2Roc_(atc_->field_2roc(TEMPERATURE)),
fieldsRoc_(atc_->fields_roc()),
filterScale_((atc_->time_filter_manager())->filter_scale())
{
heatSourceRoc_.reset(nNodes_,1);
fluxRoc_[TEMPERATURE].reset(nNodes_,1);
}
//--------------------------------------------------------
// compute_boundary_flux
// also sets time derivatives of boundary flux and
// heat sources
//--------------------------------------------------------
void ThermostatPowerVerletFiltered::compute_boundary_flux(FIELDS & fields)
{
ThermostatPowerVerlet::compute_boundary_flux(fields);
// compute boundary flux rate of change
fluxRoc_[TEMPERATURE] = 0.;
atc_->compute_boundary_flux(fieldMask_,
fieldsRoc_,
fluxRoc_,
atomMaterialGroups_,
shpFcnDerivs_);
// compute extrinsic model rate of change
(atc_->extrinsic_model_manager()).set_sources(fieldsRoc_,fluxRoc_);
heatSourceRoc_ = fluxRoc_[TEMPERATURE].quantity();
}
//--------------------------------------------------------
// add_to_rhs:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the thermostat
//--------------------------------------------------------
void ThermostatPowerVerletFiltered::add_to_rhs(FIELDS & rhs)
{
rhs[TEMPERATURE] += lambdaPowerFiltered_->quantity() + boundaryFlux_[TEMPERATURE].quantity();
}
//--------------------------------------------------------
// set_thermostat_rhs:
// sets up the right-hand side including boundary
// fluxes (coupling & prescribed), heat sources, and
// fixed (uncoupled) nodes
//--------------------------------------------------------
void ThermostatPowerVerletFiltered::set_thermostat_rhs(DENS_MAT & rhs_nodes)
{
// (a) for flux based :
// form rhs : \int N_I r dV
// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
// fluxes are set in ATC transfer
rhs_nodes = heatSource_.quantity() + filterScale_*heatSourceRoc_.quantity();
// (b) for ess. bcs
// form rhs : {sum_a (N_Ia * v_ia * f_ia) - 0.5*(dtheta/dt)_I}
const DENS_MAT & myNodalAtomicPower(nodalAtomicPower_->quantity());
const DIAG_MAT & myMdMassMatrix(mdMassMatrix_.quantity());
const DENS_MAT & myNodalTemperatureRoc(nodalTemperatureRoc_.quantity());
const DENS_MAT & myNodalTemperature2Roc(nodalTemperature2Roc_.quantity());
for (int i = 0; i < nNodes_; i++) {
if (prescribedDataMgr_->is_fixed(i,TEMPERATURE,0)) {
rhs_nodes(i,0) = 0.5*(myNodalAtomicPower(i,0) - myMdMassMatrix(i,i)*(myNodalTemperatureRoc(i,0)+ filterScale_*myNodalTemperature2Roc(i,0)));
}
}
}
//--------------------------------------------------------
// output:
// adds all relevant output to outputData
//--------------------------------------------------------
void ThermostatPowerVerletFiltered::output(OUTPUT_LIST & outputData)
{
outputData["lambda"] = &(lambda_->set_quantity());
outputData["nodalLambdaPower"] = &(lambdaPowerFiltered_->set_quantity());
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ThermostatHooverVerletFiltered
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
// Grab references to ATC and thermostat data
//--------------------------------------------------------
ThermostatHooverVerletFiltered::ThermostatHooverVerletFiltered(AtomicRegulator * thermostat) :
ThermostatPowerVerletFiltered(thermostat),
lambdaHoover_(NULL),
nodalAtomicHooverLambdaPower_(NULL)
{
// consistent with stage 3 of ATC_Method::initialize
lambdaHoover_ = atomicRegulator_->regulator_data("LambdaHoover",1);
}
//--------------------------------------------------------
// constructor_transfers
// instantiates or obtains all dependency managed data
//--------------------------------------------------------
void ThermostatHooverVerletFiltered::construct_transfers()
{
ThermostatPowerVerletFiltered::construct_transfers();
InterscaleManager & interscaleManager(atc_->interscale_manager());
FtaShapeFunctionProlongation * atomHooverLambdas = new FtaShapeFunctionProlongation(atc_,
lambdaHoover_,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_per_atom_quantity(atomHooverLambdas,
regulatorPrefix_+"AtomHooverLambda");
AtomicThermostatForce * atomHooverThermostatForces = new AtomicThermostatForce(atc_,atomHooverLambdas);
interscaleManager.add_per_atom_quantity(atomHooverThermostatForces,
regulatorPrefix_+"AtomHooverThermostatForce");
SummedAtomicQuantity<double> * atomTotalThermostatForces =
new SummedAtomicQuantity<double>(atc_,atomThermostatForces_,atomHooverThermostatForces);
interscaleManager.add_per_atom_quantity(atomTotalThermostatForces,
regulatorPrefix_+"AtomTotalThermostatForce");
atomThermostatForces_ = atomTotalThermostatForces;
// transfers dependent on time integration method
DotTwiceKineticEnergy * atomicHooverLambdaPower =
new DotTwiceKineticEnergy(atc_,atomHooverThermostatForces);
interscaleManager.add_per_atom_quantity(atomicHooverLambdaPower,
regulatorPrefix_+"AtomicHooverLambdaPower");
nodalAtomicHooverLambdaPower_ = new AtfShapeFunctionRestriction(atc_,
atomicHooverLambdaPower,
interscaleManager.per_atom_sparse_matrix("Interpolant"));
interscaleManager.add_dense_matrix(nodalAtomicHooverLambdaPower_,
regulatorPrefix_+"NodalAtomicHooverLambdaPower");
}
//--------------------------------------------------------
// compute_thermostat:
// sets up and solves the thermostat equations since
// they are the same at different parts of the time
// step
//--------------------------------------------------------
void ThermostatHooverVerletFiltered::compute_thermostat(double dt)
{
// apply prescribed/extrinsic sources and fixed nodes
ThermostatPowerVerletFiltered::compute_thermostat(0.5*dt);
_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); // save power from lambda in power-based thermostat
// set up Hoover rhs
set_hoover_rhs(_rhs_);
// solve linear system for lambda
DENS_MAT & myLambda(lambdaHoover_->set_quantity());
solve_for_lambda(_rhs_,myLambda);
// compute force applied by lambda
// compute nodal atomic power from Hoover coupling
// only add in contribution to uncoupled nodes
if (atomicRegulator_->use_localized_lambda())
add_to_lambda_power(atomThermostatForces_->quantity(),0.5*dt);
}
//--------------------------------------------------------
// set_hoover_rhs:
// sets up the right-hand side for fixed value,
// i.e. Hoover coupling
//--------------------------------------------------------
void ThermostatHooverVerletFiltered::set_hoover_rhs(DENS_MAT & rhs)
{
// form rhs : sum_a (N_Ia * v_ia * f_ia) - 0.5*M_MD*(dtheta/dt)_I
rhs = nodalAtomicPower_->quantity();
rhs -= mdMassMatrix_.quantity()*(nodalTemperatureRoc_.quantity() + filterScale_*nodalTemperature2Roc_.quantity());
rhs /= 2.;
}
//--------------------------------------------------------
// add_to_nodal_lambda_power:
// determines the power exerted by the Hoover
// thermostat at each FE node
//--------------------------------------------------------
void ThermostatHooverVerletFiltered::add_to_lambda_power(const DENS_MAT & myLambdaForce,
double dt)
{
_myNodalLambdaPower_ = nodalAtomicHooverLambdaPower_->quantity();
const INT_ARRAY nodeToOverlapMap(nodeToOverlapMap_->quantity());
for (int i = 0; i < nNodes_; ++i) {
if (nodeToOverlapMap(i,0)==-1)
_nodalAtomicLambdaPowerOut_(i,0) += _myNodalLambdaPower_(i,0);
else
_myNodalLambdaPower_(i,0) = 0.;
}
timeFilter_->apply_post_step1(lambdaPowerFiltered_->set_quantity(),_myNodalLambdaPower_,dt);
}
//--------------------------------------------------------
// add_to_rhs:
// determines what if any contributions to the
// finite element equations are needed for
// consistency with the thermostat
//--------------------------------------------------------
void ThermostatHooverVerletFiltered::add_to_rhs(FIELDS & rhs)
{
rhs[TEMPERATURE] += lambdaPowerFiltered_->quantity();
}
};
c4science · Help