1 // ATC headers 2 #include "ATC_CouplingMomentumEnergy.h" 3 #include "KinetoThermostat.h" 4 #include "ATC_Error.h" 5 #include "PrescribedDataManager.h" 6 7 // Other Headers 8 #include <vector> 9 #include <map> 10 #include <set> 11 #include <utility> 12 #include <typeinfo> 13 #include <iostream> 14 15 using std::string; 16 17 namespace ATC { 18 19 //-------------------------------------------------------- 20 //-------------------------------------------------------- 21 // Class ATC_CouplingMomentumEnergy 22 //-------------------------------------------------------- 23 //-------------------------------------------------------- 24 25 //-------------------------------------------------------- 26 // Constructor 27 //-------------------------------------------------------- ATC_CouplingMomentumEnergy(string groupName,double ** & perAtomArray,LAMMPS_NS::Fix * thisFix,string matParamFile,ExtrinsicModelType extrinsicModel)28 ATC_CouplingMomentumEnergy::ATC_CouplingMomentumEnergy(string groupName, 29 double ** & perAtomArray, 30 LAMMPS_NS::Fix * thisFix, 31 string matParamFile, 32 ExtrinsicModelType extrinsicModel) 33 : ATC_Coupling(groupName,perAtomArray,thisFix), 34 nodalAtomicKineticTemperature_(NULL), 35 nodalAtomicConfigurationalTemperature_(NULL), 36 refPE_(0) 37 { 38 // Allocate PhysicsModel 39 create_physics_model(THERMO_ELASTIC, matParamFile); 40 41 // create extrinsic physics model 42 if (extrinsicModel != NO_MODEL) { 43 extrinsicModelManager_.create_model(extrinsicModel,matParamFile); 44 } 45 46 // set up field data based on physicsModel 47 physicsModel_->num_fields(fieldSizes_,fieldMask_); 48 49 // Defaults 50 set_time(); 51 bndyIntType_ = FE_INTERPOLATION; 52 trackCharge_ = false; 53 54 // set up atomic regulator 55 atomicRegulator_ = new KinetoThermostat(this); 56 57 // set up physics specific time integrator and thermostat 58 trackDisplacement_ = true; 59 fieldSizes_[DISPLACEMENT] = fieldSizes_[VELOCITY]; 60 timeIntegrators_[VELOCITY] = new MomentumTimeIntegrator(this,TimeIntegrator::FRACTIONAL_STEP); 61 timeIntegrators_[TEMPERATURE] = new ThermalTimeIntegrator(this,TimeIntegrator::FRACTIONAL_STEP); 62 ghostManager_.set_boundary_dynamics(GhostManager::PRESCRIBED); 63 64 // default physics 65 temperatureDef_ = KINETIC; 66 67 // output variable vector info: 68 // output[1] = total coarse scale mechanical kinetic energy 69 // output[2] = total coarse scale mechanical potential energy 70 // output[3] = total coarse scale mechanical energy 71 // output[1] = total coarse scale thermal energy 72 // output[2] = average temperature 73 scalarFlag_ = 1; 74 vectorFlag_ = 1; 75 sizeVector_ = 5; 76 scalarVectorFreq_ = 1; 77 extVector_ = 1; 78 if (extrinsicModel != NO_MODEL) 79 sizeVector_ += extrinsicModelManager_.size_vector(sizeVector_); 80 } 81 82 //-------------------------------------------------------- 83 // Destructor 84 //-------------------------------------------------------- ~ATC_CouplingMomentumEnergy()85 ATC_CouplingMomentumEnergy::~ATC_CouplingMomentumEnergy() 86 { 87 // clear out all managed memory to avoid conflicts with dependencies on class member data 88 interscaleManager_.clear(); 89 } 90 91 //-------------------------------------------------------- 92 // initialize 93 // sets up all the necessary data 94 //-------------------------------------------------------- initialize()95 void ATC_CouplingMomentumEnergy::initialize() 96 { 97 // clear displacement entries if requested 98 if (!trackDisplacement_) { 99 fieldSizes_.erase(DISPLACEMENT); 100 for (int i = 0; i < NUM_FLUX; i++) 101 fieldMask_(DISPLACEMENT,i) = false; 102 } 103 104 // Base class initalizations 105 ATC_Coupling::initialize(); 106 107 // reset integration field mask 108 intrinsicMask_.reset(NUM_FIELDS,NUM_FLUX); 109 intrinsicMask_ = false; 110 for (int i = 0; i < NUM_FLUX; i++) 111 intrinsicMask_(VELOCITY,i) = fieldMask_(VELOCITY,i); 112 for (int i = 0; i < NUM_FLUX; i++) 113 intrinsicMask_(TEMPERATURE,i) = fieldMask_(TEMPERATURE,i); 114 115 refPE_=0; 116 refPE_=potential_energy(); 117 } 118 119 //-------------------------------------------------------- 120 // construct_transfers 121 // constructs needed transfer operators 122 //-------------------------------------------------------- construct_transfers()123 void ATC_CouplingMomentumEnergy::construct_transfers() 124 { 125 ATC_Coupling::construct_transfers(); 126 127 // momentum of each atom 128 AtomicMomentum * atomicMomentum = new AtomicMomentum(this); 129 interscaleManager_.add_per_atom_quantity(atomicMomentum, 130 "AtomicMomentum"); 131 132 // nodal momentum for RHS 133 AtfShapeFunctionRestriction * nodalAtomicMomentum = new AtfShapeFunctionRestriction(this, 134 atomicMomentum, 135 shpFcn_); 136 interscaleManager_.add_dense_matrix(nodalAtomicMomentum, 137 "NodalAtomicMomentum"); 138 139 // nodal forces 140 FundamentalAtomQuantity * atomicForce = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE); 141 AtfShapeFunctionRestriction * nodalAtomicForce = new AtfShapeFunctionRestriction(this, 142 atomicForce, 143 shpFcn_); 144 interscaleManager_.add_dense_matrix(nodalAtomicForce, 145 "NodalAtomicForce"); 146 147 // nodal velocity derived only from atoms 148 AtfShapeFunctionMdProjection * nodalAtomicVelocity = new AtfShapeFunctionMdProjection(this, 149 nodalAtomicMomentum, 150 VELOCITY); 151 interscaleManager_.add_dense_matrix(nodalAtomicVelocity, 152 "NodalAtomicVelocity"); 153 154 if (trackDisplacement_) { 155 // mass-weighted (center-of-mass) displacement of each atom 156 AtomicMassWeightedDisplacement * atomicMassWeightedDisplacement; 157 if (needXrefProcessorGhosts_ || groupbitGhost_) { // explicit construction on internal group 158 PerAtomQuantity<double> * atomReferencePositions = interscaleManager_.per_atom_quantity("AtomicInternalReferencePositions"); 159 atomicMassWeightedDisplacement = new AtomicMassWeightedDisplacement(this,atomPositions_, 160 atomMasses_, 161 atomReferencePositions, 162 INTERNAL); 163 } 164 else 165 atomicMassWeightedDisplacement = new AtomicMassWeightedDisplacement(this); 166 interscaleManager_.add_per_atom_quantity(atomicMassWeightedDisplacement, 167 "AtomicMassWeightedDisplacement"); 168 169 // nodal (RHS) mass-weighted displacement 170 AtfShapeFunctionRestriction * nodalAtomicMassWeightedDisplacement = new AtfShapeFunctionRestriction(this, 171 atomicMassWeightedDisplacement, 172 shpFcn_); 173 interscaleManager_.add_dense_matrix(nodalAtomicMassWeightedDisplacement, 174 "NodalAtomicMassWeightedDisplacement"); 175 176 // nodal displacement derived only from atoms 177 AtfShapeFunctionMdProjection * nodalAtomicDisplacement = new AtfShapeFunctionMdProjection(this, 178 nodalAtomicMassWeightedDisplacement, 179 VELOCITY); 180 interscaleManager_.add_dense_matrix(nodalAtomicDisplacement, 181 "NodalAtomicDisplacement"); 182 } 183 184 // always need kinetic energy 185 FtaShapeFunctionProlongation * atomicMeanVelocity = new FtaShapeFunctionProlongation(this,&fields_[VELOCITY],shpFcn_); 186 interscaleManager_.add_per_atom_quantity(atomicMeanVelocity, 187 "AtomicMeanVelocity"); 188 AtomicEnergyForTemperature * atomicTwiceKineticEnergy = new TwiceFluctuatingKineticEnergy(this); 189 AtomicEnergyForTemperature * atomEnergyForTemperature = NULL; 190 191 // Appropriate per-atom quantity based on desired temperature definition 192 if (temperatureDef_==KINETIC) { 193 atomEnergyForTemperature = atomicTwiceKineticEnergy; 194 } 195 else if (temperatureDef_==TOTAL) { 196 if (timeIntegrators_[TEMPERATURE]->time_integration_type() != TimeIntegrator::FRACTIONAL_STEP) 197 throw ATC_Error("ATC_CouplingMomentumEnergy:construct_transfers() on the fractional step time integrator can be used with non-kinetic defitions of the temperature"); 198 199 // kinetic energy 200 interscaleManager_.add_per_atom_quantity(atomicTwiceKineticEnergy, 201 "AtomicTwiceKineticEnergy"); 202 203 // atomic potential energy 204 ComputedAtomQuantity * atomicPotentialEnergy = new ComputedAtomQuantity(this,lammpsInterface_->compute_pe_name(), 205 1./(lammpsInterface_->mvv2e())); 206 interscaleManager_.add_per_atom_quantity(atomicPotentialEnergy, 207 "AtomicPotentialEnergy"); 208 209 // reference potential energy 210 AtcAtomQuantity<double> * atomicReferencePotential; 211 if (!initialized_) { 212 atomicReferencePotential = new AtcAtomQuantity<double>(this); 213 interscaleManager_.add_per_atom_quantity(atomicReferencePotential, 214 "AtomicReferencePotential"); 215 atomicReferencePotential->set_memory_type(PERSISTENT); 216 } 217 else { 218 atomicReferencePotential = static_cast<AtcAtomQuantity<double> * >(interscaleManager_.per_atom_quantity("AtomicReferencePotential")); 219 } 220 nodalRefPotentialEnergy_ = new AtfShapeFunctionRestriction(this, 221 atomicReferencePotential, 222 shpFcn_); 223 interscaleManager_.add_dense_matrix(nodalRefPotentialEnergy_, 224 "NodalAtomicReferencePotential"); 225 226 // fluctuating potential energy 227 AtomicEnergyForTemperature * atomicFluctuatingPotentialEnergy = new FluctuatingPotentialEnergy(this, 228 atomicPotentialEnergy, 229 atomicReferencePotential); 230 interscaleManager_.add_per_atom_quantity(atomicFluctuatingPotentialEnergy, 231 "AtomicFluctuatingPotentialEnergy"); 232 233 // atomic total energy 234 atomEnergyForTemperature = new MixedKePeEnergy(this,1,1); 235 236 // kinetic temperature measure for post-processing 237 // nodal restriction of the atomic energy quantity for the temperature definition 238 AtfShapeFunctionRestriction * nodalAtomicTwiceKineticEnergy = new AtfShapeFunctionRestriction(this, 239 atomicTwiceKineticEnergy, 240 shpFcn_); 241 interscaleManager_.add_dense_matrix(nodalAtomicTwiceKineticEnergy, 242 "NodalAtomicTwiceKineticEnergy"); 243 nodalAtomicKineticTemperature_ = new AtfShapeFunctionMdProjection(this, 244 nodalAtomicTwiceKineticEnergy, 245 TEMPERATURE); 246 interscaleManager_.add_dense_matrix(nodalAtomicKineticTemperature_, 247 "NodalAtomicKineticTemperature"); 248 249 // potential temperature measure for post-processing (must multiply by 2 for configurational temperature 250 // nodal restriction of the atomic energy quantity for the temperature definition 251 AtfShapeFunctionRestriction * nodalAtomicFluctuatingPotentialEnergy = new AtfShapeFunctionRestriction(this, 252 atomicFluctuatingPotentialEnergy, 253 shpFcn_); 254 interscaleManager_.add_dense_matrix(nodalAtomicFluctuatingPotentialEnergy, 255 "NodalAtomicFluctuatingPotentialEnergy"); 256 nodalAtomicConfigurationalTemperature_ = new AtfShapeFunctionMdProjection(this, 257 nodalAtomicFluctuatingPotentialEnergy, 258 TEMPERATURE); 259 interscaleManager_.add_dense_matrix(nodalAtomicConfigurationalTemperature_, 260 "NodalAtomicConfigurationalTemperature"); 261 } 262 263 // register the per-atom quantity for the temperature definition 264 interscaleManager_.add_per_atom_quantity(atomEnergyForTemperature, 265 "AtomicEnergyForTemperature"); 266 267 // nodal restriction of the atomic energy quantity for the temperature definition 268 AtfShapeFunctionRestriction * nodalAtomicEnergy = new AtfShapeFunctionRestriction(this, 269 atomEnergyForTemperature, 270 shpFcn_); 271 interscaleManager_.add_dense_matrix(nodalAtomicEnergy, 272 "NodalAtomicEnergy"); 273 274 // nodal atomic temperature field 275 276 AtfShapeFunctionMdProjection * nodalAtomicTemperature = new AtfShapeFunctionMdProjection(this, 277 nodalAtomicEnergy, 278 TEMPERATURE); 279 interscaleManager_.add_dense_matrix(nodalAtomicTemperature, 280 "NodalAtomicTemperature"); 281 282 for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) { 283 (_tiIt_->second)->construct_transfers(); 284 } 285 atomicRegulator_->construct_transfers(); 286 } 287 288 //--------------------------------------------------------- 289 // init_filter 290 // sets up the time filtering operations in all objects 291 //--------------------------------------------------------- init_filter()292 void ATC_CouplingMomentumEnergy::init_filter() 293 { 294 if (timeIntegrators_[TEMPERATURE]->time_integration_type() != TimeIntegrator::FRACTIONAL_STEP) { 295 throw ATC_Error("ATC_CouplingMomentumEnergy::initialize - method only valid with fractional step time integration"); 296 } 297 298 299 ATC_Coupling::init_filter(); 300 301 302 303 304 if (timeFilterManager_.end_equilibrate() && equilibriumStart_) { 305 if (atomicRegulator_->coupling_mode(VELOCITY)==AtomicRegulator::FLUX || atomicRegulator_->coupling_mode(VELOCITY)==AtomicRegulator::GHOST_FLUX) 306 // nothing needed in other cases since kinetostat force is balanced by boundary flux in FE equations 307 atomicRegulator_->reset_lambda_contribution(nodalAtomicFieldsRoc_[VELOCITY].quantity(),VELOCITY); 308 309 DENS_MAT powerMat(-1.*(nodalAtomicFields_[TEMPERATURE].quantity())); 310 atomicRegulator_->reset_lambda_contribution(powerMat,TEMPERATURE); 311 } 312 } 313 314 //-------------------------------------------------------- 315 // modify 316 // parses inputs and modifies state of the filter 317 //-------------------------------------------------------- modify(int narg,char ** arg)318 bool ATC_CouplingMomentumEnergy::modify(int narg, char **arg) 319 { 320 bool foundMatch = false; 321 int argIndex = 0; 322 return foundMatch; 323 } 324 325 //-------------------------------------------------------------------- 326 // compute_scalar : added energy 327 //-------------------------------------------------------------------- compute_scalar(void)328 double ATC_CouplingMomentumEnergy::compute_scalar(void) 329 { 330 double energy = 0.0; 331 energy += extrinsicModelManager_.compute_scalar(); 332 return energy; 333 } 334 335 //-------------------------------------------------------------------- 336 // total kinetic energy 337 //-------------------------------------------------------------------- kinetic_energy(void)338 double ATC_CouplingMomentumEnergy::kinetic_energy(void) 339 { 340 const MATRIX & M = massMats_[VELOCITY].quantity(); 341 342 const DENS_MAT & velocity(fields_[VELOCITY].quantity()); 343 double mvv2e = lammpsInterface_->mvv2e(); 344 double kineticEnergy = 0; 345 DENS_VEC velocitySquared(nNodes_); 346 for (int i = 0; i < nNodes_; i++) 347 for (int j = 0; j < nsd_; j++) 348 velocitySquared(i) += velocity(i,j)*velocity(i,j); 349 kineticEnergy = (M*velocitySquared).sum(); 350 kineticEnergy *= mvv2e; // convert to LAMMPS units 351 return kineticEnergy; 352 } 353 //-------------------------------------------------------------------- 354 // total potential energy 355 //-------------------------------------------------------------------- potential_energy(void)356 double ATC_CouplingMomentumEnergy::potential_energy(void) 357 { 358 Array<FieldName> mask(1); 359 mask(0) = VELOCITY; 360 FIELD_MATS energy; 361 feEngine_->compute_energy(mask, 362 fields_, 363 physicsModel_, 364 elementToMaterialMap_, 365 energy, 366 &(elementMask_->quantity())); 367 double potentialEnergy = energy[VELOCITY].col_sum(); 368 double mvv2e = lammpsInterface_->mvv2e(); 369 potentialEnergy *= mvv2e; // convert to LAMMPS units 370 return potentialEnergy-refPE_; 371 } 372 373 //-------------------------------------------------------------------- 374 // compute_vector 375 //-------------------------------------------------------------------- 376 // this is for direct output to lammps thermo compute_vector(int n)377 double ATC_CouplingMomentumEnergy::compute_vector(int n) 378 { 379 // output[1] = total coarse scale kinetic energy 380 // output[2] = total coarse scale potential energy 381 // output[3] = total coarse scale energy 382 // output[4] = total coarse scale thermal energy 383 // output[5] = average temperature 384 385 double mvv2e = lammpsInterface_->mvv2e(); // convert to lammps energy units 386 387 if (n == 0) { 388 return kinetic_energy(); 389 } 390 else if (n == 1) { 391 return potential_energy(); 392 } 393 else if (n == 2) { 394 return kinetic_energy()+potential_energy(); 395 } 396 else if (n == 4) { 397 Array<FieldName> mask(1); 398 FIELD_MATS energy; 399 mask(0) = TEMPERATURE; 400 401 feEngine_->compute_energy(mask, 402 fields_, 403 physicsModel_, 404 elementToMaterialMap_, 405 energy, 406 &(elementMask_->quantity())); 407 408 double phononEnergy = mvv2e * energy[TEMPERATURE].col_sum(); 409 return phononEnergy; 410 } 411 else if (n == 5) { 412 double aveT = (fields_[TEMPERATURE].quantity()).col_sum()/nNodes_; 413 return aveT; 414 } 415 else if (n > 5) { 416 double extrinsicValue = extrinsicModelManager_.compute_vector(n); 417 return extrinsicValue; 418 } 419 420 return 0.; 421 422 } 423 424 //-------------------------------------------------------------------- 425 // output 426 //-------------------------------------------------------------------- output()427 void ATC_CouplingMomentumEnergy::output() 428 { 429 if (output_now()) { 430 feEngine_->departition_mesh(); 431 // avoid possible mpi calls 432 if (nodalAtomicKineticTemperature_) 433 _keTemp_ = nodalAtomicKineticTemperature_->quantity(); 434 if (nodalAtomicConfigurationalTemperature_) 435 _peTemp_ = nodalAtomicConfigurationalTemperature_->quantity(); 436 437 OUTPUT_LIST outputData; 438 439 // base class output 440 ATC_Method::output(); 441 442 // push atc fields time integrator modifies into output arrays 443 for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) { 444 (_tiIt_->second)->post_process(); 445 } 446 447 // auxilliary data 448 449 for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) { 450 (_tiIt_->second)->output(outputData); 451 } 452 atomicRegulator_->output(outputData); 453 extrinsicModelManager_.output(outputData); 454 455 DENS_MAT & velocity(nodalAtomicFields_[VELOCITY].set_quantity()); 456 DENS_MAT & rhs(rhs_[VELOCITY].set_quantity()); 457 DENS_MAT & temperature(nodalAtomicFields_[TEMPERATURE].set_quantity()); 458 DENS_MAT & dotTemperature(dot_fields_[TEMPERATURE].set_quantity()); 459 DENS_MAT & ddotTemperature(ddot_fields_[TEMPERATURE].set_quantity()); 460 DENS_MAT & rocTemperature(nodalAtomicFieldsRoc_[TEMPERATURE].set_quantity()); 461 DENS_MAT & fePower(rhs_[TEMPERATURE].set_quantity()); 462 if (lammpsInterface_->rank_zero()) { 463 // global data 464 double T_mean = (fields_[TEMPERATURE].quantity()).col_sum(0)/nNodes_; 465 feEngine_->add_global("temperature_mean", T_mean); 466 double T_stddev = (fields_[TEMPERATURE].quantity()).col_stdev(0); 467 feEngine_->add_global("temperature_std_dev", T_stddev); 468 double Ta_mean = (nodalAtomicFields_[TEMPERATURE].quantity()).col_sum(0)/nNodes_; 469 feEngine_->add_global("atomic_temperature_mean", Ta_mean); 470 double Ta_stddev = (nodalAtomicFields_[TEMPERATURE].quantity()).col_stdev(0); 471 feEngine_->add_global("atomic_temperature_std_dev", Ta_stddev); 472 473 // different temperature measures, if appropriate 474 if (nodalAtomicKineticTemperature_) 475 outputData["kinetic_temperature"] = & _keTemp_; 476 477 if (nodalAtomicConfigurationalTemperature_) { 478 _peTemp_ *= 2; // account for full temperature 479 outputData["configurational_temperature"] = & _peTemp_; 480 } 481 482 // mesh data 483 outputData["NodalAtomicVelocity"] = &velocity; 484 outputData["FE_Force"] = &rhs; 485 if (trackDisplacement_) 486 outputData["NodalAtomicDisplacement"] = & nodalAtomicFields_[DISPLACEMENT].set_quantity(); 487 outputData["NodalAtomicTemperature"] = &temperature; 488 outputData["dot_temperature"] = &dotTemperature; 489 outputData["ddot_temperature"] = &ddotTemperature; 490 outputData["NodalAtomicPower"] = &rocTemperature; 491 outputData["fePower"] = &fePower; 492 493 feEngine_->write_data(output_index(), fields_, & outputData); 494 } 495 496 // hence propagation is performed on proc 0 but not others. 497 // The real fix is to have const data in the output list 498 // force optional variables to reset to keep in sync 499 if (trackDisplacement_) { 500 nodalAtomicFields_[DISPLACEMENT].force_reset(); 501 } 502 fields_[VELOCITY].propagate_reset(); 503 504 feEngine_->partition_mesh(); 505 } 506 } 507 508 }; 509