ATC version 2.0, date: Aug22

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10658 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
rjones 2013-08-23 00:14:36 +00:00
parent 9f64c9b12b
commit dddbef699d
17 changed files with 292 additions and 991 deletions

View File

@ -34,6 +34,9 @@ namespace ATC {
elementMask_(NULL),
elementMaskMass_(NULL),
elementMaskMassMd_(NULL),
nodalAtomicMass_(NULL),
nodalAtomicCount_(NULL),
nodalAtomicHeatCapacity_(NULL),
internalToMask_(NULL),
internalElement_(NULL),
ghostElement_(NULL),
@ -732,6 +735,10 @@ namespace ATC {
}
}
//--------------------------------------------------------
// construct_methods
// have managers instantiate requested algorithms
// and methods
//--------------------------------------------------------
void ATC_Coupling::construct_methods()
{
@ -748,6 +755,11 @@ namespace ATC {
massMatTimeFilters_[thisField] = timeFilterManager_.construct(TimeFilterManager::INSTANTANEOUS);
}
}
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_methods();
}
atomicRegulator_->construct_methods();
}
//-------------------------------------------------------------------
void ATC_Coupling::init_filter()
@ -833,8 +845,8 @@ namespace ATC {
field != fields.end(); field++) {
FieldName thisFieldName = field->first;
FIELDS::const_iterator fieldItr = fields.find(thisFieldName);
const DENS_MAT & field = (fieldItr->second).quantity();
atomicSources[thisFieldName].reset(field.nRows(),field.nCols());
const DENS_MAT & f = (fieldItr->second).quantity();
atomicSources[thisFieldName].reset(f.nRows(),f.nCols());
}
}
}
@ -981,8 +993,8 @@ namespace ATC {
field != fields.end(); field++) {
FieldName thisFieldName = field->first;
FIELDS::const_iterator fieldItr = fields.find(thisFieldName);
const DENS_MAT & field = (fieldItr->second).quantity();
(rhs[thisFieldName].set_quantity()).reset(field.nRows(),field.nCols());
const DENS_MAT & f = (fieldItr->second).quantity();
(rhs[thisFieldName].set_quantity()).reset(f.nRows(),f.nCols());
}
}
}
@ -1411,8 +1423,42 @@ namespace ATC {
void ATC_Coupling::construct_transfers()
{
ATC_Method::construct_transfers();
if (!useFeMdMassMatrix_) {
// transfer for MD mass matrices based on requested intrinsic fields
if (fieldSizes_.find(TEMPERATURE) != fieldSizes_.end()) {
// classical thermodynamic heat capacity of the atoms
HeatCapacity * heatCapacity = new HeatCapacity(this);
interscaleManager_.add_per_atom_quantity(heatCapacity,
"AtomicHeatCapacity");
// atomic thermal mass matrix
nodalAtomicHeatCapacity_ = new AtfShapeFunctionRestriction(this,
heatCapacity,
shpFcn_);
interscaleManager_.add_dense_matrix(nodalAtomicHeatCapacity_,
"NodalAtomicHeatCapacity");
}
if ((fieldSizes_.find(VELOCITY) != fieldSizes_.end()) || (fieldSizes_.find(DISPLACEMENT) != fieldSizes_.end())) {
// atomic momentum mass matrix
FundamentalAtomQuantity * atomicMass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
nodalAtomicMass_ = new AtfShapeFunctionRestriction(this,
atomicMass,
shpFcn_);
interscaleManager_.add_dense_matrix(nodalAtomicMass_,
"AtomicMomentumMassMat");
}
if (fieldSizes_.find(MASS_DENSITY) != fieldSizes_.end()) {
// atomic dimensionless mass matrix
ConstantQuantity<double> * atomicOnes = new ConstantQuantity<double>(this,1);
interscaleManager_.add_per_atom_quantity(atomicOnes,"AtomicOnes");
nodalAtomicCount_ = new AtfShapeFunctionRestriction(this,
atomicOnes,
shpFcn_);
interscaleManager_.add_dense_matrix(nodalAtomicCount_,
"AtomicDimensionlessMassMat");
}
}
extrinsicModelManager_.construct_transfers();
}
@ -1573,6 +1619,57 @@ namespace ATC {
myMassMatInv(iNode,iNode) = 0.;
}
}
//---------------------------------------------------------
// compute_md_mass_matrix
// compute the mass matrix arising from only atomistic
// quadrature and contributions as a summation
//---------------------------------------------------------
void ATC_Coupling::compute_md_mass_matrix(FieldName thisField,
DIAG_MAT & massMat)
{
if (thisField == TEMPERATURE) {
massMat.shallowreset(nodalAtomicHeatCapacity_->quantity());
}
else if (thisField == DISPLACEMENT || thisField == VELOCITY) {
massMat.shallowreset(nodalAtomicMass_->quantity());
}
else if (thisField == MASS_DENSITY || thisField == SPECIES_CONCENTRATION) {
massMat.shallowreset(nodalAtomicVolume_->quantity());
}
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_Coupling::write_restart_data(string fileName, RESTART_LIST & data)
{
atomicRegulator_->pack_fields(data);
ATC_Method::write_restart_data(fileName,data);
}
//--------------------------------------------------
// read_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_Coupling::read_restart_data(string fileName, RESTART_LIST & data)
{
atomicRegulator_->pack_fields(data);
ATC_Method::read_restart_data(fileName,data);
}
//--------------------------------------------------
void ATC_Coupling::reset_nlocal()
{
ATC_Method::reset_nlocal();
atomicRegulator_->reset_nlocal();
}
//--------------------------------------------------------
void ATC_Coupling::reset_atom_materials()
{
@ -1603,6 +1700,9 @@ namespace ATC {
}
}
}
atomicRegulator_->reset_atom_materials(elementToMaterialMap_,
atomElement_);
}
//--------------------------------------------------------
@ -1678,7 +1778,7 @@ namespace ATC {
ATC_Method::post_init_integrate();
// Apply time filtering to mass matrices, if needed
if (timeFilterManager_.filter_dynamics() && !useFeMdMassMatrix_) {
if ((atomToElementMapType_ == EULERIAN) && timeFilterManager_.filter_dynamics() && !useFeMdMassMatrix_) {
map<FieldName,int>::const_iterator field;
for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
FieldName thisField = field->first;
@ -1706,6 +1806,16 @@ namespace ATC {
ATC_Method::pre_exchange();
}
//--------------------------------------------------------
// pre_force
// prior to calculation of forces
//--------------------------------------------------------
void ATC_Coupling::pre_force()
{
ATC_Method::pre_force();
atomicRegulator_->pre_force();
}
//--------------------------------------------------------
void ATC_Coupling::post_force()
{
@ -1746,6 +1856,93 @@ namespace ATC {
extrinsicModelManager_.post_force();
}
//--------------------------------------------------------
// post_final_integrate
// integration after the second stage lammps atomic
// update of Verlet step 2
//--------------------------------------------------------
void ATC_Coupling::post_final_integrate()
{
double dt = lammpsInterface_->dt();
// update of atomic contributions for fractional step methods
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->pre_final_integrate1(dt);
}
// Set sources
prescribedDataMgr_->set_sources(time()+0.5*dt,sources_);
extrinsicModelManager_.pre_final_integrate();
bool needsSources = false;
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
if ((_tiIt_->second)->has_final_predictor()) {
needsSources = true;
break;
}
}
if (needsSources) {
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(intrinsicMask_,fields_,atomicSources_);
}
atomicRegulator_->apply_pre_corrector(dt,lammpsInterface_->ntimestep());
// Compute atom-integrated rhs
// parallel communication happens within FE_Engine
compute_rhs_vector(intrinsicMask_,fields_,rhs_,FE_DOMAIN);
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->add_to_rhs();
}
atomicRegulator_->add_to_rhs(rhs_);
// Compute and add atomic contributions to FE equations
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate1(dt);
}
// fix nodes, non-group bcs applied through FE
set_fixed_nodes();
// corrector step extrinsic model
extrinsicModelManager_.post_final_integrate();
// set state-based sources
needsSources = false;
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
if ((_tiIt_->second)->has_final_corrector()) {
needsSources = true;
break;
}
}
if (needsSources) {
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(intrinsicMask_,fields_,atomicSources_);
}
// Finish update of FE velocity
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate2(dt);
}
// apply corrector phase of thermostat
atomicRegulator_->apply_post_corrector(dt,lammpsInterface_->ntimestep());
// final phase of time integration
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate3(dt);
}
// Fix nodes, non-group bcs applied through FE
set_fixed_nodes();
update_time(0.5);
output();
lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+1); // adds next step to computes
}
//=================================================================
//
//=================================================================
@ -1756,6 +1953,7 @@ namespace ATC {
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->finish();
}
atomicRegulator_->finish();
}

