Thermostat.cpp
0977af1299aalammps-icms

Thermostat.cpp
View Options

			// ATC_Transfer Headers
			#include "Thermostat.h"
			#include "CG.h"
			#include "ATC_Error.h"
			#include "PrescribedDataManager.h"
			#include "TimeIntegrator.h"

			namespace ATC {

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class Thermostat
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			//--------------------------------------------------------
			Thermostat::Thermostat(ATC_Transfer * atcTransfer) :
			AtomicRegulator(atcTransfer),
			thermostatType_(NONE)
			{
			// 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;

			// thermostat type
			/*! \page man_thermal_control fix_modify AtC transfer thermal control
			\section syntax
			fix_modify AtC transfer thermal control <control_type> <optional args>\n
			- control_type (string) = none \| rescale \| hoover \| flux\n

			fix_modify AtC transfer thermal control rescale <frequency>\n
			- frequency (int) = time step frequency for applying velocity rescaling \n

			fix_modify AtC transfer thermal control hoover \n

			fix_modify AtC transfer thermal control flux <boundary_integration_type> <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 transfer thermal control none </TT> \n
			<TT> fix_modify AtC transfer thermal control rescale 10 </TT> \n
			<TT> fix_modify AtC transfer thermal control hoover </TT> \n
			<TT> fix_modify AtC transfer thermal control flux 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
			*/
			int argIndex = 0;
			if (strcmp(arg[argIndex],"none")==0) { // restore defaults
			thermostatType_ = NONE;
			howOften_ = 1;
			boundaryIntegrationType_ = ATC_Transfer::NO_QUADRATURE;
			foundMatch = true;
			}
			else if (strcmp(arg[argIndex],"rescale")==0) {
			argIndex++;
			howOften_ = atoi(arg[argIndex]);
			if (howOften_ < 1) {
			throw ATC_Error(0,"Bad rescaling thermostat frequency");
			}
			else {
			thermostatType_ = RESCALE;
			boundaryIntegrationType_ = ATC_Transfer::NO_QUADRATURE;
			foundMatch = true;
			}
			}
			else if (strcmp(arg[argIndex],"hoover")==0) {
			thermostatType_ = HOOVER;
			howOften_ = 1;
			boundaryIntegrationType_ = ATC_Transfer::NO_QUADRATURE;
			foundMatch = true;
			}
			else if (strcmp(arg[argIndex],"flux")==0) {
			thermostatType_ = FLUX;
			howOften_ = 1;
			argIndex++;
			boundaryIntegrationType_ = atcTransfer_->parse_boundary_integration(narg-argIndex,&arg[argIndex],boundaryFaceSet_);
			foundMatch = true;
			}

			if (!foundMatch)
			foundMatch = AtomicRegulator::modify(narg,arg);
			if (foundMatch)
			needReset_ = true;
			return foundMatch;
			}

			//--------------------------------------------------------
			// reset_nodal_lambda_power:
			// resets the thermostat generated power to a
			// prescribed value
			//--------------------------------------------------------
			void Thermostat::reset_lambda_power(DENS_MAT & target)
			{
			for (int i = 0; i < nNodes_; ++i)
			lambdaPowerFiltered_(i) = target(i,0);
			}

			//--------------------------------------------------------
			// initialize:
			// sets up methods before a run
			// NOTE we currently only include time filter
			// dependence, but in general there is also a
			// time integrator dependence. In general the
			// precedence order is:
			// time filter -> time integrator -> kinetostat
			// In the future this may need to be added if
			// different types of time integrators can be
			// specified.
			//--------------------------------------------------------
			void Thermostat::initialize()
			{
			//TimeIntegrationType myIntegrationType = (atcTransfer_->get_time_integrator())->get_time_integration_type();
			// HACK until thermal time integrator is implemented
			TimeIntegrator::TimeIntegrationType myIntegrationType = TimeIntegrator::VERLET;
			TimeFilterManager * timeFilterManager = atcTransfer_->get_time_filter_manager();

			// reset data if needed and perform any error/conflict checking
			if (resetData_) {
			AtomicRegulator::reset_data();

			// set up storage
			lambda_.reset(nNodes_,1);
			nodalAtomicLambdaPower_.reset(nNodes_);
			lambdaPowerFiltered_.reset(nNodes_);
			}

			if (needReset_ \|\| timeFilterManager->need_reset() \|\| timeFilterManager->end_equilibrate()) {
			// eliminate existing methods
			destroy();
			if (timeFilterManager->need_reset()) {
			if (timeFilter_)
			delete timeFilter_;
			if (myIntegrationType == TimeIntegrator::VERLET)
			timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT);
			else if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP)
			timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT_IMPLICIT);
			}

			if (timeFilterManager->filter_dynamics()) {
			switch (thermostatType_) {
			case NONE: {
			break;
			}
			case RESCALE: { // error check, rescale and filtering not supported together
			throw ATC_Error(0,"Cannot use rescaling thermostat with time filtering");
			break;
			}
			case HOOVER: {
			if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP)
			//regulatorMethod_ = new ThermostatHooverFiltered(this);
			cout << "temporary line\n";
			else
			regulatorMethod_ = new ThermostatHooverVerletFiltered(this);
			break;
			}
			case FLUX: {
			if (atcTransfer_->use_localized_lambda())
			throw ATC_Error(0,"Cannot use flux thermostating with localized lambda");
			if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP)
			regulatorMethod_ = new ThermostatPowerFiltered(this);
			else
			regulatorMethod_ = new ThermostatPowerVerletFiltered(this);
			break;
			}
			default:
			throw ATC_Error(0,"Unknown thermostat type in Thermostat::initialize");
			}
			}
			else {
			switch (thermostatType_) {
			case NONE: {
			break;
			}
			case RESCALE: {
			regulatorMethod_ = new ThermostatRescale(this);
			break;
			}
			case HOOVER: {
			if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP)
			regulatorMethod_ = new ThermostatHoover(this);
			else
			regulatorMethod_ = new ThermostatHooverVerlet(this);
			break;
			}
			case FLUX: {
			if (atcTransfer_->use_localized_lambda())
			throw ATC_Error(0,"Cannot use flux thermostating with localized lambda");
			if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP)
			regulatorMethod_ = new ThermostatPower(this);
			else
			regulatorMethod_ = new ThermostatPowerVerlet(this);
			break;
			}
			default:
			throw ATC_Error(0,"Unknown thermostat type in Thermostat::initialize");
			}
			}

			AtomicRegulator::reset_method();
			}

			AtomicRegulator::initialize();
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatShapeFunction
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			//--------------------------------------------------------
			ThermostatShapeFunction::ThermostatShapeFunction(Thermostat * thermostat) :
			RegulatorShapeFunction(thermostat),
			thermostat_(thermostat),
			v_(atcTransfer_->get_v()),
			mdMassMatrix_(atcTransfer_->get_mass_mat_md(TEMPERATURE))
			{
			fieldMask_(TEMPERATURE,FLUX) = true;
			}

			//---------------------------------------------------------
			// set_weights:
			// set the diagonal weighting matrix to be the atomic
			// temperatures
			//---------------------------------------------------------
			void ThermostatShapeFunction::set_weights(DIAG_MAT & weights)
			{
			if (nLocalLambda_>0) {
			DENS_VEC weightVector(nLocal_);
			for (int i = 0; i < nLocal_; i++)
			for (int j = 0; j < nsd_; j++)
			weightVector(i) += v_[internalToAtom_(i)][j]*v_[internalToAtom_(i)][j];
			DENS_VEC maskedWeightVector(nLocalLambda_);
			maskedWeightVector = internalToOverlapMap_*weightVector;
			weights.reset(maskedWeightVector);
			}
			}

			//--------------------------------------------------------
			// reset_nlocal:
			// resets data dependent on local atom count
			//--------------------------------------------------------
			void ThermostatShapeFunction::reset_nlocal()
			{
			RegulatorShapeFunction::reset_nlocal();
			if (nLocal_ > 0) {
			atcTransfer_->compute_atomic_mass(atomicMass_);
			}
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatRescale
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			//--------------------------------------------------------
			ThermostatRescale::ThermostatRescale(Thermostat * thermostat) :
			ThermostatShapeFunction(thermostat),
			nodalTemperature_(atcTransfer_->get_field(TEMPERATURE))
			{
			// do nothing
			}

			//--------------------------------------------------------
			// apply_post_corrector:
			// apply the thermostat in the post corrector phase
			//--------------------------------------------------------
			void ThermostatRescale::apply_post_corrector(double dt)
			{
			// compute right-hand side
			DENS_MAT rhs(nNodes_,1);
			rhs = mdMassMatrix_*nodalTemperature_;

			// solve equations
			solve_for_lambda(rhs);

			// application of rescaling lambda due
			apply_to_atoms(v_,dt);
			}

			//--------------------------------------------------------
			// apply_lambda_to_atoms:
			// applies the velocity rescale with an existing lambda
			// note oldAtomicQuantity and dt are not used
			//--------------------------------------------------------
			void ThermostatRescale::apply_to_atoms(double ** atomicQuantity,
			const double dt)
			{
			if (nLocal_>0) {
			DENS_MAT lambdaAtom(nLocal_,1);
			atcTransfer_->prolong(lambda_,lambdaAtom);

			double xi;
			for (int i = 0; i < nLocal_; i++) {
			xi = sqrt(lambdaAtom(i,0)/atomicMass_(i));
			for (int j = 0; j < nsd_; j++) {
			atomicQuantity[internalToAtom_(i)][j] *= xi;
			}
			}
			}
			}

			//--------------------------------------------------------
			// output:
			// adds all relevant output to outputData
			//--------------------------------------------------------
			void ThermostatRescale::output(double dt, OUTPUT_LIST & outputData)
			{
			outputData["lambda"] = &(lambda_);
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatGlc
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			//--------------------------------------------------------
			ThermostatGlc::ThermostatGlc(Thermostat * thermostat) :
			ThermostatShapeFunction(thermostat),
			timeFilter_(atomicRegulator_->get_time_filter()),
			nodalAtomicLambdaPower_(thermostat_->get_nodal_atomic_lambda_power()),
			lambdaPowerFiltered_(thermostat_->get_filtered_lambda_power()),
			lambdaForce_(atomicRegulator_->get_lambda_force()),
			isFixedNode_(atcTransfer_->get_fixed_node_flags(TEMPERATURE))
			{
			// do nothing
			}

			//--------------------------------------------------------
			// apply_to_atoms:
			// determines what if any contributions to the
			// atomic moition is needed for
			// consistency with the thermostat
			//--------------------------------------------------------
			void ThermostatGlc::apply_to_atoms(double ** atomicQuantity,
			const DENS_MAT & lambdaForce,
			double dt)
			{
			if (nLocal_>0) {
			// explicit update
			for (int i = 0; i < nLocal_; i++) {
			for (int j = 0; j < nsd_; j++) {
			atomicQuantity[internalToAtom_(i)][j] += dt*lambdaForce(i,j)/atomicMass_(i);
			}
			}
			}
			}

			//---------------------------------------------------------
			// compute_lambda_force:
			// computes the atomic force applied by lambda
			//---------------------------------------------------------
			void ThermostatGlc::compute_lambda_force(double * const * atomicQuantity,
			DENS_MAT & lambdaForce)
			{
			if (nLocal_>0) {
			// scaled prolongation to (unique) nodes
			DENS_MAT lambdaAtom(nLocal_,1);
			atcTransfer_->prolong(lambda_,lambdaAtom);

			// explicit update
			for (int i = 0; i < nLocal_; i++) {
			for (int j = 0; j < nsd_; j++) {
			lambdaForce(i,j) = -lambdaAtom(i,0)*atomicQuantity[internalToAtom_(i)][j];
			}
			}
			}
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatPower
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			// Grab references to ATC and thermostat data
			//--------------------------------------------------------
			ThermostatPower::ThermostatPower(Thermostat * thermostat) :
			ThermostatGlc(thermostat),
			nodalTemperature_(atcTransfer_->get_field(TEMPERATURE)),
			nodalAtomicTemperature_(atcTransfer_->get_atomic_field(TEMPERATURE)),
			heatSource_(atcTransfer_->get_atomic_source(TEMPERATURE)),
			nodalWeights_(atcTransfer_->get_shape_function_weights()),
			myAtomicVelocity_(NULL),
			f_(atcTransfer_->get_f()),
			nLocalTotal_(0),
			lambdaMaxIterations_(3),
			lambdaTolerance_(1.e-14),
			filterCoefficient_(1.),
			Nstar_(atcTransfer_->get_shape_function_ghost_overlap()),
			ghostToAtom_(atcTransfer_->get_ghost_to_atom_map())
			{
			deltaTemperature_.reset(nNodes_,1);
			deltaNodalAtomicTemperature_.reset(nNodes_,1);
			nodalTemperatureRoc_.reset(nNodes_,1);

			// sets up time filter for cases where variables temporally filtered
			TimeFilterManager * timeFilterManager = atcTransfer_->get_time_filter_manager();
			if (!timeFilterManager->end_equilibrate()) {
			// NOTE do we need a good initial guess for lambdaForce_?
			nodalAtomicLambdaPower_ = 0.;
			lambdaPowerFiltered_ = 0.;
			timeFilter_->initialize(lambdaPowerFiltered_);
			}
			}

			//--------------------------------------------------------
			// Destructor
			// delete thermostat data
			//--------------------------------------------------------
			ThermostatPower::~ThermostatPower()
			{
			destroy();
			}

			//--------------------------------------------------------
			// destroy
			// deletes allocated memory
			//--------------------------------------------------------
			void ThermostatPower::destroy()
			{
			if (myAtomicVelocity_) {
			for (int i = 0; i < nLocalTotal_; i++)
			delete [] myAtomicVelocity_[i];
			delete [] myAtomicVelocity_;
			}
			}

			//--------------------------------------------------------
			// reset_nlocal:
			// resizes myAtomicVelocity force if necessary
			//--------------------------------------------------------
			void ThermostatPower::reset_nlocal()
			{
			destroy();
			ThermostatGlc::reset_nlocal();
			nLocalTotal_ = atcTransfer_->get_nlocal_total();
			if (nLocal_ > 0) {
			myAtomicVelocity_ = new double*[nLocalTotal_];
			for (int i = 0; i < nLocalTotal_; i++)
			myAtomicVelocity_[i] = new double[nsd_];
			}
			}

			//--------------------------------------------------------
			// apply_predictor:
			// apply the thermostat to the atoms in the first step
			// of the Verlet algorithm
			//--------------------------------------------------------
			void ThermostatPower::apply_pre_predictor(double dt)
			{
			// initialize delta temperatures
			deltaTemperature_ = -1.*nodalTemperature_;
			deltaNodalAtomicTemperature_ = -1.*nodalAtomicTemperature_;

			// apply lambda force to atoms and compute instantaneous lambda power
			apply_to_atoms(v_,lambdaForce_,dt);

			// update nodal variables
			update_nodal_quantities_predictor(dt);
			}

			//--------------------------------------------------------
			// update_nodal_quantities_predictor:
			// updates all needed nodal quantities after the
			// predictor step
			//--------------------------------------------------------
			void ThermostatPower::update_nodal_quantities_predictor(double dt)
			{
			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);
			DENS_MAT boundaryFluxRoc(nNodes_,1);
			atcTransfer_->apply_inverse_mass_matrix(boundaryFlux_[TEMPERATURE],
			boundaryFluxRoc,
			TEMPERATURE);
			nodalTemperatureRoc_ += boundaryFluxRoc;

			// update FE temperature
			nodalTemperature_ += dt*nodalTemperatureRoc_;

			// finalize filtered lambda power
			timeFilter_->apply_pre_step1(lambdaPowerFiltered_,nodalAtomicLambdaPower_,dt);
			}

			//--------------------------------------------------------
			// apply_pre_corrector:
			// apply the thermostat to the atoms in the first part
			// of the corrector step of the Verlet algorithm
			//--------------------------------------------------------
			void ThermostatPower::apply_pre_corrector(double dt)
			{
			// copy velocities over into temporary storage
			if (nLocal_ > 0) {
			for (int i = 0; i < nLocalTotal_; i++)
			for (int j = 0; j < nsd_; j++)
			myAtomicVelocity_[i][j] = v_[i][j];
			}

			// compute predicted changes in FE and nodal atomic temperatures
			compute_delta_nodal_atomic_temperature(dt);
			DENS_MAT myDeltaTemperature(deltaTemperature_);
			deltaTemperature_ += nodalTemperature_;

			// set up rhs for lambda equation
			DENS_MAT rhs(nNodes_,1);
			set_thermostat_rhs(rhs,dt);

			// solve linear system for lambda guess
			solve_for_lambda(rhs);

			// iterate to solution
			iterate_lambda(rhs,v_,dt);

			// compute force applied by lambda
			compute_lambda_force(v_,lambdaForce_);

			// apply lambda force to atoms and compute instantaneous lambda power
			apply_to_atoms(myAtomicVelocity_,lambdaForce_,dt);

			// update nodal variables
			update_nodal_quantities_pre_corrector(dt);

			// reset temporary variables
			deltaTemperature_ = myDeltaTemperature;
			}

			//--------------------------------------------------------
			// compute_delta_nodal_atomic_temperature:
			// computes the change in the nodal atomic temperature
			// that has occured over the past timestep
			//--------------------------------------------------------
			void ThermostatPower::compute_delta_nodal_atomic_temperature(double dt)
			{
			// set delta temperatures based on predicted atomic velocities
			DENS_MAT atomicTemperature(nLocal_,1);
			atcTransfer_->compute_atomic_temperature(atomicTemperature, v_);
			DENS_MAT nodalAtomicTemperature(nNodes_,1);
			// NOTE change this when we use physical mass matrices
			atcTransfer_->restrict(atomicTemperature,nodalAtomicTemperature);
			deltaNodalAtomicTemperature_ += nodalAtomicTemperature;
			}

			//--------------------------------------------------------
			// update_nodal_quantities_pre_corrector:
			// updates all needed nodal quantities after the
			// pre-corrector step
			//--------------------------------------------------------
			void ThermostatPower::update_nodal_quantities_pre_corrector(double dt)
			{
			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);
			DENS_MAT boundaryFluxRoc(nNodes_,1);
			atcTransfer_->apply_inverse_mass_matrix(boundaryFlux_[TEMPERATURE],
			boundaryFluxRoc,
			TEMPERATURE);
			nodalTemperatureRoc_ += boundaryFluxRoc;

			// update FE temperature
			nodalTemperature_ += dt*nodalTemperatureRoc_;
			}

			//--------------------------------------------------------
			// undo_update_nodal_quantities_pre_corrector:
			// undoes the nodal quantities update from after the
			// pre-corrector step
			//--------------------------------------------------------
			void ThermostatPower::undo_update_nodal_quantities_pre_corrector(double dt)
			{
			// remove predicted power effects
			nodalTemperature_ -= dt*nodalTemperatureRoc_;
			}

			//--------------------------------------------------------
			// apply_post_corrector:
			// apply the thermostat to the atoms in the second part
			// of the corrector step of the Verlet algorithm
			//--------------------------------------------------------
			void ThermostatPower::apply_post_corrector(double dt)
			{
			// remove predicted power effects
			undo_update_nodal_quantities_pre_corrector(dt);

			// set delta FE temperature
			deltaTemperature_ += nodalTemperature_;

			// set up rhs for lambda equation
			DENS_MAT rhs(nNodes_,1);
			set_thermostat_rhs(rhs,dt);

			// solve linear system for lambda guess
			solve_for_lambda(rhs);

			// iterate to solution
			iterate_lambda(rhs,myAtomicVelocity_,dt);

			// compute force applied by lambda
			compute_lambda_force(myAtomicVelocity_,lambdaForce_);

			// apply lambda force to atoms and compute instantaneous lambda power
			apply_to_atoms(v_,lambdaForce_,dt);

			// update nodal variables
			update_nodal_quantities_post_corrector(dt);
			}

			//--------------------------------------------------------
			// update_nodal_quantities_post_corrector:
			// updates all needed nodal quantities after the
			// post-corrector step
			//--------------------------------------------------------
			void ThermostatPower::update_nodal_quantities_post_corrector(double dt)
			{
			// finalize filtered lambda power
			timeFilter_->apply_post_step1(lambdaPowerFiltered_,nodalAtomicLambdaPower_,dt);

			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);
			DENS_MAT boundaryFluxRoc(nNodes_,1);
			atcTransfer_->apply_inverse_mass_matrix(boundaryFlux_[TEMPERATURE],
			boundaryFluxRoc,
			TEMPERATURE);
			nodalTemperatureRoc_ += boundaryFluxRoc;

			// update FE temperature
			nodalTemperature_ += dt*nodalTemperatureRoc_;
			}

			//--------------------------------------------------------
			// apply_to_atoms:
			// add the thermostat force to the atoms and compute
			// the instantaneous induced power
			//--------------------------------------------------------
			void ThermostatPower::apply_to_atoms(double ** atomicVelocity,
			const DENS_MAT & lambdaForce,
			double dt)
			{
			// compute initial contributions to lambda power
			DENS_VEC atomicTemperature(nLocal_);
			DENS_VEC nodalAtomicTemperature(nNodes_);
			atcTransfer_->compute_atomic_temperature(atomicTemperature, atomicVelocity);
			atcTransfer_->restrict_unscaled(atomicTemperature,nodalAtomicTemperature);
			nodalAtomicLambdaPower_ = -1.*nodalAtomicTemperature;

			ThermostatGlc::apply_to_atoms(atomicVelocity,lambdaForce,dt);

			// finalize lambda power
			atcTransfer_->compute_atomic_temperature(atomicTemperature, atomicVelocity);
			atcTransfer_->restrict_unscaled(atomicTemperature,nodalAtomicTemperature);
			nodalAtomicLambdaPower_ += nodalAtomicTemperature;
			nodalAtomicLambdaPower_ = (1./dt)*nodalAtomicLambdaPower_;
			}

			//---------------------------------------------------------
			// set_weights:
			// set the diagonal weighting matrix to be the atomic
			// temperatures
			//---------------------------------------------------------
			void ThermostatPower::set_weights(DIAG_MAT & weights)
			{
			if (nLocalLambda_>0) {
			DENS_VEC weightVector(nLocal_);
			DENS_VEC maskedWeightVector(nLocalLambda_);
			atcTransfer_->compute_atomic_temperature(weightVector,v_,myAtomicVelocity_);
			maskedWeightVector = internalToOverlapMap_*weightVector;
			weights.reset(maskedWeightVector);
			}
			}

			//--------------------------------------------------------
			// set_thermostat_rhs:
			// sets up the right-hand side including boundary
			// fluxes (coupling & prescribed), heat sources, and
			// fixed (uncoupled) nodes
			//--------------------------------------------------------
			void ThermostatPower::set_thermostat_rhs(DENS_MAT & rhs_nodes,
			double dt)
			{
			// NOTE is the scaling right? (a) E/t/kB & (b) T/t/kB
			// 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
			// NOTE change coef and nodalWeights when physical mass matrices are used
			// and remove rhoCp
			double rhoCp = LammpsInterface::instance()->heat_capacity();
			double k_boltzmann = LammpsInterface::instance()->kBoltzmann();
			double coef = filterCoefficient_rhoCp2./(nsd_*k_boltzmann);
			rhs_nodes = coefheatSource_nodalWeights_;

			// (b) for essential bcs (fixed nodes) :
			// form rhs : (delThetaV - delTheta)/dt
			DENS_MAT rhs_nodes_prescribed(nNodes_,1);
			rhs_nodes_prescribed = (1./dt)filterCoefficient_(deltaNodalAtomicTemperature_-deltaTemperature_);

			// replace rhs for prescribed nodes
			// NOTE can we make is_fixed member data?
			for (int i = 0; i < nNodes_; i++) {
			if (isFixedNode_(i,0)) {
			rhs_nodes(i,0) = rhs_nodes_prescribed(i,0);
			}
			}
			}

			//--------------------------------------------------------
			// iterate_lambda:
			// iteratively solves the equations for lambda
			// for the higher order dt corrections, assuming
			// an initial guess for lambda
			//--------------------------------------------------------
			void ThermostatPower::iterate_lambda(const MATRIX & rhs,
			double * const * atomicVelocity,
			double dt)
			{
			DENS_VEC lambdaOverlap(nNodeOverlap_);
			atcTransfer_->map_unique_to_overlap(lambda_, lambdaOverlap);

			DENS_VEC atomicTemperature;
			DENS_MAT lambdaAtom;
			DENS_MAT lambdaSquaredT;
			if (nLocalLambda_>0) {
			lambdaAtom.reset(nLocal_,1);
			lambdaSquaredT.reset(nLocal_,1);
			atomicTemperature.reset(nLocal_);
			atcTransfer_->compute_atomic_temperature(atomicTemperature, atomicVelocity);
			}

			DENS_MAT lambdaOld(nNodes_,1);
			DENS_MAT rhsOverlap(nNodeOverlap_,1);
			DENS_MAT rhsTotalLocal(nNodes_,1);
			DENS_MAT rhsTotal(nNodes_,1);
			// solve assuming we get initial guess for lambda
			for (int i = 0; i < lambdaMaxIterations_; ++i) {
			lambdaOld = lambda_;

			// assemble quadratic term
			if (nLocalLambda_ > 0) {
			atcTransfer_->prolong_scaled(lambda_, lambdaAtom);
			for (int i = 0; i < nLocal_; i++) {
			lambdaSquaredT(i,0) = lambdaAtom(i,0)lambdaAtom(i,0)atomicTemperature(i);
			}
			rhsOverlap = (dt/4.)shapeFunctionMatrix_.transMat(internalToOverlapMap_lambdaSquaredT);
			atcTransfer_->map_overlap_to_unique(rhsOverlap,rhsTotalLocal);
			}
			LammpsInterface::instance()->allsum(rhsTotalLocal.get_ptr(),rhsTotal.get_ptr(),rhsTotal.size());

			// the system with the new rhs
			rhsTotal += rhs;
			solve_for_lambda(rhsTotal);

			// check convergence
			DENS_MAT difference = lambda_-lambdaOld;
			if (difference.col_norm()/lambdaOld.col_norm() < lambdaTolerance_) break;

			}
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatHoover
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			// Grab references to ATC and thermostat data
			//--------------------------------------------------------
			ThermostatHoover::ThermostatHoover(Thermostat * thermostat) :
			ThermostatPower(thermostat),
			nodeToOverlapMap_(atcTransfer_->get_node_to_overlap_map())
			{
			lambdaForceHoover_.reset(nLocal_,nsd_);
			}

			//--------------------------------------------------------
			// reset_nlocal:
			// resizes Hoover lambda force if necessary
			//--------------------------------------------------------
			void ThermostatHoover::reset_nlocal()
			{
			ThermostatPower::reset_nlocal();
			if (nLocal_ > 0)
			lambdaForceHoover_.reset(nLocal_,nsd_);
			}

			//--------------------------------------------------------
			// apply_predictor:
			// apply the thermostat to the atoms in the first step
			// of the Verlet algorithm
			//--------------------------------------------------------
			void ThermostatHoover::apply_pre_predictor(double dt)
			{
			ThermostatPower::apply_pre_predictor(dt);

			// apply lambda force to atoms and compute instantaneous lambda power
			DENS_VEC myNodalAtomicLambdaPower(nodalAtomicLambdaPower_);
			apply_to_atoms(v_,lambdaForceHoover_,dt);

			// compute nodal atomic power from Hoover coupling
			// only add in contribution to uncoupled nodes
			if (atcTransfer_->use_localized_lambda())
			add_to_lambda_power_predictor(dt);
			else
			nodalAtomicLambdaPower_ = 0.;

			// finish updating lambda power
			nodalAtomicLambdaPower_ += myNodalAtomicLambdaPower;
			}

			//--------------------------------------------------------
			// add_to_lambda_power_predictor:
			// adds Hoover contributions to the lambda power term
			//--------------------------------------------------------
			void ThermostatHoover::add_to_lambda_power_predictor(double dt)
			{
			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);

			// update FE temperature
			for (int i = 0; i < nNodes_; ++i) {
			if (nodeToOverlapMap_(i)==-1)
			nodalTemperature_(i,0) += dt*nodalTemperatureRoc_(i,0);
			else
			nodalAtomicLambdaPower_(i) = 0.;
			}

			// finalize filtered lambda power
			timeFilter_->apply_pre_step2(lambdaPowerFiltered_,nodalAtomicLambdaPower_,dt);
			}

			//--------------------------------------------------------
			// apply_pre_corrector:
			// apply the thermostat to the atoms in the first part
			// of the corrector step of the Verlet algorithm
			//--------------------------------------------------------
			void ThermostatHoover::apply_pre_corrector(double dt)
			{
			ThermostatPower::apply_pre_corrector(dt);

			// set up Hoover rhs
			DENS_MAT myDeltaTemperature(deltaTemperature_);
			deltaTemperature_ += nodalTemperature_;
			DENS_MAT rhs(nNodes_,1);
			set_hoover_rhs(rhs,dt);

			// solve linear system for lambda
			solve_for_lambda(rhs);

			// iterate to solution
			iterate_lambda(rhs,v_,dt);

			// compute force applied by lambda
			compute_lambda_force(v_,lambdaForceHoover_);

			// apply the force to the atoms
			DENS_VEC myNodalAtomicLambdaPower(nodalAtomicLambdaPower_);
			apply_to_atoms(myAtomicVelocity_,lambdaForceHoover_,dt);

			// compute nodal atomic power from Hoover coupling
			// only add in contribution to uncoupled nodes
			if (atcTransfer_->use_localized_lambda())
			add_to_lambda_power_pre_corrector(dt);
			else
			nodalAtomicLambdaPower_ = 0.;

			// finish updating lambda power
			nodalAtomicLambdaPower_ += myNodalAtomicLambdaPower;

			// reset temporary variables
			deltaTemperature_ = myDeltaTemperature;
			}

			//--------------------------------------------------------
			// add_to_lambda_power_pre_corrector:
			// adds Hoover contributions to the lambda power term
			//--------------------------------------------------------
			void ThermostatHoover::add_to_lambda_power_pre_corrector(double dt)
			{
			DENS_MAT myNodalTemperatureRoc(nodalTemperatureRoc_);
			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);

			// update FE temperature
			for (int i = 0; i < nNodes_; ++i) {
			if (nodeToOverlapMap_(i)==-1)
			nodalTemperature_(i,0) += dt*nodalTemperatureRoc_(i,0);
			else {
			nodalTemperatureRoc_(i,0) = 0.;
			nodalAtomicLambdaPower_(i) = 0.;
			}
			}

			// correct nodal temperature rate of change
			nodalTemperatureRoc_ += myNodalTemperatureRoc;
			}

			//--------------------------------------------------------
			// apply_post_corrector:
			// apply the thermostat to the atoms in the second part
			// of the corrector step of the Verlet algorithm
			//--------------------------------------------------------
			void ThermostatHoover::apply_post_corrector(double dt)
			{
			ThermostatPower::apply_post_corrector(dt);

			// set up Hoover rhs
			DENS_MAT rhs(nNodes_,1);
			set_hoover_rhs(rhs,dt);

			// solve linear system for lambda
			solve_for_lambda(rhs);

			// iterate to solution
			iterate_lambda(rhs,myAtomicVelocity_,dt);

			// compute force applied by lambda
			compute_lambda_force(myAtomicVelocity_,lambdaForceHoover_);

			// apply the force to the atoms
			DENS_VEC myNodalAtomicLambdaPower(nodalAtomicLambdaPower_);
			apply_to_atoms(v_,lambdaForceHoover_,dt);

			// compute nodal atomic power from Hoover coupling
			// only add in contribution to uncoupled nodes
			if (atcTransfer_->use_localized_lambda())
			add_to_lambda_power_post_corrector(dt);
			else
			nodalAtomicLambdaPower_ = 0.;

			// finish updating lambda power
			nodalAtomicLambdaPower_ += myNodalAtomicLambdaPower;
			}

			//--------------------------------------------------------
			// add_to_lambda_power_post_corrector:
			// adds Hoover contributions to the lambda power term
			//--------------------------------------------------------
			void ThermostatHoover::add_to_lambda_power_post_corrector(double dt)
			{
			// finalize filtered lambda power
			timeFilter_->apply_post_step2(lambdaPowerFiltered_,nodalAtomicLambdaPower_,dt);

			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);

			// update FE temperature
			for (int i = 0; i < nNodes_; ++i) {
			if (nodeToOverlapMap_(i)==-1)
			nodalTemperature_(i,0) += dt*nodalTemperatureRoc_(i,0);
			else
			nodalAtomicLambdaPower_(i) = 0.;
			}
			}

			//--------------------------------------------------------
			// set_hoover_rhs:
			// sets up the right-hand side for Hoover coupling with
			// boundary nodes
			//--------------------------------------------------------
			void ThermostatHoover::set_hoover_rhs(DENS_MAT & rhs_nodes,
			double dt)
			{
			// form rhs : (delThetaV - delTheta)/dt
			rhs_nodes = (1./dt)filterCoefficient_(deltaNodalAtomicTemperature_-deltaTemperature_);
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatPowerFiltered
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			// Grab references to ATC and thermostat data
			//--------------------------------------------------------
			ThermostatPowerFiltered::ThermostatPowerFiltered(Thermostat * thermostat) :
			ThermostatPower(thermostat)
			{
			// do nothing
			}

			//--------------------------------------------------------
			// update_nodal_quantities_predictor:
			// updates all needed nodal quantities after the
			// predictor step
			//--------------------------------------------------------
			void ThermostatPowerFiltered::update_nodal_quantities_predictor(double dt)
			{
			// update nodal atomic temperature
			//atcTransfer_->apply_inverse_mass_matrix(lambdaPowerFiltered_,
			// nodalTemperatureRoc_,
			// TEMPERATURE);
			// NOTE temporary until new mass matrices are used
			DIAG_MAT & atomicWeights(atcTransfer_->get_atomic_weights());
			nodalTemperatureRoc_ = (1./atomicWeights(0,0))nodalWeights_lambdaPowerFiltered_;
			nodalAtomicTemperature_ += dt*nodalTemperatureRoc_;

			// update FE temperature
			atcTransfer_->apply_inverse_mass_matrix(lambdaPowerFiltered_,
			nodalTemperatureRoc_,
			TEMPERATURE);
			//DENS_MAT boundaryFluxRoc(nNodes_,1);
			//atcTransfer_->apply_inverse_mass_matrix(boundaryFlux_[TEMPERATURE],
			// boundaryFluxRoc,
			// TEMPERATURE);
			//nodalTemperatureRoc_ += boundaryFluxRoc;
			nodalTemperature_ += 0.5dtnodalTemperatureRoc_;

			// finalize filtered lambda power
			timeFilter_->apply_pre_step1(lambdaPowerFiltered_,nodalAtomicLambdaPower_,dt);
			}

			//--------------------------------------------------------
			// compute_delta_nodal_atomic_temperature:
			// computes the change in the nodal atomic temperature
			// that has occured over the past timestep
			//--------------------------------------------------------
			void ThermostatPowerFiltered::compute_delta_nodal_atomic_temperature(double dt)
			{
			// compute next lambda power assuming zero new force
			//DENS_VEC myLambdaPowerFiltered(lambdaPowerFiltered_);
			DENS_VEC zeroNodalLambdaPower(nNodes_);
			timeFilter_->apply_post_step1(lambdaPowerFiltered_,zeroNodalLambdaPower,dt);

			// update nodal atomic temperature
			//atcTransfer_->apply_inverse_mass_matrix(lambdaPowerFiltered_,
			// nodalTemperatureRoc_,
			// TEMPERATURE);
			// NOTE temporary until new mass matrices are used
			DIAG_MAT & atomicWeights(atcTransfer_->get_atomic_weights());
			nodalTemperatureRoc_ = (1./atomicWeights(0,0))nodalWeights_lambdaPowerFiltered_;
			nodalAtomicTemperature_ += dt*nodalTemperatureRoc_;

			// update FE temperature
			atcTransfer_->apply_inverse_mass_matrix(lambdaPowerFiltered_,
			nodalTemperatureRoc_,
			TEMPERATURE);
			nodalTemperature_ += .5dtnodalTemperatureRoc_;

			deltaNodalAtomicTemperature_ += nodalAtomicTemperature_;
			}

			//--------------------------------------------------------
			// update_nodal_quantities_pre_corrector:
			// updates all needed nodal quantities after the
			// pre-corrector step
			//--------------------------------------------------------
			void ThermostatPowerFiltered::update_nodal_quantities_pre_corrector(double dt)
			{
			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);

			double filterCoefficient = timeFilter_->get_unfiltered_coefficient_post_s1(dt);

			//DENS_MAT boundaryFluxRoc(nNodes_,1);
			//atcTransfer_->apply_inverse_mass_matrix(boundaryFlux_[TEMPERATURE],
			// boundaryFluxRoc,
			// TEMPERATURE);
			//nodalTemperatureRoc_ = boundaryFluxRoc + filterCoefficient*nodalTemperatureRoc_;
			nodalTemperatureRoc_ = 0.5filterCoefficientnodalTemperatureRoc_;
			// update FE temperature
			nodalTemperature_ += dt*nodalTemperatureRoc_;
			}

			// //--------------------------------------------------------
			// // undo_update_nodal_quantities_pre_corrector:
			// // undoes the update of all nodal quantities from after
			// // the pre-corrector step
			// //--------------------------------------------------------
			// void ThermostatPowerFiltered::undo_update_nodal_quantities_pre_corrector(double dt)
			// {
			// double filterCoefficient = timeFilter_->get_unfiltered_coefficient_pre_s2(dt);

			// // update FE temperature
			// nodalTemperature_ -= dtfilterCoefficientnodalTemperatureRoc_;
			// }

			//--------------------------------------------------------
			// update_nodal_quantities_post_corrector:
			// updates all needed nodal quantities after the
			// post-corrector step
			//--------------------------------------------------------
			void ThermostatPowerFiltered::update_nodal_quantities_post_corrector(double dt)
			{
			// finalize filtered lambda power
			timeFilter_->apply_post_step2(lambdaPowerFiltered_,nodalAtomicLambdaPower_,dt);
			double filterCoefficient = timeFilter_->get_unfiltered_coefficient_post_s1(dt);

			// update nodal atomic temperature
			//atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			// nodalTemperatureRoc_,
			// TEMPERATURE);
			// NOTE temporary until new mass matrices are used
			DIAG_MAT & atomicWeights(atcTransfer_->get_atomic_weights());
			nodalTemperatureRoc_ = (1./atomicWeights(0,0))nodalWeights_nodalAtomicLambdaPower_;
			//nodalTemperatureRoc_.print("NTR");
			nodalAtomicTemperature_ += filterCoefficientdtnodalTemperatureRoc_;

			// update FE temperature
			atcTransfer_->apply_inverse_mass_matrix(nodalAtomicLambdaPower_,
			nodalTemperatureRoc_,
			TEMPERATURE);
			//DENS_MAT boundaryFluxRoc(nNodes_,1);
			//atcTransfer_->apply_inverse_mass_matrix(boundaryFlux_[TEMPERATURE],
			// boundaryFluxRoc,
			// TEMPERATURE);
			//nodalTemperatureRoc_ = boundaryFluxRoc + filterCoefficient*nodalTemperatureRoc_;
			nodalTemperatureRoc_ = .5filterCoefficientnodalTemperatureRoc_;
			nodalTemperature_ += dt*nodalTemperatureRoc_;
			}

			//--------------------------------------------------------
			// set_thermostat_rhs:
			// sets up the right-hand side including boundary
			// fluxes (coupling & prescribed), heat sources, and
			// fixed (uncoupled) nodes
			//--------------------------------------------------------
			void ThermostatPowerFiltered::set_thermostat_rhs(DENS_MAT & rhs_nodes,
			double dt)
			{
			filterCoefficient_ = 1./(timeFilter_->get_unfiltered_coefficient_post_s1(dt));
			// NOTE change the 0.5 when we change how the lambda power is accounted
			// NOTE change rhoCp^-1 when physical mass matrices are used
			double rhoCp = LammpsInterface::instance()->heat_capacity();
			heatSource_ += (0.5)*lambdaPowerFiltered_;
			ThermostatPower::set_thermostat_rhs(rhs_nodes,dt);
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatPowerVerlet
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			// Grab references to ATC and thermostat data
			//--------------------------------------------------------
			ThermostatPowerVerlet::ThermostatPowerVerlet(Thermostat * thermostat) :
			ThermostatGlc(thermostat),
			nodalTemperatureRoc_(atcTransfer_->get_field_roc(TEMPERATURE)),
			heatSource_(atcTransfer_->get_atomic_source(TEMPERATURE)),
			f_(atcTransfer_->get_f())
			{
			// sets up time filter for cases where variables temporally filtered
			TimeFilterManager * timeFilterManager = atcTransfer_->get_time_filter_manager();
			if (!timeFilterManager->end_equilibrate()) {
			nodalAtomicLambdaPower_ = 0.;
			lambdaPowerFiltered_ = 0.;
			timeFilter_->initialize(lambdaPowerFiltered_);
			}
			}

			//--------------------------------------------------------
			// apply_predictor:
			// apply the thermostat to the atoms in the first step
			// of the Verlet algorithm
			//--------------------------------------------------------
			void ThermostatPowerVerlet::apply_pre_predictor(double dt)
			{
			compute_thermostat(dt);

			// apply lambda force to atoms
			apply_to_atoms(v_,lambdaForce_,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)
			{
			compute_thermostat(dt);

			// apply lambda force to atoms
			apply_to_atoms(v_,lambdaForce_,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_ + boundaryFlux_[TEMPERATURE];
			}

			//--------------------------------------------------------
			// 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_;

			// (b) for ess. bcs
			// form rhs : {sum_a (2 * N_Ia * v_ia * f_ia) - (dtheta/dt)_I}
			DENS_MAT atomicPower;
			atcTransfer_->compute_atomic_power(atomicPower, v_, f_);
			DENS_MAT rhs_nodes_prescribed(nNodes_,1);
			atcTransfer_->restrict_volumetric_quantity(atomicPower, rhs_nodes_prescribed);
			rhs_nodes_prescribed -= 0.5(mdMassMatrix_nodalTemperatureRoc_);

			// replace rhs for prescribed nodes
			for (int i = 0; i < nNodes_; i++) {
			if (isFixedNode_(i,0)) {
			rhs_nodes(i,0) = rhs_nodes_prescribed(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
			DENS_MAT rhs(nNodes_,1);
			set_thermostat_rhs(rhs);

			// solve linear system for lambda
			solve_for_lambda(rhs);

			// compute force applied by lambda
			compute_lambda_force(v_,lambdaForce_);

			// compute nodal atomic power
			compute_nodal_lambda_power(dt);
			}

			//--------------------------------------------------------
			// compute_nodal_lambda_power:
			// determines the power exerted by the thermostat at
			// each FE node
			//--------------------------------------------------------
			void ThermostatPowerVerlet::compute_nodal_lambda_power(double dt)
			{
			DENS_VEC atomicLambdaPower;
			atcTransfer_->compute_atomic_power(atomicLambdaPower,v_,lambdaForce_);
			atcTransfer_->restrict_volumetric_quantity(atomicLambdaPower,nodalAtomicLambdaPower_);
			nodalAtomicLambdaPower_ *= 2.; // accounts for kinetic definition of temperature
			timeFilter_->apply_pre_step1(lambdaPowerFiltered_,nodalAtomicLambdaPower_,dt);
			}

			//--------------------------------------------------------
			// output:
			// adds all relevant output to outputData
			//--------------------------------------------------------
			void ThermostatPowerVerlet::output(double dt, OUTPUT_LIST & outputData)
			{
			outputData["lambda"] = &(lambda_);
			outputData["nodalLambdaPower"] = &(nodalAtomicLambdaPower_);
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatHooverVerlet
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			// Grab references to ATC and thermostat data
			//--------------------------------------------------------
			ThermostatHooverVerlet::ThermostatHooverVerlet(Thermostat * thermostat) :
			ThermostatPowerVerlet(thermostat),
			nodeToOverlapMap_(atcTransfer_->get_node_to_overlap_map())
			{
			// do nothing
			}

			//--------------------------------------------------------
			// 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(dt);

			// set up Hoover rhs
			DENS_MAT rhs(nNodes_);
			set_hoover_rhs(rhs);

			// solve linear system for lambda
			solve_for_lambda(rhs);

			// compute force applied by lambda
			DENS_MAT lambdaForceHoover(nLocal_,nsd_);
			compute_lambda_force(v_,lambdaForceHoover);
			lambdaForce_ += lambdaForceHoover;

			// compute nodal atomic power from Hoover coupling
			// only add in contribution to uncoupled nodes
			if (atcTransfer_->use_localized_lambda())
			add_to_lambda_power(lambdaForceHoover,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_nodes)
			{
			// form rhs : sum_a ( N_Ia * v_ia * f_ia) - 0.5M_MD(dtheta/dt)_I
			DENS_MAT atomicPower(atcTransfer_->get_nlocal(),1);
			atcTransfer_->compute_atomic_power(atomicPower, v_, f_);
			atcTransfer_->restrict_volumetric_quantity(atomicPower, rhs_nodes);
			rhs_nodes -= 0.5(mdMassMatrix_nodalTemperatureRoc_);
			}

			//--------------------------------------------------------
			// add_to_nodal_lambda_power:
			// determines the power exerted by the Hoover
			// thermostat at each FE node
			//--------------------------------------------------------
			void ThermostatHooverVerlet::add_to_lambda_power(DENS_MAT & myLambdaForce,
			double dt)
			{
			DENS_VEC atomicLambdaPower;
			atcTransfer_->compute_atomic_power(atomicLambdaPower,v_,myLambdaForce);
			DENS_VEC myNodalLambdaPower(nNodes_);
			atcTransfer_->restrict_volumetric_quantity(atomicLambdaPower,myNodalLambdaPower);
			myNodalLambdaPower *= 2.;
			for (int i = 0; i < nNodes_; ++i) {
			if (nodeToOverlapMap_(i)==-1)
			nodalAtomicLambdaPower_(i) += myNodalLambdaPower(i);
			else
			myNodalLambdaPower(i) = 0.;
			}
			timeFilter_->apply_post_step1(lambdaPowerFiltered_,myNodalLambdaPower,dt);
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class ThermostatPowerVerletFiltered
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			// Grab references to ATC and thermostat data
			//--------------------------------------------------------
			ThermostatPowerVerletFiltered::ThermostatPowerVerletFiltered(Thermostat * thermostat) :
			ThermostatPowerVerlet(thermostat),
			nodalTemperature2Roc_(atcTransfer_->get_field_2roc(TEMPERATURE)),
			fieldsRoc_(atcTransfer_->get_fields_roc()),
			filterScale_((atcTransfer_->get_time_filter_manager())->get_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.;
			atcTransfer_->compute_boundary_flux(fieldMask_,
			fieldsRoc_,
			fluxRoc_);

			// NOTE add time derivatives from sources, currently not implemented

			// compute extrinsic model rate of change
			(atcTransfer_->get_extrinsic_model_manager()).set_sources(fieldsRoc_,fluxRoc_);
			heatSourceRoc_ = fluxRoc_[TEMPERATURE];


			}

			//--------------------------------------------------------
			// 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(std::map<FieldName,DenseMatrix<double> > & rhs)
			{
			rhs[TEMPERATURE] += lambdaPowerFiltered_ + boundaryFlux_[TEMPERATURE];
			}

			//--------------------------------------------------------
			// 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
			// for (int i = 0; i < nNodes_; i++) {
			// rhs_nodes(i) = heatSource_(i,0) + filterScale_*heatSourceRoc_(i,0);
			// }
			rhs_nodes = heatSource_ + filterScale_*heatSourceRoc_;

			// (b) for ess. bcs
			// form rhs : {sum_a (N_Ia * v_ia * f_ia) - 0.5*(dtheta/dt)_I}
			DENS_MAT atomicPower;
			atcTransfer_->compute_atomic_power(atomicPower, v_, f_);
			DENS_MAT rhs_nodes_prescribed(nNodes_,1);
			atcTransfer_->restrict_volumetric_quantity(atomicPower, rhs_nodes_prescribed);
			// DENS_VEC myTemperatureRoc(nNodes_);
			// for (int i = 0; i < nNodes_; i++) {
			// myTemperatureRoc(i) = nodalTemperatureRoc_(i,0) + filterScale_*nodalTemperature2Roc_(i,0);
			// }
			rhs_nodes_prescribed -= 0.5(mdMassMatrix_(nodalTemperatureRoc_ + filterScale_*nodalTemperature2Roc_));

			// replace rhs for prescribed nodes
			for (int i = 0; i < nNodes_; i++) {
			if (isFixedNode_(i,0)) {
			rhs_nodes(i,0) = rhs_nodes_prescribed(i,0);
			}
			}
			}

			//--------------------------------------------------------
			// output:
			// adds all relevant output to outputData
			//--------------------------------------------------------
			void ThermostatPowerVerletFiltered::output(double dt, OUTPUT_LIST & outputData)
			{
			outputData["lambda"] = &(lambda_);
			outputData["nodalLambdaPower"] = &(lambdaPowerFiltered_);
			}

			//--------------------------------------------------------
			//--------------------------------------------------------
			// Class HooverVerletFiltered
			//--------------------------------------------------------
			//--------------------------------------------------------

			//--------------------------------------------------------
			// Constructor
			// Grab references to ATC and thermostat data
			//--------------------------------------------------------
			ThermostatHooverVerletFiltered::ThermostatHooverVerletFiltered(Thermostat * thermostat) :
			ThermostatPowerVerletFiltered(thermostat),
			nodeToOverlapMap_(atcTransfer_->get_node_to_overlap_map())
			{
			// do nothing
			}

			//--------------------------------------------------------
			// 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(dt);

			// set up Hoover rhs
			DENS_MAT rhs(nNodes_,1);
			set_hoover_rhs(rhs);

			// solve linear system for lambda
			solve_for_lambda(rhs);

			// NOTE modify these sections to be for hoover
			// compute force applied by lambda
			DENS_MAT lambdaForceHoover(nLocal_,nsd_);
			compute_lambda_force(v_,lambdaForceHoover);
			lambdaForce_ += lambdaForceHoover;

			// compute nodal atomic power from Hoover coupling
			// only add in contribution to uncoupled nodes
			if (atcTransfer_->use_localized_lambda())
			add_to_lambda_power(lambdaForceHoover,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_nodes)
			{
			// form rhs : sum_a (N_Ia * v_ia * f_ia) - 0.5M_MD(dtheta/dt)_I
			DENS_MAT atomicPower(atcTransfer_->get_nlocal(),1);
			atcTransfer_->compute_atomic_power(atomicPower, v_, f_);
			atcTransfer_->restrict_volumetric_quantity(atomicPower, rhs_nodes);
			// DENS_VEC myTemperatureRoc(nNodes_);
			// for (int i = 0; i < nNodes_; i++) {
			// myTemperatureRoc(i) = nodalTemperatureRoc_(i,0) + filterScale_*nodalTemperature2Roc_(i,0);
			// }
			rhs_nodes -= 0.5(mdMassMatrix_(nodalTemperatureRoc_ + filterScale_*nodalTemperature2Roc_));
			}

			//--------------------------------------------------------
			// add_to_nodal_lambda_power:
			// determines the power exerted by the Hoover
			// thermostat at each FE node
			//--------------------------------------------------------
			void ThermostatHooverVerletFiltered::add_to_lambda_power(DENS_MAT & myLambdaForce,
			double dt)
			{
			DENS_VEC atomicLambdaPower;
			atcTransfer_->compute_atomic_power(atomicLambdaPower,v_,myLambdaForce);
			DENS_VEC myNodalLambdaPower(nNodes_);
			atcTransfer_->restrict_volumetric_quantity(atomicLambdaPower,myNodalLambdaPower);
			myNodalLambdaPower *= 2.;
			for (int i = 0; i < nNodes_; ++i) {
			if (nodeToOverlapMap_(i)==-1)
			nodalAtomicLambdaPower_(i) += myNodalLambdaPower(i);
			else
			myNodalLambdaPower(i) = 0.;
			}
			timeFilter_->apply_post_step1(lambdaPowerFiltered_,myNodalLambdaPower,dt);
			}

			};

lammps/lib/atc/Thermostat.cpp0977af1299aalammps-icms

Thermostat.cppView Options

lammps/lib/atc/Thermostat.cpp
0977af1299aalammps-icms

Thermostat.cpp
View Options