From dddbef699d3e31a253baf0546c9cf45f1f644234 Mon Sep 17 00:00:00 2001 From: rjones Date: Fri, 23 Aug 2013 00:14:36 +0000 Subject: [PATCH] ATC version 2.0, date: Aug22 git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10658 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/atc/ATC_Coupling.cpp | 210 ++++++++++++++++++++- lib/atc/ATC_Coupling.h | 27 ++- lib/atc/ATC_CouplingEnergy.cpp | 234 +---------------------- lib/atc/ATC_CouplingEnergy.h | 33 +--- lib/atc/ATC_CouplingMass.cpp | 170 +---------------- lib/atc/ATC_CouplingMass.h | 24 +-- lib/atc/ATC_CouplingMomentum.cpp | 246 +------------------------ lib/atc/ATC_CouplingMomentum.h | 38 +--- lib/atc/ATC_CouplingMomentumEnergy.cpp | 231 ++--------------------- lib/atc/ATC_CouplingMomentumEnergy.h | 42 +---- lib/atc/ATC_Error.h | 1 - lib/atc/CbEam.h | 4 +- lib/atc/DiagonalMatrix.h | 2 +- lib/atc/GhostManager.cpp | 5 + lib/atc/GhostManager.h | 2 +- lib/atc/LammpsInterface.cpp | 13 +- lib/atc/LammpsInterface.h | 1 + 17 files changed, 292 insertions(+), 991 deletions(-) diff --git a/lib/atc/ATC_Coupling.cpp b/lib/atc/ATC_Coupling.cpp index e9c98612fd..c6e239cc40 100644 --- a/lib/atc/ATC_Coupling.cpp +++ b/lib/atc/ATC_Coupling.cpp @@ -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 * atomicOnes = new ConstantQuantity(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::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(); } diff --git a/lib/atc/ATC_Coupling.h b/lib/atc/ATC_Coupling.h index 464db33074..7b9ff14983 100644 --- a/lib/atc/ATC_Coupling.h +++ b/lib/atc/ATC_Coupling.h @@ -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 &field_mask() {return fieldMask_;}; /** create field mask */ void reset_flux_mask(); + /** field mask for intrinsic integration */ + Array2D intrinsicMask_; /** wrapper for FE_Engine's compute_flux functions */ void compute_flux(const Array2D & 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 * elementMask_; MatrixDependencyManager * elementMaskMass_; MatrixDependencyManager * 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 * create_full_element_mask(); MatrixDependencyManager * create_element_set_mask(const std::string & elementSetName); LargeToSmallAtomMap * internalToMask_; diff --git a/lib/atc/ATC_CouplingEnergy.cpp b/lib/atc/ATC_CouplingEnergy.cpp index 066bf0b43e..4592614782 100644 --- a/lib/atc/ATC_CouplingEnergy.cpp +++ b/lib/atc/ATC_CouplingEnergy.cpp @@ -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 //-------------------------------------------------------------------- diff --git a/lib/atc/ATC_CouplingEnergy.h b/lib/atc/ATC_CouplingEnergy.h index ba26961187..7ce6063dce 100644 --- a/lib/atc/ATC_CouplingEnergy.h +++ b/lib/atc/ATC_CouplingEnergy.h @@ -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 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); - }; }; diff --git a/lib/atc/ATC_CouplingMass.cpp b/lib/atc/ATC_CouplingMass.cpp index 54fb065413..7bc7e7e3c8 100644 --- a/lib/atc/ATC_CouplingMass.cpp +++ b/lib/atc/ATC_CouplingMass.cpp @@ -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 diff --git a/lib/atc/ATC_CouplingMass.h b/lib/atc/ATC_CouplingMass.h index 38535b1479..582845a8a2 100644 --- a/lib/atc/ATC_CouplingMass.h +++ b/lib/atc/ATC_CouplingMass.h @@ -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 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_; diff --git a/lib/atc/ATC_CouplingMomentum.cpp b/lib/atc/ATC_CouplingMomentum.cpp index f81c6ec62e..6eb8529a75 100644 --- a/lib/atc/ATC_CouplingMomentum.cpp +++ b/lib/atc/ATC_CouplingMomentum.cpp @@ -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 diff --git a/lib/atc/ATC_CouplingMomentum.h b/lib/atc/ATC_CouplingMomentum.h index 8f83f527f1..d6ade45335 100644 --- a/lib/atc/ATC_CouplingMomentum.h +++ b/lib/atc/ATC_CouplingMomentum.h @@ -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 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_; diff --git a/lib/atc/ATC_CouplingMomentumEnergy.cpp b/lib/atc/ATC_CouplingMomentumEnergy.cpp index c82ecffb45..6ffaed5213 100644 --- a/lib/atc/ATC_CouplingMomentumEnergy.cpp +++ b/lib/atc/ATC_CouplingMomentumEnergy.cpp @@ -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 //-------------------------------------------------------------------- diff --git a/lib/atc/ATC_CouplingMomentumEnergy.h b/lib/atc/ATC_CouplingMomentumEnergy.h index 60410e75df..d0c78bdaee 100644 --- a/lib/atc/ATC_CouplingMomentumEnergy.h +++ b/lib/atc/ATC_CouplingMomentumEnergy.h @@ -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 velocityMask_; - - /** field mask for temperature integration */ - Array2D 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_; diff --git a/lib/atc/ATC_Error.h b/lib/atc/ATC_Error.h index 91d3b4dcff..fc356dbec8 100644 --- a/lib/atc/ATC_Error.h +++ b/lib/atc/ATC_Error.h @@ -17,7 +17,6 @@ #define HACK(l,m) - namespace ATC { /** * @class ATC_Error diff --git a/lib/atc/CbEam.h b/lib/atc/CbEam.h index 542aec9fa0..1ff7178c00 100644 --- a/lib/atc/CbEam.h +++ b/lib/atc/CbEam.h @@ -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 diff --git a/lib/atc/DiagonalMatrix.h b/lib/atc/DiagonalMatrix.h index f490816928..2ec096cb50 100644 --- a/lib/atc/DiagonalMatrix.h +++ b/lib/atc/DiagonalMatrix.h @@ -321,7 +321,7 @@ template void DiagonalMatrix::shallowreset(const DenseMatrix &c) { _delete(); - _data = new CloneVector(*(c._data)); + _data = new CloneVector(c,CLONE_COL); } //----------------------------------------------------------------------------- // reference indexing operator - must throw an error if i!=j diff --git a/lib/atc/GhostManager.cpp b/lib/atc/GhostManager.cpp index f722a2fa27..1f7fc8c838 100644 --- a/lib/atc/GhostManager.cpp +++ b/lib/atc/GhostManager.cpp @@ -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); } //-------------------------------------------------------- diff --git a/lib/atc/GhostManager.h b/lib/atc/GhostManager.h index dcacdbd3a0..c146644db2 100644 --- a/lib/atc/GhostManager.h +++ b/lib/atc/GhostManager.h @@ -246,7 +246,7 @@ namespace ATC { PerAtomQuantity * atomForces_; /** spring constant */ - double kappa_; + double kappa_, k0_; /** damping constant */ double gamma_; diff --git a/lib/atc/LammpsInterface.cpp b/lib/atc/LammpsInterface.cpp index b9b1aa8d31..565a3c9872 100644 --- a/lib/atc/LammpsInterface.cpp +++ b/lib/atc/LammpsInterface.cpp @@ -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 { diff --git a/lib/atc/LammpsInterface.h b/lib/atc/LammpsInterface.h index bb7c0f4825..e306a37aa1 100644 --- a/lib/atc/LammpsInterface.h +++ b/lib/atc/LammpsInterface.h @@ -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*/