View File

@ -53,7 +53,7 @@ namespace ATC {
virtual void reset_atoms(){};
/** pre force */
virtual void pre_force() {};
virtual void pre_force();
/** post force */
virtual void post_force();
@ -72,14 +72,16 @@ namespace ATC {
/** Predictor phase, executed after Verlet */
virtual void post_init_integrate();
#ifdef OBSOLETE
/** Corrector phase, executed before Verlet */
virtual void pre_final_integrate(){};
#endif
/** Corrector phase, executed after Verlet*/
virtual void post_final_integrate() {lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+1);};
virtual void post_final_integrate();
#ifdef OBSOLETE
{lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+1);};
#endif
/** pre/post atomic force calculation in minimize */
virtual void min_pre_force(){};
virtual void min_post_force(){};
@ -224,6 +226,8 @@ namespace ATC {
//---------------------------------------------------------------
/*@{*/
void pack_fields(RESTART_LIST & data);
virtual void read_restart_data(std::string fileName_, RESTART_LIST & data);
virtual void write_restart_data(std::string fileName_, RESTART_LIST & data);
void output() { ATC_Method::output(); }
/*@}*/
@ -271,6 +275,8 @@ namespace ATC {
Array2D<bool> &field_mask() {return fieldMask_;};
/** create field mask */
void reset_flux_mask();
/** field mask for intrinsic integration */
Array2D<bool> intrinsicMask_;
/** wrapper for FE_Engine's compute_flux functions */
void compute_flux(const Array2D<bool> & rhs_mask,
const FIELDS &fields,
@ -307,6 +313,9 @@ namespace ATC {
void compute_mass_matrix(FieldName thisField, PhysicsModel * physicsModel = NULL);
/** updates filtering of MD contributions */
void update_mass_matrix(FieldName thisField);
/** compute the mass matrix components coming from MD integration */
virtual void compute_md_mass_matrix(FieldName thisField,
DIAG_MAT & massMats);
private: /** methods */
ATC_Coupling(); // do not define
@ -328,6 +337,8 @@ namespace ATC {
virtual void construct_molecule_transfers();
/** sets up accumulant & interpolant */
virtual void construct_interpolant();
/** reset number of local atoms */
virtual void reset_nlocal();
//---------------------------------------------------------------
/** status */
@ -383,6 +394,12 @@ namespace ATC {
MatrixDependencyManager<DenseMatrix, bool> * elementMask_;
MatrixDependencyManager<DenseMatrix, bool> * elementMaskMass_;
MatrixDependencyManager<DenseMatrix, bool> * elementMaskMassMd_;
/** operator to compute the mass matrix for the momentum equation from MD integration */
AtfShapeFunctionRestriction * nodalAtomicMass_;
/** operator to compute the dimensionless mass matrix from MD integration */
AtfShapeFunctionRestriction * nodalAtomicCount_;
/** operator to compute mass matrix from MD */
AtfShapeFunctionRestriction * nodalAtomicHeatCapacity_;
MatrixDependencyManager<DenseMatrix, bool> * create_full_element_mask();
MatrixDependencyManager<DenseMatrix, int> * create_element_set_mask(const std::string & elementSetName);
LargeToSmallAtomMap * internalToMask_;

View File

@ -30,7 +30,9 @@ namespace ATC {
string matParamFile,
ExtrinsicModelType extrinsicModel)
: ATC_Coupling(groupName,perAtomArray,thisFix),
#ifdef OBSOLETE
nodalAtomicHeatCapacity_(NULL),
#endif
nodalAtomicKineticTemperature_(NULL),
nodalAtomicConfigurationalTemperature_(NULL)
{
@ -145,25 +147,10 @@ namespace ATC {
}
// reset integration field mask
temperatureMask_.reset(NUM_FIELDS,NUM_FLUX);
temperatureMask_ = false;
intrinsicMask_.reset(NUM_FIELDS,NUM_FLUX);
intrinsicMask_ = false;
for (int i = 0; i < NUM_FLUX; i++)
temperatureMask_(TEMPERATURE,i) = fieldMask_(TEMPERATURE,i);
}
//--------------------------------------------------------
// construct_methods
// have managers instantiate requested algorithms
// and methods
//--------------------------------------------------------
void ATC_CouplingEnergy::construct_methods()
{
ATC_Coupling::construct_methods();
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_methods();
}
atomicRegulator_->construct_methods();
intrinsicMask_(TEMPERATURE,i) = fieldMask_(TEMPERATURE,i);
}
//--------------------------------------------------------
@ -270,7 +257,7 @@ namespace ATC {
TEMPERATURE);
interscaleManager_.add_dense_matrix(nodalAtomicTemperature,
"NodalAtomicTemperature");
#ifdef OBSOLETE
if (!useFeMdMassMatrix_) {
// classical thermodynamic heat capacity of the atoms
HeatCapacity * heatCapacity = new HeatCapacity(this);
@ -284,6 +271,7 @@ namespace ATC {
interscaleManager_.add_dense_matrix(nodalAtomicHeatCapacity_,
"NodalAtomicHeatCapacity");
}
#endif
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_transfers();
}
@ -326,7 +314,7 @@ namespace ATC {
}
}
}
#ifdef OBSOLETE
//---------------------------------------------------------
// compute_md_mass_matrix
// compute the mass matrix arising from only atomistic
@ -338,20 +326,8 @@ namespace ATC {
if (thisField == TEMPERATURE)
massMat.reset(nodalAtomicHeatCapacity_->quantity());
}
//--------------------------------------------------------
// finish
// final clean up after a run
//--------------------------------------------------------
void ATC_CouplingEnergy::finish()
{
// base class
ATC_Coupling::finish();
atomicRegulator_->finish();
}
#endif
//--------------------------------------------------------
// modify
// parses inputs and modifies state of the filter
@ -412,198 +388,6 @@ namespace ATC {
}
//--------------------------------------------------
// pack_fields
// bundle all allocated field matrices into a list
// for output needs
//--------------------------------------------------
void ATC_CouplingEnergy::pack_thermal_fields(RESTART_LIST & data)
{
atomicRegulator_->pack_fields(data);
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingEnergy::write_restart_data(string fileName, RESTART_LIST & data)
{
pack_thermal_fields(data);
ATC_Method::write_restart_data(fileName,data);
}
//--------------------------------------------------
// read_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingEnergy::read_restart_data(string fileName, RESTART_LIST & data)
{
pack_thermal_fields(data);
ATC_Method::read_restart_data(fileName,data);
}
//--------------------------------------------------
void ATC_CouplingEnergy::reset_nlocal()
{
ATC_Coupling::reset_nlocal();
atomicRegulator_->reset_nlocal();
}
//--------------------------------------------------
// reset_atom_materials
// update the atom materials map
//--------------------------------------------------
void ATC_CouplingEnergy::reset_atom_materials()
{
ATC_Coupling::reset_atom_materials();
atomicRegulator_->reset_atom_materials(elementToMaterialMap_,
atomElement_);
}
#ifdef OBSOLETE
//--------------------------------------------------------
// mid_init_integrate
// time integration between the velocity update and
// the position lammps update of Verlet step 1
//--------------------------------------------------------
void ATC_CouplingEnergy::mid_init_integrate()
{
// CONTINUOUS VELOCITY UPDATE
ATC_Coupling::mid_init_integrate();
double dt = lammpsInterface_->dt();
// Compute nodal velocity at n+1/2
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->mid_initial_integrate1(dt);
}
atomicRegulator_->apply_mid_predictor(dt,lammpsInterface_->ntimestep());
extrinsicModelManager_.mid_init_integrate();
}
//--------------------------------------------------------
// post_init_integrate
// time integration after the lammps atomic updates of
// Verlet step 1
//--------------------------------------------------------
void ATC_CouplingEnergy::post_init_integrate()
{
double dt = lammpsInterface_->dt();
// Compute nodal velocity at n+1
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_initial_integrate1(dt);
}
// Update kinetostat quantities if displacement is being regulated
atomicRegulator_->apply_post_predictor(dt,lammpsInterface_->ntimestep());
// Update extrisic model
extrinsicModelManager_.post_init_integrate();
// fixed values, non-group bcs handled through FE
set_fixed_nodes();
update_time(0.5);
ATC_Coupling::post_init_integrate();
}
#endif
//--------------------------------------------------------
// post_final_integrate
// integration after the second stage lammps atomic
// update of Verlet step 2
//--------------------------------------------------------
void ATC_CouplingEnergy::post_final_integrate()
{
double dt = lammpsInterface_->dt();
// update changes in atomic energy or from atomic work, if needed
// this is here to simplify computing changes in total atomic energy
// even though all the data needed is available by pre_final_integrate
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->pre_final_integrate1(dt);
}
// Set prescribed sources for current time
prescribedDataMgr_->set_sources(time()+0.5*dt,sources_);
// predictor step in extrinsic model
extrinsicModelManager_.pre_final_integrate();
// predict thermostat contributions
// compute sources based on predicted FE temperature
if (timeIntegrators_[TEMPERATURE]->has_final_predictor()) {
// set state-based sources
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(temperatureMask_,fields_,atomicSources_);
}
// Compute thermostat forces
atomicRegulator_->apply_pre_corrector(dt,lammpsInterface_->ntimestep());
// Determine FE contributions to d theta/dt
// Compute atom-integrated rhs
// parallel communication happens within FE_Engine
// Determine FE contributions to dT/dt-----------------------
compute_rhs_vector(temperatureMask_,fields_,rhs_,FE_DOMAIN);
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->add_to_rhs();
}
// For flux matching, add appropriate fraction of "drag" power
atomicRegulator_->add_to_rhs(rhs_);
// final phase predictor step
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate1(dt);
}
// fix nodes, non-group bcs applied through FE
set_fixed_nodes();
// corrector step extrinsic model
extrinsicModelManager_.post_final_integrate();
// correct thermostat and finish
if (timeIntegrators_[TEMPERATURE]->has_final_corrector()) {
// set state-based sources
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(temperatureMask_,fields_,atomicSources_);
}
// finish FE temperature update
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate2(dt);
}
// apply corrector phase of thermostat
atomicRegulator_->apply_post_corrector(dt,lammpsInterface_->ntimestep());
// finalalize time filtering
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate3(dt);
}
// Fix nodes, non-group bcs applied through FE
set_fixed_nodes();
update_time(0.5);
output();
ATC_Coupling::post_final_integrate(); // adds next step to computes
}
//--------------------------------------------------------------------
// compute_vector
//--------------------------------------------------------------------

View File

@ -51,28 +51,12 @@ namespace ATC {
return resetMethods || (_ctiIt_->second)->need_reset();
};
/** post time integration */
virtual void finish();
#ifdef OBSOLETE
/** first time, after atomic velocity but before position integration */
virtual void mid_init_integrate();
/** first time, after atomic integration */
virtual void post_init_integrate();
#endif
/** second time substep routine */
virtual void post_final_integrate();
/** compute vector for output */
virtual double compute_vector(int n);
/** output */
virtual void output();
/** set up atom to material identification */
virtual void reset_atom_materials();
protected:
//---------------------------------------------------------------
@ -80,30 +64,22 @@ namespace ATC {
//---------------------------------------------------------------
/** constructs all data which is updated with time integration, i.e. fields */
//virtual void construct_time_integration_data();
/** create methods, e.g. time integrators, filters */
virtual void construct_methods();
/** set up data which is dependency managed */
virtual void construct_transfers();
/** sets the position/velocity of the ghost atoms */
virtual void set_ghost_atoms(){};
/** adds resetting of any thermostat arrays associated with local atom count */
virtual void reset_nlocal();
#ifdef OBSOLETE
/** compute the mass matrix components coming from MD integration */
virtual void compute_md_mass_matrix(FieldName thisField,
DIAG_MAT & massMats);
/** operator to compute mass matrix from MD */
AtfShapeFunctionRestriction * nodalAtomicHeatCapacity_;
#endif
/** physics specific filter initialization */
void init_filter();
/** field mask for velocity integration */
Array2D<bool> temperatureMask_;
double compute_lambda_power(int gid);
@ -116,11 +92,6 @@ namespace ATC {
/** workspace matrices for output */
DENS_MAT _keTemp_, _peTemp_;
// Add in fields for restarting
virtual void read_restart_data(std::string fileName_, RESTART_LIST & data);
virtual void write_restart_data(std::string fileName_, RESTART_LIST & data);
void pack_thermal_fields(RESTART_LIST & data);
};
};

View File

@ -186,25 +186,11 @@ namespace ATC {
}
// reset integration field mask
speciesMask_.reset(NUM_FIELDS,NUM_FLUX);
speciesMask_ = false;
intrinsicMask_.reset(NUM_FIELDS,NUM_FLUX);
intrinsicMask_ = false;
}
//--------------------------------------------------------
// construct_methods
// have managers instantiate requested algorithms
// and methods
//--------------------------------------------------------
void ATC_CouplingMass::construct_methods()
{
ATC_Coupling::construct_methods();
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_methods();
}
atomicRegulator_->construct_methods();
}
void ATC_CouplingMass::construct_transfers()
{
ATC_Coupling::construct_transfers();
@ -259,7 +245,7 @@ namespace ATC {
ATC_Coupling::init_filter();
}
#ifdef OBSOLETE
void ATC_CouplingMass::compute_md_mass_matrix(FieldName thisField,
DIAG_MAT & massMat)
{
@ -269,37 +255,8 @@ namespace ATC {
massMat.reset(nodalAtomicVolume_->quantity());
}
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingMass::write_restart_data(string fileName, RESTART_LIST & data)
{
ATC_Method::write_restart_data(fileName,data);
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingMass::read_restart_data(string fileName, RESTART_LIST & data)
{
ATC_Method::read_restart_data(fileName,data);
}
//--------------------------------------------------------
// pre_force
// prior to calculation of forces
//--------------------------------------------------------
void ATC_CouplingMass::pre_force()
{
ATC_Coupling::pre_force();
atomicRegulator_->pre_force();
}
#endif
//WIP_JAT consolidate to coupling when we handle the temperature correctly
//--------------------------------------------------------
// pre_exchange
// prior to exchange of atoms
@ -318,123 +275,6 @@ namespace ATC {
}
}
#ifdef OBSOLETE
//--------------------------------------------------------
// mid_init_integrate
// time integration between the velocity update and
// the position lammps update of Verlet step 1
//--------------------------------------------------------
void ATC_CouplingMass::mid_init_integrate()
{
ATC_Coupling::mid_init_integrate();
double dt = lammpsInterface_->dt();
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->mid_initial_integrate1(dt);
}
atomicRegulator_->apply_mid_predictor(dt,lammpsInterface_->ntimestep());
extrinsicModelManager_.mid_init_integrate();
}
//--------------------------------------------------------
// post_init_integrate
// time integration after the lammps atomic updates of
// Verlet step 1
//--------------------------------------------------------
void ATC_CouplingMass::post_init_integrate()
{
double dt = lammpsInterface_->dt();
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_initial_integrate1(dt);
}
atomicRegulator_->apply_post_predictor(dt,lammpsInterface_->ntimestep());
extrinsicModelManager_.post_init_integrate();
set_fixed_nodes();
update_time(0.5); // half step
ATC_Coupling::post_init_integrate();
}
#endif
//--------------------------------------------------------
// post_final_integrate
// integration after the second stage lammps atomic
// update of Verlet step 2
//--------------------------------------------------------
void ATC_CouplingMass::post_final_integrate()
{
double dt = lammpsInterface_->dt();
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->pre_final_integrate1(dt);
}
prescribedDataMgr_->set_sources(time()+0.5*dt,sources_);
extrinsicModelManager_.pre_final_integrate();
if (timeIntegrators_[MASS_DENSITY]->has_final_predictor()) {
// set state-based sources
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
compute_atomic_sources(speciesMask_,fields_,atomicSources_);
}
// set state-based RHS
// Determine FE contributions to dv/dt-----------------------
// Compute atom-integrated rhs
// parallel communication happens within FE_Engine
compute_rhs_vector(speciesMask_,fields_,rhs_,FE_DOMAIN);
// Compute and add atomic contributions to FE equations
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->add_to_rhs();
}
atomicRegulator_->add_to_rhs(rhs_);
// final phase predictor step
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate1(dt);
}
set_fixed_nodes();
// corrector step extrinsic model
extrinsicModelManager_.post_final_integrate();
if (timeIntegrators_[MASS_DENSITY]->has_final_corrector()) {
// set state-based sources
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
compute_atomic_sources(speciesMask_,fields_,atomicSources_);
}
// finish FE temperature update
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate2(dt);
}
// apply corrector phase of thermostat
atomicRegulator_->apply_post_corrector(dt,lammpsInterface_->ntimestep());
// finalize time integration
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate3(dt);
}
// Fix nodes, non-group bcs applied through FE
set_fixed_nodes();
update_time(0.5);
output();
ATC_Coupling::post_final_integrate(); // addstep for computes
}
//--------------------------------------------------------
// output
// does post-processing steps and outputs data

View File

@ -49,20 +49,9 @@ namespace ATC {
/** pre time integration */
virtual void initialize();
/** prior to force */
virtual void pre_force();
/** prior to exchange */
virtual void pre_exchange();
virtual void reset_atoms() { resetNlocal_=true;}
#ifdef OBSOLETE
/** first time, after atomic velocity but before position integration */
virtual void mid_init_integrate();
/** first time, after atomic integration */
virtual void post_init_integrate();
#endif
/** second time, after atomic integration */
virtual void post_final_integrate();
/** compute vector for output */
virtual double compute_vector(int n);
@ -78,28 +67,19 @@ namespace ATC {
//---------------------------------------------------------------
/** constructs all data which is updated with time integration, i.e. fields */
//virtual void construct_time_integration_data();
/** create methods, e.g. time integrators, filters */
virtual void construct_methods();
/** set up data which is dependency managed */
virtual void construct_transfers();
/** sets the position/velocity of the ghost atoms */
virtual void set_ghost_atoms(){};
#ifdef OBSOLETE
/** compute the mass matrix components coming from MD integration */
virtual void compute_md_mass_matrix(FieldName thisField,
DIAG_MAT & massMats);
#endif
/** physics specific filter initialization */
void init_filter();
/** field mask for velocity integration */
Array2D<bool> speciesMask_;
// Add in fields for restarting
virtual void read_restart_data(std::string fileName_, RESTART_LIST & data);
virtual void write_restart_data(std::string fileName_, RESTART_LIST & data);
// DATA structures for tracking individual species and molecules
FIELD_POINTERS atomicFields_;

View File

@ -33,8 +33,10 @@ namespace ATC {
PhysicsType intrinsicModel,
ExtrinsicModelType extrinsicModel)
: ATC_Coupling(groupName,perAtomArray,thisFix),
#ifdef OBSOLETE
nodalAtomicMass_(NULL),
nodalAtomicCount_(NULL),
#endif
refPE_(0)
{
// Allocate PhysicsModel
@ -176,31 +178,16 @@ namespace ATC {
}
// reset integration field mask
velocityMask_.reset(NUM_FIELDS,NUM_FLUX);
velocityMask_ = false;
intrinsicMask_.reset(NUM_FIELDS,NUM_FLUX);
intrinsicMask_ = false;
for (int i = 0; i < NUM_FLUX; i++)
velocityMask_(VELOCITY,i) = fieldMask_(VELOCITY,i);
intrinsicMask_(VELOCITY,i) = fieldMask_(VELOCITY,i);
refPE_=0;
refPE_=potential_energy();
}
//--------------------------------------------------------
// construct_methods
// have managers instantiate requested algorithms
// and methods
//--------------------------------------------------------
void ATC_CouplingMomentum::construct_methods()
{
ATC_Coupling::construct_methods();
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_methods();
}
atomicRegulator_->construct_methods();
}
//--------------------------------------------------------
// construct_transfers
// constructs needed transfer operators
@ -265,7 +252,7 @@ namespace ATC {
interscaleManager_.add_dense_matrix(nodalAtomicDisplacement,
"NodalAtomicDisplacement");
}
#ifdef OBSOLETE
// atomic mass matrix data
if (!useFeMdMassMatrix_) {
// atomic momentum mass matrix
@ -285,7 +272,7 @@ namespace ATC {
interscaleManager_.add_dense_matrix(nodalAtomicCount_,
"AtomicDimensionlessMassMat");
}
#endif
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_transfers();
}
@ -306,7 +293,7 @@ namespace ATC {
// nothing needed in other cases since kinetostat force is balanced by boundary flux in FE equations
atomicRegulator_->reset_lambda_contribution(nodalAtomicFieldsRoc_[VELOCITY].quantity());
}
#ifdef OBSOLETE
//---------------------------------------------------------
// compute_md_mass_matrix
// compute the mass matrix arising from only atomistic
@ -323,19 +310,7 @@ namespace ATC {
massMat.reset(nodalAtomicCount_->quantity());
}
}
//--------------------------------------------------------
// finish
// final clean up after a run
//--------------------------------------------------------
void ATC_CouplingMomentum::finish()
{
// base class
ATC_Coupling::finish();
atomicRegulator_->finish();
}
#endif
//--------------------------------------------------------
// modify
// parses inputs and modifies state of the filter
@ -402,209 +377,6 @@ namespace ATC {
}
//--------------------------------------------------
// pack_fields
// bundle all allocated field matrices into a list
// for output needs
//--------------------------------------------------
void ATC_CouplingMomentum::pack_elastic_fields(RESTART_LIST & data)
{
atomicRegulator_->pack_fields(data);
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingMomentum::write_restart_data(string fileName, RESTART_LIST & data)
{
pack_elastic_fields(data);
ATC_Method::write_restart_data(fileName,data);
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingMomentum::read_restart_data(string fileName, RESTART_LIST & data)
{
pack_elastic_fields(data);
ATC_Method::read_restart_data(fileName,data);
}
//--------------------------------------------------------
void ATC_CouplingMomentum::reset_nlocal()
{
ATC_Coupling::reset_nlocal();
atomicRegulator_->reset_nlocal();
}
//--------------------------------------------------
// reset_atom_materials
// update the atom materials map
//--------------------------------------------------
void ATC_CouplingMomentum::reset_atom_materials()
{
ATC_Coupling::reset_atom_materials();
atomicRegulator_->reset_atom_materials(elementToMaterialMap_,
atomElement_);
}
#ifdef OBSOLETE
//--------------------------------------------------------
// mid_init_integrate
// time integration between the velocity update and
// the position lammps update of Verlet step 1
//--------------------------------------------------------
void ATC_CouplingMomentum::mid_init_integrate()
{
// CONTINUOUS VELOCITY UPDATE
ATC_Coupling::mid_init_integrate();
double dt = lammpsInterface_->dt();
// Compute nodal velocity at n+1/2
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->mid_initial_integrate1(dt);
}
atomicRegulator_->apply_mid_predictor(dt,lammpsInterface_->ntimestep());
extrinsicModelManager_.mid_init_integrate();
}
//--------------------------------------------------------
// post_init_integrate
// time integration after the lammps atomic updates of
// Verlet step 1
//--------------------------------------------------------
void ATC_CouplingMomentum::post_init_integrate()
{
// CONTINUOUS DISPLACEMENT UPDATE
double dt = lammpsInterface_->dt();
// Compute nodal velocity at n+1
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_initial_integrate1(dt);
}
// Update kinetostat quantities if displacement is being regulated
atomicRegulator_->apply_post_predictor(dt,lammpsInterface_->ntimestep());
// Update extrisic model
extrinsicModelManager_.post_init_integrate();
// fixed values, non-group bcs handled through FE
set_fixed_nodes();
// update time by a half dt
update_time(0.5);
ATC_Coupling::post_init_integrate();
}
#endif
//--------------------------------------------------------
// pre_final_integrate
// integration before the second stage lammps atomic
// update of Verlet step 2
//--------------------------------------------------------
void ATC_CouplingMomentum::pre_final_integrate()
{
ATC_Coupling::pre_final_integrate();
}
//--------------------------------------------------------
// post_final_integrate
// integration after the second stage lammps atomic
// update of Verlet step 2
//--------------------------------------------------------
void ATC_CouplingMomentum::post_final_integrate()
{
// COMPUTE FORCES FOR FE VELOCITY RHS
double dt = lammpsInterface_->dt();
// updating of data based on atomic forces
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->pre_final_integrate1(dt);
}
// Set prescribed sources for current time
prescribedDataMgr_->set_sources(time()+0.5*dt,sources_);
// predictor step in extrinsic model
extrinsicModelManager_.pre_final_integrate();
if (timeIntegrators_[VELOCITY]->has_final_predictor()) {
// set state-based sources
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(velocityMask_,fields_,atomicSources_);
}
// Compute kinetostat forces and add kinetostat contributions to FE equations
atomicRegulator_->apply_pre_corrector(dt,lammpsInterface_->ntimestep()); // computes but does not apply kstat, and only for StressFlux
// set state-based RHS
// Determine FE contributions to dv/dt-----------------------
// Compute atom-integrated rhs
// parallel communication happens within FE_Engine
compute_rhs_vector(velocityMask_,fields_,rhs_,FE_DOMAIN);
// Compute and add atomic contributions to FE equations
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->add_to_rhs();
}
// add in kinetostat contributions to FE equations
atomicRegulator_->add_to_rhs(rhs_);
// final phase predictor step
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate1(dt);
}
// fix nodes, non-group bcs applied through FE
set_fixed_nodes();
// CONTINUOUS VELOCITY RHS UPDATE
// corrector step extrinsic model
extrinsicModelManager_.post_final_integrate();
if (timeIntegrators_[VELOCITY]->has_final_corrector()) {
// set state-based sources
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(velocityMask_,fields_,atomicSources_);
}
// Finish update of FE velocity
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate2(dt);
}
// Apply kinetostat to atoms
atomicRegulator_->apply_post_corrector(dt,lammpsInterface_->ntimestep());
// finalize time integration
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate3(dt);
}
// Fix nodes, non-group bcs applied through FE
set_fixed_nodes();
// update time by a half dt
update_time(0.5);
output();
ATC_Coupling::post_final_integrate(); // addstep for computes
}
//--------------------------------------------------------
// min_pre_force
// add to interatomic forces for minimize

View File

@ -52,21 +52,6 @@ namespace ATC {
return resetMethods || (_ctiIt_->second)->need_reset();
};
/** post time integration */
virtual void finish();
#ifdef OBSOLETE
/** first time, after atomic velocity but before position integration */
virtual void mid_init_integrate();
/** first time, after atomic integration */
virtual void post_init_integrate();
#endif
/** second time, before atomic integration */
virtual void pre_final_integrate();
/** second time, after atomic integration */
virtual void post_final_integrate();
/** pre/post atomic force calculation in minimize */
virtual void min_pre_force();
virtual void min_post_force();
@ -82,10 +67,6 @@ namespace ATC {
/** output routines */
virtual void output(void);
/** set up atom to material identification */
virtual void reset_atom_materials();
protected:
//---------------------------------------------------------------
@ -93,35 +74,22 @@ namespace ATC {
//---------------------------------------------------------------
/** constructs all data which is updated with time integration, i.e. fields */
//virtual void construct_time_integration_data();
/** create methods, e.g. time integrators, filters */
virtual void construct_methods();
/** set up data which is dependency managed */
virtual void construct_transfers();
/** adds resetting of any kinetostat arrays associated with local atom count */
virtual void reset_nlocal();
#ifdef OBSOLETE
/** compute the mass matrix components coming from MD integration */
virtual void compute_md_mass_matrix(FieldName thisField,
DIAG_MAT & massMats);
//
/** operator to compute the mass matrix for the momentum equation from MD integration */
AtfShapeFunctionRestriction * nodalAtomicMass_;
/** operator to compute the dimensionless mass matrix from MD integration */
AtfShapeFunctionRestriction * nodalAtomicCount_;
#endif
/** physics specific filter initialization */
void init_filter();
/** field mask for velocity integration */
Array2D<bool> velocityMask_;
// Add in fields for restarting
virtual void read_restart_data(std::string fileName_, RESTART_LIST & data);
virtual void write_restart_data(std::string fileName_, RESTART_LIST & data);
void pack_elastic_fields(RESTART_LIST & data);
// data
double refPE_;

View File

@ -31,9 +31,11 @@ namespace ATC {
string matParamFile,
ExtrinsicModelType extrinsicModel)
: ATC_Coupling(groupName,perAtomArray,thisFix),
#ifdef OBSOLETE
nodalAtomicMass_(NULL),
nodalAtomicCount_(NULL),
nodalAtomicHeatCapacity_(NULL),
#endif
nodalAtomicKineticTemperature_(NULL),
nodalAtomicConfigurationalTemperature_(NULL),
refPE_(0)
@ -116,6 +118,7 @@ namespace ATC {
// which is why the time integrator is initialized first
// other initializations
#ifdef OBSOLETE
if (reset_methods()) {
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
@ -123,6 +126,7 @@ namespace ATC {
}
atomicRegulator_->force_reset();
}
#endif
if (reset_methods()) {
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->initialize();
@ -191,34 +195,17 @@ namespace ATC {
}
// reset integration field mask
velocityMask_.reset(NUM_FIELDS,NUM_FLUX);
velocityMask_ = false;
intrinsicMask_.reset(NUM_FIELDS,NUM_FLUX);
intrinsicMask_ = false;
for (int i = 0; i < NUM_FLUX; i++)
velocityMask_(VELOCITY,i) = fieldMask_(VELOCITY,i);
temperatureMask_.reset(NUM_FIELDS,NUM_FLUX);
temperatureMask_ = false;
intrinsicMask_(VELOCITY,i) = fieldMask_(VELOCITY,i);
for (int i = 0; i < NUM_FLUX; i++)
temperatureMask_(TEMPERATURE,i) = fieldMask_(TEMPERATURE,i);
intrinsicMask_(TEMPERATURE,i) = fieldMask_(TEMPERATURE,i);
refPE_=0;
refPE_=potential_energy();
}
//--------------------------------------------------------
// construct_methods
// have managers instantiate requested algorithms
// and methods
//--------------------------------------------------------
void ATC_CouplingMomentumEnergy::construct_methods()
{
ATC_Coupling::construct_methods();
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_methods();
}
atomicRegulator_->construct_methods();
}
//--------------------------------------------------------
// construct_transfers
// constructs needed transfer operators
@ -381,7 +368,7 @@ namespace ATC {
TEMPERATURE);
interscaleManager_.add_dense_matrix(nodalAtomicTemperature,
"NodalAtomicTemperature");
#ifdef OBSOLETE
if (!useFeMdMassMatrix_) {
// atomic momentum mass matrix
FundamentalAtomQuantity * atomicMass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
@ -412,7 +399,7 @@ namespace ATC {
interscaleManager_.add_dense_matrix(nodalAtomicHeatCapacity_,
"NodalAtomicHeatCapacity");
}
#endif
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->construct_transfers();
}
@ -444,7 +431,7 @@ namespace ATC {
atomicRegulator_->reset_lambda_contribution(powerMat,TEMPERATURE);
}
}
#ifdef OBSOLETE
//---------------------------------------------------------
// compute_md_mass_matrix
// compute the mass matrix arising from only atomistic
@ -465,19 +452,7 @@ namespace ATC {
massMat.reset(nodalAtomicHeatCapacity_->quantity());
}
}
//--------------------------------------------------------
// finish
// final clean up after a run
//--------------------------------------------------------
void ATC_CouplingMomentumEnergy::finish()
{
// base class
ATC_Coupling::finish();
atomicRegulator_->finish();
}
#endif
//--------------------------------------------------------
// modify
// parses inputs and modifies state of the filter
@ -489,188 +464,6 @@ namespace ATC {
return foundMatch;
}
//--------------------------------------------------
// pack_fields
// bundle all allocated field matrices into a list
// for output needs
//--------------------------------------------------
void ATC_CouplingMomentumEnergy::pack_quantity_fields(RESTART_LIST & data)
{
atomicRegulator_->pack_fields(data);
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingMomentumEnergy::write_restart_data(string fileName, RESTART_LIST & data)
{
pack_quantity_fields(data);
ATC_Method::write_restart_data(fileName,data);
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void ATC_CouplingMomentumEnergy::read_restart_data(string fileName, RESTART_LIST & data)
{
pack_quantity_fields(data);
ATC_Method::read_restart_data(fileName,data);
}
//--------------------------------------------------
void ATC_CouplingMomentumEnergy::reset_nlocal()
{
ATC_Coupling::reset_nlocal();
atomicRegulator_->reset_nlocal();
}
//--------------------------------------------------
// reset_atom_materials
// update the atom materials map
//--------------------------------------------------
void ATC_CouplingMomentumEnergy::reset_atom_materials()
{
ATC_Coupling::reset_atom_materials();
atomicRegulator_->reset_atom_materials(elementToMaterialMap_,
atomElement_);
}
#ifdef OBSOLETE
//--------------------------------------------------------
// mid_init_integrate
// time integration between the velocity update and
// the position lammps update of Verlet step 1
//--------------------------------------------------------
void ATC_CouplingMomentumEnergy::mid_init_integrate()
{
// CONTINUOUS VELOCITY UPDATE
ATC_Coupling::mid_init_integrate();
double dt = lammpsInterface_->dt();
// Compute nodal velocity at n+1/2
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->mid_initial_integrate1(dt);
}
atomicRegulator_->apply_mid_predictor(dt,lammpsInterface_->ntimestep());
extrinsicModelManager_.mid_init_integrate();
}
//--------------------------------------------------------
// post_init_integrate
// time integration after the lammps atomic updates of
// Verlet step 1
//--------------------------------------------------------
void ATC_CouplingMomentumEnergy::post_init_integrate()
{
// CONTINUOUS DISPLACEMENT UPDATE
double dt = lammpsInterface_->dt();
// Compute nodal velocity at n+1
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_initial_integrate1(dt);
}
// Update kinetostat quantities if displacement is being regulated
atomicRegulator_->apply_post_predictor(dt,lammpsInterface_->ntimestep());
// Update extrisic model
extrinsicModelManager_.post_init_integrate();
// fixed values, non-group bcs handled through FE
set_fixed_nodes();
// update time by a half dt
update_time(0.5);
ATC_Coupling::post_init_integrate();
}
#endif
//--------------------------------------------------------
// pre_final_integrate
// integration before the second stage lammps atomic
// update of Verlet step 2
//--------------------------------------------------------
void ATC_CouplingMomentumEnergy::pre_final_integrate()
{
ATC_Coupling::pre_final_integrate();
}
//--------------------------------------------------
void ATC_CouplingMomentumEnergy::post_final_integrate()
{
// CONTINUOUS VELOCITY RHS UPDATE
double dt = lammpsInterface_->dt();
// update of atomic contributions for fractional step methods
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->pre_final_integrate1(dt);
}
// Set sources
prescribedDataMgr_->set_sources(time()+0.5*dt,sources_);
extrinsicModelManager_.pre_final_integrate();
if (timeIntegrators_[TEMPERATURE]->has_final_predictor() || timeIntegrators_[MOMENTUM]->has_final_predictor()) {
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(velocityMask_,fields_,atomicSources_);
}
atomicRegulator_->apply_pre_corrector(dt,lammpsInterface_->ntimestep());
// Compute atom-integrated rhs
// parallel communication happens within FE_Engine
compute_rhs_vector(velocityMask_,fields_,rhs_,FE_DOMAIN);
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->add_to_rhs();
}
atomicRegulator_->add_to_rhs(rhs_);
// Compute and add atomic contributions to FE equations
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate1(dt);
}
// fix nodes, non-group bcs applied through FE
set_fixed_nodes();
// corrector step extrinsic model
extrinsicModelManager_.post_final_integrate();
if (timeIntegrators_[TEMPERATURE]->has_final_corrector() || timeIntegrators_[MOMENTUM]->has_final_corrector()) {
// set state-based sources
extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
atomicRegulator_->compute_boundary_flux(fields_);
compute_atomic_sources(velocityMask_,fields_,atomicSources_);
}
// Finish update of FE velocity
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate2(dt);
}
// apply corrector phase of thermostat
atomicRegulator_->apply_post_corrector(dt,lammpsInterface_->ntimestep());
// final phase of time integration
for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
(_tiIt_->second)->post_final_integrate3(dt);
}
// Fix nodes, non-group bcs applied through FE
set_fixed_nodes();
update_time(0.5);
output();
ATC_Coupling::post_final_integrate(); // adds next step to computes
}
//--------------------------------------------------------------------
// compute_scalar : added energy
//--------------------------------------------------------------------

View File

@ -52,21 +52,6 @@ namespace ATC {
return resetMethods;
};
/** post time integration */
virtual void finish();
#ifdef OBSOLETE
/** first time, after atomic velocity but before position integration */
virtual void mid_init_integrate();
/** first time, after atomic integration */
virtual void post_init_integrate();
#endif
/** second time, before atomic integration */
virtual void pre_final_integrate();
/** second time, after atomic integration */
virtual void post_final_integrate();
/** compute scalar for output - added energy */
virtual double compute_scalar(void);
@ -78,10 +63,6 @@ namespace ATC {
/** output */
virtual void output();
/** set up atom to material identification */
virtual void reset_atom_materials();
protected:
//---------------------------------------------------------------
@ -89,18 +70,13 @@ namespace ATC {
//---------------------------------------------------------------
/** constructs all data which is updated with time integration, i.e. fields */
//virtual void construct_time_integration_data();
/** create methods, e.g. time integrators, filters */
virtual void construct_methods();
/** set up data which is dependency managed */
virtual void construct_transfers();
/** adds resetting of any kinetostat/thermostat arrays associated with local atom count */
virtual void reset_nlocal();
#ifdef OBSOLETE
/** compute the mass matrix components coming from MD integration */
virtual void compute_md_mass_matrix(FieldName thisField,
DIAG_MAT & massMats);
//
/** operator to compute the mass matrix for the momentum equation from MD integration */
AtfShapeFunctionRestriction * nodalAtomicMass_;
@ -109,16 +85,9 @@ namespace ATC {
/** operator to compute mass matrix from MD */
AtfShapeFunctionRestriction * nodalAtomicHeatCapacity_;
#endif
/** physics specific filter initialization */
void init_filter();
/** field mask for velocity integration */
Array2D<bool> velocityMask_;
/** field mask for temperature integration */
Array2D<bool> temperatureMask_;
/** kinetic temperature for post-processing */
AtfShapeFunctionMdProjection * nodalAtomicKineticTemperature_;
@ -128,11 +97,6 @@ namespace ATC {
/** workspace matrices for output */
DENS_MAT _keTemp_, _peTemp_;
// Add in fields for restarting
virtual void read_restart_data(std::string fileName_, RESTART_LIST & data);
virtual void write_restart_data(std::string fileName_, RESTART_LIST & data);
void pack_quantity_fields(RESTART_LIST & data);
// data
double refPE_;

View File

@ -17,7 +17,6 @@
#define HACK(l,m)
namespace ATC {
/**
* @class ATC_Error

View File

@ -9,9 +9,7 @@
namespace ATC
{
// forward declares
class LAMMPS_NS::PairEAM;
class LAMMPS_NS::PairEAM; // necessary for non-lib build of ATC (vs USER-ATC)
/**
* @class CbEam

View File

@ -321,7 +321,7 @@ template<typename T>
void DiagonalMatrix<T>::shallowreset(const DenseMatrix<T> &c)
{
_delete();
_data = new CloneVector<T>(*(c._data));
_data = new CloneVector<T>(c,CLONE_COL);
}
//-----------------------------------------------------------------------------
// reference indexing operator - must throw an error if i!=j

View File

@ -389,6 +389,11 @@ namespace ATC {
atomShapeFunctions,
GHOST);
interscaleManager.add_per_atom_quantity(atomFeVelocity_,field_to_prolongation_name(VELOCITY)+"Ghost");
// calculate nominal bond stiffness
int i = 0, j = 1;// HACk should be an atom and its neighbor in the boundary
double rsq = 0.0;
//k0_ = LammpsInterface_->bond_stiffness(i,j,rsq);
k0_ = LammpsInterface::instance()->bond_stiffness(i,j,rsq);
}
//--------------------------------------------------------

View File

@ -246,7 +246,7 @@ namespace ATC {
PerAtomQuantity<double> * atomForces_;
/** spring constant */
double kappa_;
double kappa_, k0_;
/** damping constant */
double gamma_;

View File

@ -786,7 +786,7 @@ double LammpsInterface::pair_force(
}
}
double LammpsInterface::pair_force(
std::pair< std::pair< int,int >,int > apair, double rsq,
std::pair< std::pair< int,int >,int > apair, double rsq,
double & fmag_over_rmag, int nbonds) const
{
int n = apair.second;
@ -800,6 +800,17 @@ double LammpsInterface::pair_force(
return pair_force(i,j, rsq,fmag_over_rmag);
}
}
double LammpsInterface::bond_stiffness(int i, int j, double rsq0) const
{
const double perturbation = 1.e-8;
double rsq1 = sqrt(rsq0)+perturbation;
rsq1 *= rsq1;
double f0,f1;
pair_force(i,j,rsq0,f0);
pair_force(i,j,rsq1,f1);
double k = (f1-f0)/perturbation;
return k;
}
double LammpsInterface::pair_cutoff() const
{

View File

@ -468,6 +468,7 @@ class LammpsInterface {
void pair_reinit() const;
int single_enable() const;
LAMMPS_NS::PairEAM * pair_eam(void) const;
double bond_stiffness(int i, int j, double rsq) const;
/*@}*/
/** \name Methods for addition/deletion of atoms*/