1 #include "Thermostat.h" 2 #include "ATC_Coupling.h" 3 #include "ATC_Error.h" 4 #include "PrescribedDataManager.h" 5 #include "ThermalTimeIntegrator.h" 6 #include "TransferOperator.h" 7 8 using namespace std; 9 10 namespace ATC { 11 12 //-------------------------------------------------------- 13 //-------------------------------------------------------- 14 // Class Thermostat 15 //-------------------------------------------------------- 16 //-------------------------------------------------------- 17 18 //-------------------------------------------------------- 19 // Constructor 20 //-------------------------------------------------------- Thermostat(ATC_Coupling * atc,const string & regulatorPrefix)21 Thermostat::Thermostat(ATC_Coupling * atc, 22 const string & regulatorPrefix) : 23 AtomicRegulator(atc,regulatorPrefix), 24 lambdaMaxIterations_(myLambdaMaxIterations) 25 { 26 // nothing to do 27 } 28 29 //-------------------------------------------------------- 30 // modify: 31 // parses and adjusts thermostat state based on 32 // user input, in the style of LAMMPS user input 33 //-------------------------------------------------------- modify(int narg,char ** arg)34 bool Thermostat::modify(int narg, char **arg) 35 { 36 bool foundMatch = false; 37 38 int argIndex = 0; 39 if (strcmp(arg[argIndex],"thermal")==0) { 40 argIndex++; 41 42 // thermostat type 43 /*! \page man_control_thermal fix_modify AtC control thermal 44 \section syntax 45 fix_modify AtC control thermal <control_type> <optional_args> 46 - control_type (string) = none | rescale | hoover | flux\n 47 48 fix_modify AtC control thermal rescale <frequency> \n 49 - frequency (int) = time step frequency for applying velocity rescaling \n 50 51 fix_modify AtC control thermal hoover \n 52 53 fix_modify AtC control thermal flux <boundary_integration_type(optional)> <face_set_id(optional)>\n 54 - boundary_integration_type (string) = faceset | interpolate\n 55 - face_set_id (string), optional = id of boundary face set, if not specified 56 (or not possible when the atomic domain does not line up with 57 mesh boundaries) defaults to an atomic-quadrature approximate 58 evaulation, does not work with interpolate\n 59 \section examples 60 <TT> fix_modify AtC control thermal none </TT> \n 61 <TT> fix_modify AtC control thermal rescale 10 </TT> \n 62 <TT> fix_modify AtC control thermal hoover </TT> \n 63 <TT> fix_modify AtC control thermal flux </TT> \n 64 <TT> fix_modify AtC control thermal flux faceset bndy_faces </TT> \n 65 \section description 66 Sets the energy exchange mechansim from the finite elements to the atoms, managed through a control algorithm. Rescale computes a scale factor for each atom to match the finite element temperature. Hoover is a Gaussian least-constraint isokinetic thermostat enforces that the nodal restricted atomic temperature matches the finite element temperature. Flux is a similar mode, but rather adds energy to the atoms based on conservation of energy. Hoover and flux allows the prescription of sources or fixed temperatures on the atoms. 67 \section restrictions 68 only for be used with specific transfers : 69 thermal (rescale, hoover, flux), two_temperature (flux) \n 70 rescale not valid with time filtering activated 71 \section related 72 \section default 73 none\n 74 rescale frequency is 1\n 75 flux boundary_integration_type is interpolate 76 */ 77 if (strcmp(arg[argIndex],"none")==0) { // restore defaults 78 regulatorTarget_ = NONE; 79 couplingMode_ = UNCOUPLED; 80 howOften_ = 1; 81 boundaryIntegrationType_ = NO_QUADRATURE; 82 foundMatch = true; 83 } 84 else if (strcmp(arg[argIndex],"rescale")==0) { 85 argIndex++; 86 howOften_ = atoi(arg[argIndex]); 87 if (howOften_ < 1) { 88 throw ATC_Error("Bad rescaling thermostat frequency"); 89 } 90 else { 91 regulatorTarget_ = FIELD; 92 couplingMode_ = UNCOUPLED; 93 boundaryIntegrationType_ = NO_QUADRATURE; 94 foundMatch = true; 95 } 96 } 97 else if (strcmp(arg[argIndex],"hoover")==0) { 98 regulatorTarget_ = DYNAMICS; 99 couplingMode_ = FIXED; 100 howOften_ = 1; 101 boundaryIntegrationType_ = NO_QUADRATURE; 102 foundMatch = true; 103 } 104 else if (strcmp(arg[argIndex],"flux")==0) { 105 regulatorTarget_ = DYNAMICS; 106 couplingMode_ = FLUX; 107 howOften_ = 1; 108 argIndex++; 109 110 boundaryIntegrationType_ = atc_->parse_boundary_integration(narg-argIndex,&arg[argIndex],boundaryFaceSet_); 111 foundMatch = true; 112 } 113 // set parameters for numerical matrix solutions unique to this thermostat 114 /*! \page man_control_thermal_correction_max_iterations fix_modify AtC control thermal correction_max_iterations 115 \section syntax 116 fix_modify AtC control thermal correction_max_iterations <max_iterations> 117 - max_iterations (int) = maximum number of iterations that will be used by iterative matrix solvers\n 118 119 \section examples 120 <TT> fix_modify AtC control thermal correction_max_iterations 10 </TT> \n 121 \section description 122 Sets the maximum number of iterations to compute the 2nd order in time correction term for lambda with the fractional step method. The method uses the same tolerance as the controller's matrix solver. 123 \section restrictions 124 only for use with thermal physics using the fractional step method. 125 \section related 126 \section default 127 correction_max_iterations is 20 128 */ 129 else if (strcmp(arg[argIndex],"correction_max_iterations")==0) { 130 argIndex++; 131 lambdaMaxIterations_ = atoi(arg[argIndex]); 132 if (lambdaMaxIterations_ < 1) { 133 throw ATC_Error("Bad correction maximum iteration count"); 134 } 135 foundMatch = true; 136 } 137 } 138 139 if (!foundMatch) 140 foundMatch = AtomicRegulator::modify(narg,arg); 141 if (foundMatch) 142 needReset_ = true; 143 return foundMatch; 144 } 145 146 //-------------------------------------------------------- 147 // reset_lambda_contribution: 148 // resets the thermostat generated power to a 149 // prescribed value 150 //-------------------------------------------------------- reset_lambda_contribution(const DENS_MAT & target)151 void Thermostat::reset_lambda_contribution(const DENS_MAT & target) 152 { 153 DENS_MAN * lambdaPowerFiltered = regulator_data("LambdaPowerFiltered",1); 154 *lambdaPowerFiltered = target; 155 } 156 157 //-------------------------------------------------------- 158 // construct_methods: 159 // instantiations desired regulator method(s) 160 161 // dependence, but in general there is also a 162 // time integrator dependence. In general the 163 // precedence order is: 164 // time filter -> time integrator -> thermostat 165 // In the future this may need to be added if 166 // different types of time integrators can be 167 // specified. 168 //-------------------------------------------------------- construct_methods()169 void Thermostat::construct_methods() 170 { 171 // get data associated with stages 1 & 2 of ATC_Method::initialize 172 AtomicRegulator::construct_methods(); 173 174 if (atc_->reset_methods()) { 175 // eliminate existing methods 176 delete_method(); 177 178 // update time filter 179 TimeIntegrator::TimeIntegrationType myIntegrationType = (atc_->time_integrator(TEMPERATURE))->time_integration_type(); 180 TimeFilterManager * timeFilterManager = atc_->time_filter_manager(); 181 if (timeFilterManager->need_reset() ) { 182 if (myIntegrationType == TimeIntegrator::GEAR) { 183 timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT); 184 } 185 else if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) { 186 timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT_IMPLICIT); 187 } 188 } 189 190 if (timeFilterManager->filter_dynamics()) { 191 switch (regulatorTarget_) { 192 case NONE: { 193 regulatorMethod_ = new RegulatorMethod(this); 194 break; 195 } 196 case FIELD: { // error check, rescale and filtering not supported together 197 throw ATC_Error("Cannot use rescaling thermostat with time filtering"); 198 break; 199 } 200 case DYNAMICS: { 201 switch (couplingMode_) { 202 case FIXED: { 203 if (use_lumped_lambda_solve()) { 204 throw ATC_Error("Thermostat:construct_methods - lumped lambda solve cannot be used with Hoover thermostats"); 205 } 206 if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) { 207 if (md_flux_nodes(TEMPERATURE)) { 208 if (!md_fixed_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) { 209 // there are fluxes but no fixed or coupled nodes 210 regulatorMethod_ = new ThermostatIntegratorFluxFiltered(this,lambdaMaxIterations_); 211 } 212 else { 213 // there are both fixed and flux nodes 214 regulatorMethod_ = new ThermostatFluxFixedFiltered(this,lambdaMaxIterations_); 215 } 216 } 217 else { 218 // there are only fixed nodes 219 regulatorMethod_ = new ThermostatIntegratorFixedFiltered(this,lambdaMaxIterations_); 220 } 221 } 222 else { 223 regulatorMethod_ = new ThermostatHooverVerletFiltered(this); 224 } 225 break; 226 } 227 case FLUX: { 228 if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) { 229 if (use_lumped_lambda_solve()) { 230 throw ATC_Error("Thermostat:construct_methods - lumped lambda solve has been depricated for fractional step thermostats"); 231 } 232 if (md_fixed_nodes(TEMPERATURE)) { 233 if (!md_flux_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) { 234 // there are fixed nodes but no fluxes 235 regulatorMethod_ = new ThermostatIntegratorFixedFiltered(this,lambdaMaxIterations_); 236 } 237 else { 238 // there are both fixed and flux nodes 239 regulatorMethod_ = new ThermostatFluxFixedFiltered(this,lambdaMaxIterations_); 240 } 241 } 242 else { 243 // there are only flux nodes 244 regulatorMethod_ = new ThermostatIntegratorFluxFiltered(this,lambdaMaxIterations_); 245 } 246 } 247 else { 248 if (use_localized_lambda()) { 249 if (!((atc_->prescribed_data_manager())->no_fluxes(TEMPERATURE)) && 250 atc_->boundary_integration_type() != NO_QUADRATURE) { 251 throw ATC_Error("Cannot use flux coupling with localized lambda"); 252 } 253 } 254 regulatorMethod_ = new ThermostatPowerVerletFiltered(this); 255 } 256 break; 257 } 258 default: 259 throw ATC_Error("Unknown coupling mode in Thermostat::initialize"); 260 } 261 break; 262 } 263 default: 264 throw ATC_Error("Unknown thermostat type in Thermostat::initialize"); 265 } 266 } 267 else { 268 switch (regulatorTarget_) { 269 case NONE: { 270 regulatorMethod_ = new RegulatorMethod(this); 271 break; 272 } 273 case FIELD: { 274 if (atc_->temperature_def()==KINETIC) 275 regulatorMethod_ = new ThermostatRescale(this); 276 else if (atc_->temperature_def()==TOTAL) 277 regulatorMethod_ = new ThermostatRescaleMixedKePe(this); 278 else 279 throw ATC_Error("Unknown temperature definition"); 280 break; 281 } 282 case DYNAMICS: { 283 switch (couplingMode_) { 284 case FIXED: { 285 if (use_lumped_lambda_solve()) { 286 throw ATC_Error("Thermostat:construct_methods - lumped lambda solve cannot be used with Hoover thermostats"); 287 } 288 if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) { 289 if (md_flux_nodes(TEMPERATURE)) { 290 if (!md_fixed_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) { 291 // there are fluxes but no fixed or coupled nodes 292 regulatorMethod_ = new ThermostatIntegratorFlux(this,lambdaMaxIterations_); 293 } 294 else { 295 // there are both fixed and flux nodes 296 regulatorMethod_ = new ThermostatFluxFixed(this,lambdaMaxIterations_); 297 } 298 } 299 else { 300 // there are only fixed nodes 301 regulatorMethod_ = new ThermostatIntegratorFixed(this,lambdaMaxIterations_); 302 } 303 } 304 else { 305 regulatorMethod_ = new ThermostatHooverVerlet(this); 306 } 307 break; 308 } 309 case FLUX: { 310 if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) { 311 if (use_lumped_lambda_solve()) { 312 throw ATC_Error("Thermostat:construct_methods - lumped lambda solve has been depricated for fractional step thermostats"); 313 } 314 if (md_fixed_nodes(TEMPERATURE)) { 315 if (!md_flux_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) { 316 // there are fixed nodes but no fluxes 317 regulatorMethod_ = new ThermostatIntegratorFixed(this,lambdaMaxIterations_); 318 } 319 else { 320 // there are both fixed and flux nodes 321 regulatorMethod_ = new ThermostatFluxFixed(this,lambdaMaxIterations_); 322 } 323 } 324 else { 325 // there are only flux nodes 326 regulatorMethod_ = new ThermostatIntegratorFlux(this,lambdaMaxIterations_); 327 } 328 } 329 else { 330 if (use_localized_lambda()) { 331 if (!((atc_->prescribed_data_manager())->no_fluxes(TEMPERATURE)) && 332 atc_->boundary_integration_type() != NO_QUADRATURE) { 333 throw ATC_Error("Cannot use flux coupling with localized lambda"); 334 } 335 } 336 regulatorMethod_ = new ThermostatPowerVerlet(this); 337 } 338 break; 339 } 340 default: 341 throw ATC_Error("Unknown coupling mode in Thermostat::initialize"); 342 } 343 break; 344 } 345 default: 346 throw ATC_Error("Unknown thermostat target in Thermostat::initialize"); 347 } 348 } 349 350 AtomicRegulator::reset_method(); 351 } 352 else { 353 set_all_data_to_used(); 354 } 355 } 356 357 //-------------------------------------------------------- 358 //-------------------------------------------------------- 359 // Class ThermostatShapeFunction 360 //-------------------------------------------------------- 361 //-------------------------------------------------------- 362 363 //-------------------------------------------------------- 364 // Constructor 365 //-------------------------------------------------------- ThermostatShapeFunction(AtomicRegulator * thermostat,const string & regulatorPrefix)366 ThermostatShapeFunction::ThermostatShapeFunction(AtomicRegulator * thermostat, 367 const string & regulatorPrefix) : 368 RegulatorShapeFunction(thermostat,regulatorPrefix), 369 mdMassMatrix_(atc_->set_mass_mat_md(TEMPERATURE)), 370 atomVelocities_(nullptr) 371 { 372 fieldMask_(TEMPERATURE,FLUX) = true; 373 lambda_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1); // data associated with stage 3 in ATC_Method::initialize 374 } 375 376 //-------------------------------------------------------- 377 // constructor_transfers 378 // instantiates or obtains all dependency managed data 379 //-------------------------------------------------------- construct_transfers()380 void ThermostatShapeFunction::construct_transfers() 381 { 382 InterscaleManager & interscaleManager(atc_->interscale_manager()); 383 RegulatorShapeFunction::construct_transfers(); 384 385 // get atom velocity data from manager 386 atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY); 387 388 // construct lambda evaluated at atom locations 389 atomLambdas_ = new FtaShapeFunctionProlongation(atc_, 390 lambda_, 391 interscaleManager.per_atom_sparse_matrix("Interpolant")); 392 interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy"); 393 } 394 395 //--------------------------------------------------------- 396 // set_weights: 397 // set the diagonal weighting matrix to be the atomic 398 // temperatures 399 //--------------------------------------------------------- set_weights()400 void ThermostatShapeFunction::set_weights() 401 { 402 if (this->use_local_shape_functions()) { 403 VelocitySquaredMapped * myWeights = new VelocitySquaredMapped(atc_,lambdaAtomMap_); 404 weights_ = myWeights; 405 (atc_->interscale_manager()).add_per_atom_quantity(myWeights, 406 regulatorPrefix_+"AtomVelocitySquaredMapped"); 407 } 408 else { 409 VelocitySquared * myWeights = new VelocitySquared(atc_); 410 weights_ = myWeights; 411 (atc_->interscale_manager()).add_per_atom_quantity(myWeights, 412 regulatorPrefix_+"AtomVelocitySquared"); 413 } 414 } 415 416 //-------------------------------------------------------- 417 //-------------------------------------------------------- 418 // Class ThermostatRescale 419 //-------------------------------------------------------- 420 //-------------------------------------------------------- 421 422 //-------------------------------------------------------- 423 // Constructor 424 //-------------------------------------------------------- ThermostatRescale(AtomicRegulator * thermostat)425 ThermostatRescale::ThermostatRescale(AtomicRegulator * thermostat) : 426 ThermostatShapeFunction(thermostat), 427 nodalTemperature_(atc_->field(TEMPERATURE)), 428 atomVelocityRescalings_(nullptr) 429 { 430 // do nothing 431 } 432 433 //-------------------------------------------------------- 434 // constructor_transfers 435 // instantiates or obtains all dependency managed data 436 //-------------------------------------------------------- construct_transfers()437 void ThermostatRescale::construct_transfers() 438 { 439 InterscaleManager & interscaleManager(atc_->interscale_manager()); 440 441 // set up node mappings 442 create_node_maps(); 443 444 // set up data for linear solver 445 shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_); 446 (atc_->interscale_manager()).add_per_atom_sparse_matrix(shapeFunctionMatrix_, 447 regulatorPrefix_+"LambdaCouplingMatrixEnergy"); 448 linearSolverType_ = AtomicRegulator::CG_SOLVE; 449 450 // base class transfers 451 ThermostatShapeFunction::construct_transfers(); 452 453 // velocity rescaling factor 454 atomVelocityRescalings_ = new AtomicVelocityRescaleFactor(atc_,atomLambdas_); 455 interscaleManager.add_per_atom_quantity(atomVelocityRescalings_, 456 regulatorPrefix_+"AtomVelocityRescaling"); 457 } 458 459 //--------------------------------------------------------- 460 // set_weights: 461 // set the diagonal weighting matrix to be the atomic 462 // temperatures 463 //--------------------------------------------------------- set_weights()464 void ThermostatRescale::set_weights() 465 { 466 weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicEnergyForTemperature"); 467 } 468 469 //-------------------------------------------------------- 470 // apply_post_corrector: 471 // apply the thermostat in the post corrector phase 472 //-------------------------------------------------------- apply_post_corrector(double dt)473 void ThermostatRescale::apply_post_corrector(double dt) 474 { 475 compute_thermostat(dt); 476 477 // application of rescaling lambda due 478 apply_to_atoms(atomVelocities_); 479 } 480 481 //-------------------------------------------------------- 482 // compute_thermostat 483 // manages the solution of the 484 // thermostat equations and variables 485 //-------------------------------------------------------- compute_thermostat(double)486 void ThermostatRescale::compute_thermostat(double /* dt */) 487 { 488 // compute right-hand side 489 this->set_rhs(_rhs_); 490 491 // solve equations 492 solve_for_lambda(_rhs_,lambda_->set_quantity()); 493 } 494 495 //-------------------------------------------------------- 496 // set_rhs: 497 // constructs the RHS vector with the target 498 // temperature 499 //-------------------------------------------------------- set_rhs(DENS_MAT & rhs)500 void ThermostatRescale::set_rhs(DENS_MAT & rhs) 501 { 502 rhs = mdMassMatrix_.quantity()*nodalTemperature_.quantity(); 503 } 504 505 //-------------------------------------------------------- 506 // apply_lambda_to_atoms: 507 // applies the velocity rescale with an existing lambda 508 // note oldAtomicQuantity and dt are not used 509 //-------------------------------------------------------- apply_to_atoms(PerAtomQuantity<double> * atomVelocities)510 void ThermostatRescale::apply_to_atoms(PerAtomQuantity<double> * atomVelocities) 511 { 512 *atomVelocities *= atomVelocityRescalings_->quantity(); 513 } 514 515 //-------------------------------------------------------- 516 // output: 517 // adds all relevant output to outputData 518 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)519 void ThermostatRescale::output(OUTPUT_LIST & outputData) 520 { 521 DENS_MAT & lambda(lambda_->set_quantity()); 522 if ((atc_->lammps_interface())->rank_zero()) { 523 outputData["LambdaEnergy"] = λ 524 } 525 } 526 527 //-------------------------------------------------------- 528 //-------------------------------------------------------- 529 // Class ThermostatRescaleMixedKePe 530 //-------------------------------------------------------- 531 //-------------------------------------------------------- 532 533 //-------------------------------------------------------- 534 // Constructor 535 //-------------------------------------------------------- ThermostatRescaleMixedKePe(AtomicRegulator * thermostat)536 ThermostatRescaleMixedKePe::ThermostatRescaleMixedKePe(AtomicRegulator * thermostat) : 537 ThermostatRescale(thermostat), 538 nodalAtomicFluctuatingPotentialEnergy_(nullptr) 539 { 540 // do nothing 541 } 542 543 //-------------------------------------------------------- 544 // constructor_transfers 545 // instantiates or obtains all dependency managed data 546 //-------------------------------------------------------- construct_transfers()547 void ThermostatRescaleMixedKePe::construct_transfers() 548 { 549 ThermostatRescale::construct_transfers(); 550 InterscaleManager & interscaleManager(atc_->interscale_manager()); 551 552 // get fluctuating PE at nodes 553 nodalAtomicFluctuatingPotentialEnergy_ = 554 interscaleManager.dense_matrix("NodalAtomicFluctuatingPotentialEnergy"); 555 } 556 557 //--------------------------------------------------------- 558 // set_weights: 559 // set the diagonal weighting matrix to be the atomic 560 // temperatures 561 //--------------------------------------------------------- set_weights()562 void ThermostatRescaleMixedKePe::set_weights() 563 { 564 weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicTwiceKineticEnergy"); 565 } 566 567 //-------------------------------------------------------- 568 // initialize 569 // initializes all method data 570 //-------------------------------------------------------- initialize()571 void ThermostatRescaleMixedKePe::initialize() 572 { 573 ThermostatRescale::initialize(); 574 InterscaleManager & interscaleManager(atc_->interscale_manager()); 575 576 // multipliers for KE and PE 577 AtomicEnergyForTemperature * atomEnergyForTemperature = 578 static_cast<AtomicEnergyForTemperature * >(interscaleManager.per_atom_quantity("AtomicEnergyForTemperature")); 579 keMultiplier_ = atomEnergyForTemperature->kinetic_energy_multiplier(); 580 peMultiplier_ = 2. - keMultiplier_; 581 keMultiplier_ /= 2.; // account for use of 2 X KE in matrix equation 582 } 583 584 //-------------------------------------------------------- 585 // set_rhs: 586 // accounts for potential energy contribution to 587 // definition of atomic temperature 588 //-------------------------------------------------------- set_rhs(DENS_MAT & rhs)589 void ThermostatRescaleMixedKePe::set_rhs(DENS_MAT & rhs) 590 { 591 ThermostatRescale::set_rhs(rhs); 592 rhs -= peMultiplier_*(nodalAtomicFluctuatingPotentialEnergy_->quantity()); 593 rhs /= keMultiplier_; 594 } 595 596 //-------------------------------------------------------- 597 //-------------------------------------------------------- 598 // Class ThermostatFsSolver 599 //-------------------------------------------------------- 600 //-------------------------------------------------------- 601 602 //-------------------------------------------------------- 603 // Constructor 604 //-------------------------------------------------------- ThermostatFsSolver(AtomicRegulator * thermostat,int lambdaMaxIterations,const string & regulatorPrefix)605 ThermostatFsSolver::ThermostatFsSolver(AtomicRegulator * thermostat, 606 int lambdaMaxIterations, 607 const string & regulatorPrefix) : 608 RegulatorShapeFunction(thermostat,regulatorPrefix), 609 lambdaMaxIterations_(lambdaMaxIterations), 610 rhsLambdaSquared_(nullptr), 611 dtFactor_(1.) 612 { 613 fieldMask_(TEMPERATURE,FLUX) = true; 614 lambda_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1); // data associated with stage 3 in ATC_Method::initialize 615 } 616 617 //-------------------------------------------------------- 618 // initialize 619 // creates mapping from all nodes to those to which 620 // the thermostat applies 621 //-------------------------------------------------------- initialize()622 void ThermostatFsSolver::initialize() 623 { 624 RegulatorShapeFunction::initialize(); 625 626 rhsMap_.resize(overlapToNodeMap_->nRows(),1); 627 DENS_MAT rhsMapGlobal(nNodes_,1); 628 const set<int> & applicationNodes(applicationNodes_->quantity()); 629 for (int i = 0; i < nNodes_; i++) { 630 if (applicationNodes.find(i) != applicationNodes.end()) { 631 rhsMapGlobal(i,0) = 1.; 632 } 633 else { 634 rhsMapGlobal(i,0) = 0.; 635 } 636 } 637 map_unique_to_overlap(rhsMapGlobal,rhsMap_); 638 } 639 640 //--------------------------------------------------------- 641 // set_weights: 642 // set the diagonal weighting matrix to be the atomic 643 // temperatures 644 //--------------------------------------------------------- set_weights()645 void ThermostatFsSolver::set_weights() 646 { 647 if (this->use_local_shape_functions()) { 648 VelocitySquaredMapped * myWeights = new VelocitySquaredMapped(atc_,lambdaAtomMap_); 649 weights_ = myWeights; 650 (atc_->interscale_manager()).add_per_atom_quantity(myWeights, 651 regulatorPrefix_+"AtomVelocitySquaredMapped"); 652 } 653 else { 654 VelocitySquared * myWeights = new VelocitySquared(atc_); 655 weights_ = myWeights; 656 (atc_->interscale_manager()).add_per_atom_quantity(myWeights, 657 regulatorPrefix_+"AtomVelocitySquared"); 658 } 659 } 660 661 //-------------------------------------------------------- 662 // compute_lambda: 663 // solves linear system for lambda, if the 664 // bool is true it iterators to a non-linear solution 665 //-------------------------------------------------------- compute_lambda(const DENS_MAT & rhs,bool iterateSolution)666 void ThermostatFsSolver::compute_lambda(const DENS_MAT & rhs, 667 bool iterateSolution) 668 { 669 // solve linear system for lambda guess 670 DENS_MAT & lambda(lambda_->set_quantity()); 671 solve_for_lambda(rhs,lambda); 672 673 // iterate to solution 674 if (iterateSolution) { 675 iterate_lambda(rhs); 676 } 677 } 678 679 //-------------------------------------------------------- 680 // iterate_lambda: 681 // iteratively solves the equations for lambda 682 // for the higher order dt corrections, assuming 683 // an initial guess for lambda 684 //-------------------------------------------------------- iterate_lambda(const MATRIX & rhs)685 void ThermostatFsSolver::iterate_lambda(const MATRIX & rhs) 686 { 687 int nNodeOverlap = overlapToNodeMap_->nRows(); 688 DENS_VEC _lambdaOverlap_(nNodeOverlap); 689 DENS_MAT & lambda(lambda_->set_quantity()); 690 map_unique_to_overlap(lambda,_lambdaOverlap_); 691 double factor = 0.5*dtFactor_*atc_->dt(); 692 693 _lambdaOld_.resize(nNodes_,1); 694 _rhsOverlap_.resize(nNodeOverlap,1); 695 map_unique_to_overlap(rhs,_rhsOverlap_); 696 _rhsTotal_.resize(nNodeOverlap); 697 698 // solve assuming we get initial guess for lambda 699 double error(-1.); 700 for (int i = 0; i < lambdaMaxIterations_; ++i) { 701 _lambdaOld_ = lambda; 702 703 // solve the system with the new rhs 704 const DENS_MAT & rhsLambdaSquared(rhsLambdaSquared_->quantity()); 705 for (int i = 0; i < nNodeOverlap; i++) { 706 if (rhsMap_(i,0) == 1.) { 707 _rhsTotal_(i) = _rhsOverlap_(i,0) + factor*rhsLambdaSquared(i,0); 708 } 709 else { 710 _rhsTotal_(i) = 0.; 711 } 712 } 713 matrixSolver_->execute(_rhsTotal_,_lambdaOverlap_); 714 715 // check convergence 716 map_overlap_to_unique(_lambdaOverlap_,lambda); 717 lambda_->force_reset(); 718 DENS_MAT difference = lambda-_lambdaOld_; 719 error = difference.col_norm()/_lambdaOld_.col_norm(); 720 if (error < tolerance_) 721 break; 722 } 723 724 if (error >= tolerance_) { 725 stringstream message; 726 message << "WARNING: Iterative solve for lambda failed to converge after " << lambdaMaxIterations_ << " iterations, final tolerance was " << error << "\n"; 727 ATC::LammpsInterface::instance()->print_msg(message.str()); 728 } 729 } 730 731 //-------------------------------------------------------- 732 //-------------------------------------------------------- 733 // Class ThermostatGlcFs 734 //-------------------------------------------------------- 735 //-------------------------------------------------------- 736 737 //-------------------------------------------------------- 738 // Constructor 739 //-------------------------------------------------------- ThermostatGlcFs(AtomicRegulator * thermostat,int,const string & regulatorPrefix)740 ThermostatGlcFs::ThermostatGlcFs(AtomicRegulator * thermostat, 741 int /* lambdaMaxIterations */, 742 const string & regulatorPrefix) : 743 RegulatorMethod(thermostat,regulatorPrefix), 744 lambdaSolver_(nullptr), 745 mdMassMatrix_(atc_->set_mass_mat_md(TEMPERATURE)), 746 atomVelocities_(nullptr), 747 temperature_(atc_->field(TEMPERATURE)), 748 timeFilter_(atomicRegulator_->time_filter()), 749 nodalAtomicLambdaPower_(nullptr), 750 lambdaPowerFiltered_(nullptr), 751 atomLambdas_(nullptr), 752 atomThermostatForces_(nullptr), 753 atomMasses_(nullptr), 754 isFirstTimestep_(true), 755 nodalAtomicEnergy_(nullptr), 756 atomPredictedVelocities_(nullptr), 757 nodalAtomicPredictedEnergy_(nullptr), 758 firstHalfAtomForces_(nullptr) 759 { 760 // construct/obtain data corresponding to stage 3 of ATC_Method::initialize 761 nodalAtomicLambdaPower_ = thermostat->regulator_data(regulatorPrefix_+"NodalAtomicLambdaPower",1); 762 lambdaPowerFiltered_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaPowerFiltered",1); 763 } 764 765 //-------------------------------------------------------- 766 // constructor_transfers 767 // instantiates or obtains all dependency managed data 768 //-------------------------------------------------------- construct_transfers()769 void ThermostatGlcFs::construct_transfers() 770 { 771 lambdaSolver_->construct_transfers(); 772 InterscaleManager & interscaleManager(atc_->interscale_manager()); 773 774 // get atom velocity data from manager 775 atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY); 776 777 // construct lambda evaluated at atom locations 778 atomLambdas_ = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy"); 779 if (!atomLambdas_) { 780 DENS_MAN * lambda = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1); 781 atomLambdas_ = new FtaShapeFunctionProlongation(atc_, 782 lambda, 783 interscaleManager.per_atom_sparse_matrix("Interpolant")); 784 interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy"); 785 } 786 787 // get data from manager 788 atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS), 789 nodalAtomicEnergy_ = interscaleManager.dense_matrix("NodalAtomicEnergy"); 790 791 // thermostat forces based on lambda and the atomic velocities 792 atomThermostatForces_ = new AtomicThermostatForce(atc_,atomLambdas_); 793 interscaleManager.add_per_atom_quantity(atomThermostatForces_, 794 regulatorPrefix_+"AtomThermostatForce"); 795 796 // predicted temperature quantities: atom velocities, atom energies, and restricted atom energies 797 atomPredictedVelocities_ = new AtcAtomQuantity<double>(atc_,atc_->nsd()); 798 // MAKE THINGS WORK WITH ONLY ONE PREDICTED VELOCITY, CHECK IT EXISTS 799 interscaleManager.add_per_atom_quantity(atomPredictedVelocities_, 800 regulatorPrefix_+"AtomicPredictedVelocities"); 801 AtomicEnergyForTemperature * atomPredictedEnergyForTemperature = new TwiceKineticEnergy(atc_, 802 atomPredictedVelocities_); 803 interscaleManager.add_per_atom_quantity(atomPredictedEnergyForTemperature, 804 regulatorPrefix_+"AtomicPredictedTwiceKineticEnergy"); 805 nodalAtomicPredictedEnergy_ = new AtfShapeFunctionRestriction(atc_, 806 atomPredictedEnergyForTemperature, 807 interscaleManager.per_atom_sparse_matrix("Interpolant")); 808 interscaleManager.add_dense_matrix(nodalAtomicPredictedEnergy_, 809 regulatorPrefix_+"NodalAtomicPredictedEnergy"); 810 } 811 812 //-------------------------------------------------------- 813 // initialize 814 // initializes all method data 815 //-------------------------------------------------------- initialize()816 void ThermostatGlcFs::initialize() 817 { 818 RegulatorMethod::initialize(); 819 820 TimeFilterManager * timeFilterManager = atc_->time_filter_manager(); 821 if (!timeFilterManager->end_equilibrate()) { 822 // we should reset lambda and lambdaForce to zero in this case 823 // implies an initial condition of 0 for the filtered nodal lambda power 824 // initial conditions will always be needed when using time filtering 825 // however, the fractional step scheme must assume the instantaneous 826 // nodal lambda power is 0 initially because all quantities are in delta form 827 lambdaSolver_->initialize(); // ensures initial lambda force is zero 828 lambdaSolver_->set_lambda_to_value(0.); 829 *nodalAtomicLambdaPower_ = 0.; // energy change due to thermostats 830 *lambdaPowerFiltered_ = 0.; // filtered energy change due to thermostats 831 } 832 else { 833 lambdaSolver_->initialize(); 834 // we can grab lambda power variables using time integrator and atc transfer in cases for equilibration 835 } 836 837 // sets up time filter for cases where variables temporally filtered 838 if (timeFilterManager->need_reset()) { 839 // the form of this integrator implies no time filters that require history data can be used 840 timeFilter_->initialize(nodalAtomicLambdaPower_->quantity()); 841 } 842 843 atomThermostatForces_->quantity(); // initialize 844 atomThermostatForces_->fix_quantity(); 845 firstHalfAtomForces_ = atomThermostatForces_; // initialize 846 #ifdef OBSOLETE 847 compute_rhs_map(); 848 #endif 849 } 850 851 //-------------------------------------------------------- 852 // reset_atom_materials: 853 // resets the localized atom to material map 854 //-------------------------------------------------------- reset_atom_materials(const Array<int> & elementToMaterialMap,const MatrixDependencyManager<DenseMatrix,int> * atomElement)855 void ThermostatGlcFs::reset_atom_materials(const Array<int> & elementToMaterialMap, 856 const MatrixDependencyManager<DenseMatrix, int> * atomElement) 857 { 858 lambdaSolver_->reset_atom_materials(elementToMaterialMap, 859 atomElement); 860 } 861 862 //-------------------------------------------------------- 863 // apply_to_atoms: 864 // determines what if any contributions to the 865 // atomic moition is needed for 866 // consistency with the thermostat 867 // and computes the instantaneous induced power 868 //-------------------------------------------------------- apply_to_atoms(PerAtomQuantity<double> * atomicVelocity,const DENS_MAN * nodalAtomicEnergy,const DENS_MAT & lambdaForce,DENS_MAT & nodalAtomicLambdaPower,double dt)869 void ThermostatGlcFs::apply_to_atoms(PerAtomQuantity<double> * atomicVelocity, 870 const DENS_MAN * nodalAtomicEnergy, 871 const DENS_MAT & lambdaForce, 872 DENS_MAT & nodalAtomicLambdaPower, 873 double dt) 874 { 875 // compute initial contributions to lambda power 876 nodalAtomicLambdaPower = nodalAtomicEnergy->quantity(); 877 nodalAtomicLambdaPower *= -1.; 878 879 // apply lambda force to atoms 880 _velocityDelta_ = lambdaForce; 881 _velocityDelta_ /= atomMasses_->quantity(); 882 _velocityDelta_ *= dt; 883 (*atomicVelocity) += _velocityDelta_; 884 885 // finalize lambda power 886 nodalAtomicLambdaPower += nodalAtomicEnergy->quantity(); 887 } 888 889 //-------------------------------------------------------- 890 // full_prediction: 891 // flag to perform a full prediction calcalation 892 // for lambda rather than using the old value 893 //-------------------------------------------------------- full_prediction()894 bool ThermostatGlcFs::full_prediction() 895 { 896 if (isFirstTimestep_ || ((atc_->atom_to_element_map_type() == EULERIAN) 897 && (atc_->atom_to_element_map_frequency() > 1) 898 && (atc_->step() % atc_->atom_to_element_map_frequency() == 0 ))) { 899 return true; 900 } 901 return false; 902 } 903 904 //-------------------------------------------------------- 905 // apply_predictor: 906 // apply the thermostat to the atoms in the first step 907 // of the Verlet algorithm 908 //-------------------------------------------------------- apply_pre_predictor(double dt)909 void ThermostatGlcFs::apply_pre_predictor(double dt) 910 { 911 DENS_MAT & myLambdaPowerFiltered(lambdaPowerFiltered_->set_quantity()); 912 DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity()); 913 914 // update filtered power 915 timeFilter_->apply_pre_step1(myLambdaPowerFiltered,myNodalAtomicLambdaPower,dt); // equivalent update to measure power as change in energy due to thermostat 916 917 // apply lambda force to atoms and compute instantaneous lambda power for first half of time step 918 this->apply_to_atoms(atomVelocities_,nodalAtomicEnergy_, 919 firstHalfAtomForces_->quantity(), 920 myNodalAtomicLambdaPower,0.5*dt); 921 922 // update nodal variables for first half of time step 923 this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy1_,0.5*dt); 924 925 // start update of filtered lambda power 926 myNodalAtomicLambdaPower = 0.; // temporary power for first part of update 927 timeFilter_->apply_post_step1(myLambdaPowerFiltered,myNodalAtomicLambdaPower,dt); 928 } 929 930 //-------------------------------------------------------- 931 // apply_pre_corrector: 932 // apply the thermostat to the atoms in the first part 933 // of the corrector step of the Verlet algorithm 934 //-------------------------------------------------------- apply_pre_corrector(double dt)935 void ThermostatGlcFs::apply_pre_corrector(double dt) 936 { 937 (*atomPredictedVelocities_) = atomVelocities_->quantity(); 938 939 // do full prediction if we just redid the shape functions 940 if (full_prediction()) { 941 this->compute_lambda(dt); 942 943 atomThermostatForces_->unfix_quantity(); // allow update of atomic force applied by lambda 944 } 945 946 // apply lambda force to atoms and compute instantaneous lambda power to predict second half of time step 947 DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity()); 948 apply_to_atoms(atomPredictedVelocities_, 949 nodalAtomicPredictedEnergy_, 950 firstHalfAtomForces_->quantity(), 951 myNodalAtomicLambdaPower,0.5*dt); 952 953 if (full_prediction()) 954 atomThermostatForces_->fix_quantity(); 955 956 // update predicted nodal variables for second half of time step 957 this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy2_,0.5*dt); 958 // following manipulations performed this way for efficiency 959 deltaEnergy1_ += deltaEnergy2_; 960 atc_->apply_inverse_mass_matrix(deltaEnergy1_,TEMPERATURE); 961 temperature_ += deltaEnergy1_; 962 } 963 964 //-------------------------------------------------------- 965 // apply_post_corrector: 966 // apply the thermostat to the atoms in the second part 967 // of the corrector step of the Verlet algorithm 968 //-------------------------------------------------------- apply_post_corrector(double dt)969 void ThermostatGlcFs::apply_post_corrector(double dt) 970 { 971 // remove predicted power effects 972 DENS_MAT & myTemperature(temperature_.set_quantity()); 973 atc_->apply_inverse_mass_matrix(deltaEnergy2_,TEMPERATURE); 974 myTemperature -= deltaEnergy2_; 975 976 // set up equation and update lambda 977 this->compute_lambda(dt); 978 979 // apply lambda force to atoms and compute instantaneous lambda power for second half of time step 980 DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity()); 981 // allow computation of force applied by lambda using current velocities 982 atomThermostatForces_->unfix_quantity(); 983 atomThermostatForces_->quantity(); 984 atomThermostatForces_->fix_quantity(); 985 apply_to_atoms(atomVelocities_,nodalAtomicEnergy_, 986 atomThermostatForces_->quantity(), 987 myNodalAtomicLambdaPower,0.5*dt); 988 989 // finalize filtered lambda power by adding latest contribution 990 timeFilter_->apply_post_step2(lambdaPowerFiltered_->set_quantity(), 991 myNodalAtomicLambdaPower,dt); 992 993 // update nodal variables for second half of time step 994 this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy2_,0.5*dt); 995 atc_->apply_inverse_mass_matrix(deltaEnergy2_,TEMPERATURE); 996 myTemperature += deltaEnergy2_; 997 998 999 isFirstTimestep_ = false; 1000 } 1001 1002 //-------------------------------------------------------- 1003 // compute_lambda: 1004 // sets up and solves linear system for lambda, if the 1005 // bool is true it iterators to a non-linear solution 1006 //-------------------------------------------------------- compute_lambda(double dt,bool iterateSolution)1007 void ThermostatGlcFs::compute_lambda(double dt, 1008 bool iterateSolution) 1009 { 1010 // set up rhs for lambda equation 1011 this->set_thermostat_rhs(rhs_,0.5*dt); 1012 1013 // solve system 1014 lambdaSolver_->compute_lambda(rhs_,iterateSolution); 1015 #ifdef OBSOLETE 1016 // solve linear system for lambda guess 1017 DENS_MAT & lambda(lambda_->set_quantity()); 1018 solve_for_lambda(rhs_,lambda); 1019 1020 // iterate to solution 1021 if (iterateSolution) { 1022 iterate_lambda(rhs_); 1023 } 1024 #endif 1025 } 1026 1027 //-------------------------------------------------------- 1028 // compute_boundary_flux 1029 // default computation of boundary flux based on 1030 // finite 1031 //-------------------------------------------------------- compute_boundary_flux(FIELDS & fields)1032 void ThermostatGlcFs::compute_boundary_flux(FIELDS & fields) 1033 { 1034 1035 lambdaSolver_->compute_boundary_flux(fields); 1036 } 1037 1038 //-------------------------------------------------------- 1039 // output: 1040 // adds all relevant output to outputData 1041 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)1042 void ThermostatGlcFs::output(OUTPUT_LIST & outputData) 1043 { 1044 _lambdaPowerOutput_ = nodalAtomicLambdaPower_->quantity(); 1045 // approximate value for lambda power 1046 double dt = LammpsInterface::instance()->dt(); 1047 _lambdaPowerOutput_ *= (2./dt); 1048 DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity()); 1049 if ((atc_->lammps_interface())->rank_zero()) { 1050 outputData[regulatorPrefix_+"LambdaEnergy"] = λ 1051 outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_); 1052 } 1053 } 1054 1055 //-------------------------------------------------------- 1056 //-------------------------------------------------------- 1057 // Class ThermostatSolverFlux 1058 //-------------------------------------------------------- 1059 //-------------------------------------------------------- 1060 1061 //-------------------------------------------------------- 1062 // Constructor 1063 // Grab references to ATC and thermostat data 1064 //-------------------------------------------------------- ThermostatSolverFlux(AtomicRegulator * thermostat,int lambdaMaxIterations,const string & regulatorPrefix)1065 ThermostatSolverFlux::ThermostatSolverFlux(AtomicRegulator * thermostat, 1066 int lambdaMaxIterations, 1067 const string & regulatorPrefix) : 1068 ThermostatFsSolver(thermostat,lambdaMaxIterations,regulatorPrefix) 1069 { 1070 // do nothing 1071 } 1072 1073 //-------------------------------------------------------- 1074 // constructor_transfers 1075 // instantiates or obtains all dependency managed data 1076 //-------------------------------------------------------- construct_transfers()1077 void ThermostatSolverFlux::construct_transfers() 1078 { 1079 InterscaleManager & interscaleManager(atc_->interscale_manager()); 1080 1081 // set up node mappings 1082 create_node_maps(); 1083 1084 // set up data for linear solver 1085 shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_); 1086 interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_, 1087 regulatorPrefix_+"LambdaCouplingMatrixEnergy"); 1088 if (elementMask_) { 1089 lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_); 1090 interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_, 1091 regulatorPrefix_+"LambdaAtomMap"); 1092 } 1093 if (atomicRegulator_->use_localized_lambda()) { 1094 linearSolverType_ = AtomicRegulator::RSL_SOLVE; 1095 } 1096 else { 1097 linearSolverType_ = AtomicRegulator::CG_SOLVE; 1098 } 1099 1100 // base class transfers 1101 ThermostatFsSolver::construct_transfers(); 1102 1103 // add transfers for computation of extra RHS term accounting of O(lambda^2) 1104 // lambda squared followed by fractional step RHS contribution 1105 atomLambdas_ = (interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy")); 1106 if (!atomLambdas_) { 1107 atomLambdas_ = new FtaShapeFunctionProlongation(atc_, 1108 lambda_, 1109 interscaleManager.per_atom_sparse_matrix("Interpolant")); 1110 interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy"); 1111 } 1112 LambdaSquared * lambdaSquared = new LambdaSquared(atc_, 1113 interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS), 1114 weights_, 1115 atomLambdas_); 1116 interscaleManager.add_per_atom_quantity(lambdaSquared, 1117 regulatorPrefix_+"LambdaSquaredMapped"); 1118 rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_); 1119 interscaleManager.add_dense_matrix(rhsLambdaSquared_, 1120 regulatorPrefix_+"RhsLambdaSquared"); 1121 } 1122 1123 //-------------------------------------------------------- 1124 // construct_regulated_nodes: 1125 // constructs the set of nodes being regulated 1126 //-------------------------------------------------------- construct_regulated_nodes()1127 void ThermostatSolverFlux::construct_regulated_nodes() 1128 { 1129 InterscaleManager & interscaleManager(atc_->interscale_manager()); 1130 1131 // matrix requires all entries even if localized for correct lumping 1132 regulatedNodes_ = interscaleManager.set_int(regulatorPrefix_+"ThermostatRegulatedNodes"); 1133 if (!regulatedNodes_) { 1134 regulatedNodes_ = new RegulatedNodes(atc_); 1135 interscaleManager.add_set_int(regulatedNodes_, 1136 regulatorPrefix_+"ThermostatRegulatedNodes"); 1137 } 1138 1139 // if localized monitor nodes with applied fluxes 1140 if (atomicRegulator_->use_localized_lambda()) { 1141 if ((atomicRegulator_->coupling_mode() == Thermostat::FLUX) && (atomicRegulator_->boundary_integration_type() != NO_QUADRATURE)) { 1142 // include boundary nodes 1143 applicationNodes_ = new FluxBoundaryNodes(atc_); 1144 1145 boundaryNodes_ = new BoundaryNodes(atc_); 1146 interscaleManager.add_set_int(boundaryNodes_, 1147 regulatorPrefix_+"ThermostatBoundaryNodes"); 1148 } 1149 else { 1150 // fluxed nodes only 1151 applicationNodes_ = new FluxNodes(atc_); 1152 } 1153 interscaleManager.add_set_int(applicationNodes_, 1154 regulatorPrefix_+"ThermostatApplicationNodes"); 1155 } 1156 else { 1157 applicationNodes_ = regulatedNodes_; 1158 } 1159 1160 // special set of boundary elements for boundary flux quadrature 1161 if ((atomicRegulator_->boundary_integration_type() == FE_INTERPOLATION) 1162 && (atomicRegulator_->use_localized_lambda())) { 1163 elementMask_ = interscaleManager.dense_matrix_bool(regulatorPrefix_+"BoundaryElementMask"); 1164 if (!elementMask_) { 1165 elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_); 1166 interscaleManager.add_dense_matrix_bool(elementMask_, 1167 regulatorPrefix_+"BoundaryElementMask"); 1168 } 1169 } 1170 } 1171 1172 //-------------------------------------------------------- 1173 //-------------------------------------------------------- 1174 // Class ThermostatIntegratorFlux 1175 //-------------------------------------------------------- 1176 //-------------------------------------------------------- 1177 1178 //-------------------------------------------------------- 1179 // Constructor 1180 // Grab references to ATC and thermostat data 1181 //-------------------------------------------------------- ThermostatIntegratorFlux(AtomicRegulator * thermostat,int lambdaMaxIterations,const string & regulatorPrefix)1182 ThermostatIntegratorFlux::ThermostatIntegratorFlux(AtomicRegulator * thermostat, 1183 int lambdaMaxIterations, 1184 const string & regulatorPrefix) : 1185 ThermostatGlcFs(thermostat,lambdaMaxIterations,regulatorPrefix), 1186 heatSource_(atc_->atomic_source(TEMPERATURE)) 1187 { 1188 lambdaSolver_ = new ThermostatSolverFlux(thermostat, 1189 lambdaMaxIterations, 1190 regulatorPrefix); 1191 } 1192 1193 //-------------------------------------------------------- 1194 // add_to_temperature 1195 // add in contributions from lambda power and boundary 1196 // flux to the FE temperature 1197 //-------------------------------------------------------- add_to_energy(const DENS_MAT & nodalLambdaPower,DENS_MAT & deltaEnergy,double dt)1198 void ThermostatIntegratorFlux::add_to_energy(const DENS_MAT & nodalLambdaPower, 1199 DENS_MAT & deltaEnergy, 1200 double dt) 1201 { 1202 deltaEnergy.resize(nNodes_,1); 1203 const DENS_MAT & myBoundaryFlux(boundaryFlux_[TEMPERATURE].quantity()); 1204 for (int i = 0; i < nNodes_; i++) { 1205 deltaEnergy(i,0) = nodalLambdaPower(i,0) + dt*myBoundaryFlux(i,0); 1206 } 1207 } 1208 1209 //-------------------------------------------------------- 1210 // initialize 1211 // initializes all method data 1212 //-------------------------------------------------------- initialize()1213 void ThermostatIntegratorFlux::initialize() 1214 { 1215 ThermostatGlcFs::initialize(); 1216 1217 // timestep factor 1218 lambdaSolver_->set_timestep_factor(1.); 1219 } 1220 1221 //-------------------------------------------------------- 1222 // set_thermostat_rhs: 1223 // sets up the right-hand side including boundary 1224 // fluxes (coupling & prescribed), heat sources, and 1225 // fixed (uncoupled) nodes 1226 //-------------------------------------------------------- set_thermostat_rhs(DENS_MAT & rhs,double)1227 void ThermostatIntegratorFlux::set_thermostat_rhs(DENS_MAT & rhs, 1228 double /* dt */) 1229 { 1230 1231 // only tested with flux != 0 + ess bc = 0 1232 1233 // (a) for flux based : 1234 // form rhs : 2/3kB * W_I^-1 * \int N_I r dV 1235 // vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I 1236 // fluxes are set in ATC transfer 1237 const DENS_MAT & heatSource(heatSource_.quantity()); 1238 #if true 1239 const set<int> & applicationNodes((lambdaSolver_->application_nodes())->quantity()); 1240 rhs.resize(nNodes_,1); 1241 for (int i = 0; i < nNodes_; i++) { 1242 if (applicationNodes.find(i) != applicationNodes.end()) { 1243 rhs(i,0) = heatSource(i,0); 1244 } 1245 else { 1246 rhs(i,0) = 0.; 1247 } 1248 } 1249 #else 1250 rhs.resize(nNodes_,1); 1251 for (int i = 0; i < nNodes_; i++) { 1252 rhs(i,0) = heatSource(i,0); 1253 } 1254 #endif 1255 } 1256 1257 //-------------------------------------------------------- 1258 //-------------------------------------------------------- 1259 // Class ThermostatSolverFixed 1260 //-------------------------------------------------------- 1261 //-------------------------------------------------------- 1262 1263 //-------------------------------------------------------- 1264 // Constructor 1265 // Grab references to ATC and thermostat data 1266 //-------------------------------------------------------- ThermostatSolverFixed(AtomicRegulator * thermostat,int lambdaMaxIterations,const string & regulatorPrefix)1267 ThermostatSolverFixed::ThermostatSolverFixed(AtomicRegulator * thermostat, 1268 int lambdaMaxIterations, 1269 const string & regulatorPrefix) : 1270 ThermostatFsSolver(thermostat,lambdaMaxIterations,regulatorPrefix) 1271 { 1272 // do nothing 1273 } 1274 1275 //-------------------------------------------------------- 1276 // constructor_transfers 1277 // instantiates or obtains all dependency managed data 1278 //-------------------------------------------------------- construct_transfers()1279 void ThermostatSolverFixed::construct_transfers() 1280 { 1281 InterscaleManager & interscaleManager(atc_->interscale_manager()); 1282 1283 // set up node mappings 1284 create_node_maps(); 1285 1286 // determine if map is needed and set up if so 1287 if (this->use_local_shape_functions()) { 1288 lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_); 1289 interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_, 1290 regulatorPrefix_+"LambdaAtomMap"); 1291 shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_, 1292 lambdaAtomMap_, 1293 nodeToOverlapMap_); 1294 } 1295 else { 1296 shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_); 1297 } 1298 interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_, 1299 regulatorPrefix_+"LambdaCouplingMatrixEnergy"); 1300 linearSolverType_ = AtomicRegulator::CG_SOLVE; 1301 1302 // base class transfers, e.g. weights 1303 ThermostatFsSolver::construct_transfers(); 1304 1305 // add transfers for computation of extra RHS term accounting of O(lambda^2) 1306 // lambda squared followed by fractional step RHS contribution 1307 atomLambdas_ = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy"); 1308 if (!atomLambdas_) { 1309 atomLambdas_ = new FtaShapeFunctionProlongation(atc_, 1310 lambda_, 1311 interscaleManager.per_atom_sparse_matrix("Interpolant")); 1312 interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy"); 1313 } 1314 if (lambdaAtomMap_) { 1315 LambdaSquaredMapped * lambdaSquared = new LambdaSquaredMapped(atc_, 1316 lambdaAtomMap_, 1317 interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS), 1318 weights_, 1319 atomLambdas_); 1320 interscaleManager.add_per_atom_quantity(lambdaSquared, 1321 regulatorPrefix_+"LambdaSquared"); 1322 rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_); 1323 } 1324 else { 1325 LambdaSquared * lambdaSquared = new LambdaSquared(atc_, 1326 interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS), 1327 weights_, 1328 atomLambdas_); 1329 interscaleManager.add_per_atom_quantity(lambdaSquared, 1330 regulatorPrefix_+"LambdaSquaredMapped"); 1331 rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_); 1332 } 1333 interscaleManager.add_dense_matrix(rhsLambdaSquared_, 1334 regulatorPrefix_+"RhsLambdaSquared"); 1335 } 1336 1337 //-------------------------------------------------------- 1338 // construct_regulated_nodes: 1339 // constructs the set of nodes being regulated 1340 //-------------------------------------------------------- construct_regulated_nodes()1341 void ThermostatSolverFixed::construct_regulated_nodes() 1342 { 1343 InterscaleManager & interscaleManager(atc_->interscale_manager()); 1344 regulatedNodes_ = interscaleManager.set_int(regulatorPrefix_+"ThermostatRegulatedNodes"); 1345 1346 if (!regulatedNodes_) { 1347 if (!atomicRegulator_->use_localized_lambda()) { 1348 regulatedNodes_ = new RegulatedNodes(atc_); 1349 } 1350 else if (atomicRegulator_->coupling_mode() == AtomicRegulator::FLUX) { 1351 regulatedNodes_ = new FixedNodes(atc_); 1352 } 1353 else if (atomicRegulator_->coupling_mode() == AtomicRegulator::FIXED) { 1354 // include boundary nodes 1355 regulatedNodes_ = new FixedBoundaryNodes(atc_); 1356 } 1357 else { 1358 throw ATC_Error("ThermostatSolverFixed::construct_regulated_nodes - couldn't determine set of regulated nodes"); 1359 } 1360 1361 interscaleManager.add_set_int(regulatedNodes_, 1362 regulatorPrefix_+"ThermostatRegulatedNodes"); 1363 } 1364 1365 applicationNodes_ = regulatedNodes_; 1366 1367 // special set of boundary elements for defining regulated atoms 1368 if (atomicRegulator_->use_localized_lambda()) { 1369 elementMask_ = interscaleManager.dense_matrix_bool(regulatorPrefix_+"BoundaryElementMask"); 1370 if (!elementMask_) { 1371 elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_); 1372 interscaleManager.add_dense_matrix_bool(elementMask_, 1373 regulatorPrefix_+"BoundaryElementMask"); 1374 } 1375 } 1376 } 1377 1378 //-------------------------------------------------------- 1379 //-------------------------------------------------------- 1380 // Class ThermostatIntegratorFixed 1381 //-------------------------------------------------------- 1382 //-------------------------------------------------------- 1383 1384 //-------------------------------------------------------- 1385 // Constructor 1386 // Grab references to ATC and thermostat data 1387 //-------------------------------------------------------- ThermostatIntegratorFixed(AtomicRegulator * thermostat,int lambdaMaxIterations,const string & regulatorPrefix)1388 ThermostatIntegratorFixed::ThermostatIntegratorFixed(AtomicRegulator * thermostat, 1389 int lambdaMaxIterations, 1390 const string & regulatorPrefix) : 1391 ThermostatGlcFs(thermostat,lambdaMaxIterations,regulatorPrefix), 1392 atomThermostatForcesPredVel_(nullptr), 1393 filterCoefficient_(1.) 1394 { 1395 lambdaSolver_ = new ThermostatSolverFixed(thermostat, 1396 lambdaMaxIterations, 1397 regulatorPrefix); 1398 } 1399 1400 //-------------------------------------------------------- 1401 // construct_regulated_nodes: 1402 // constructs the set of nodes being regulated 1403 //-------------------------------------------------------- construct_transfers()1404 void ThermostatIntegratorFixed::construct_transfers() 1405 { 1406 ThermostatGlcFs::construct_transfers(); 1407 1408 InterscaleManager & interscaleManager(atc_->interscale_manager()); 1409 1410 // predicted forces for halving update 1411 atomThermostatForcesPredVel_ = new AtomicThermostatForce(atc_,atomLambdas_,atomPredictedVelocities_); 1412 interscaleManager.add_per_atom_quantity(atomThermostatForcesPredVel_, 1413 regulatorPrefix_+"AtomThermostatForcePredictedVelocity"); 1414 } 1415 1416 //-------------------------------------------------------- 1417 // initialize 1418 // initializes all method data 1419 //-------------------------------------------------------- initialize()1420 void ThermostatIntegratorFixed::initialize() 1421 { 1422 ThermostatGlcFs::initialize(); 1423 InterscaleManager & interscaleManager(atc_->interscale_manager()); 1424 1425 // set KE multiplier 1426 AtomicEnergyForTemperature * atomEnergyForTemperature = 1427 static_cast<AtomicEnergyForTemperature * >(interscaleManager.per_atom_quantity("AtomicEnergyForTemperature")); 1428 keMultiplier_ = atomEnergyForTemperature->kinetic_energy_multiplier(); 1429 1430 // reset data to zero 1431 deltaFeEnergy_.reset(nNodes_,1); 1432 deltaNodalAtomicEnergy_.reset(nNodes_,1); 1433 1434 // initialize filtered energy 1435 TimeFilterManager * timeFilterManager = atc_->time_filter_manager(); 1436 if (!timeFilterManager->end_equilibrate()) { 1437 nodalAtomicEnergyFiltered_ = nodalAtomicEnergy_->quantity(); 1438 } 1439 1440 // timestep factor 1441 lambdaSolver_->set_timestep_factor(0.5); 1442 } 1443 1444 //-------------------------------------------------------- 1445 // halve_force: 1446 // flag to halve the lambda force for improved 1447 // accuracy 1448 //-------------------------------------------------------- halve_force()1449 bool ThermostatIntegratorFixed::halve_force() 1450 { 1451 if (isFirstTimestep_ || ((atc_->atom_to_element_map_type() == EULERIAN) 1452 && (atc_->atom_to_element_map_frequency() > 1) 1453 && (atc_->step() % atc_->atom_to_element_map_frequency() == 1))) { 1454 return true; 1455 } 1456 return false; 1457 } 1458 1459 //-------------------------------------------------------- 1460 // initialize_delta_nodal_atomic_energy: 1461 // initializes storage for the variable tracking 1462 // the change in the nodal atomic energy 1463 // that has occurred over the past timestep 1464 //-------------------------------------------------------- initialize_delta_nodal_atomic_energy(double dt)1465 void ThermostatIntegratorFixed::initialize_delta_nodal_atomic_energy(double dt) 1466 { 1467 // initialize delta energy 1468 const DENS_MAT & myNodalAtomicEnergy(nodalAtomicEnergy_->quantity()); 1469 initialNodalAtomicEnergy_ = myNodalAtomicEnergy; 1470 initialNodalAtomicEnergy_ *= -1.; // initially stored as negative for efficiency 1471 timeFilter_->apply_pre_step1(nodalAtomicEnergyFiltered_.set_quantity(), 1472 myNodalAtomicEnergy,dt); 1473 } 1474 1475 //-------------------------------------------------------- 1476 // compute_delta_nodal_atomic_energy: 1477 // computes the change in the nodal atomic energy 1478 // that has occurred over the past timestep 1479 //-------------------------------------------------------- compute_delta_nodal_atomic_energy(double dt)1480 void ThermostatIntegratorFixed::compute_delta_nodal_atomic_energy(double dt) 1481 { 1482 // set delta energy based on predicted atomic velocities 1483 const DENS_MAT & myNodalAtomicEnergy(nodalAtomicEnergy_->quantity()); 1484 timeFilter_->apply_post_step1(nodalAtomicEnergyFiltered_.set_quantity(), 1485 myNodalAtomicEnergy,dt); 1486 deltaNodalAtomicEnergy_ = initialNodalAtomicEnergy_; 1487 deltaNodalAtomicEnergy_ += myNodalAtomicEnergy; 1488 } 1489 1490 //-------------------------------------------------------- 1491 // compute_lambda: 1492 // sets up and solves linear system for lambda 1493 //-------------------------------------------------------- compute_lambda(double dt,bool iterateSolution)1494 void ThermostatIntegratorFixed::compute_lambda(double dt, 1495 bool iterateSolution) 1496 { 1497 // compute predicted changes in nodal atomic energy 1498 compute_delta_nodal_atomic_energy(dt); 1499 1500 // change in finite element energy 1501 deltaFeEnergy_ = initialFeEnergy_; 1502 deltaFeEnergy_ += (mdMassMatrix_.quantity())*(temperature_.quantity()); 1503 1504 ThermostatGlcFs::compute_lambda(dt,iterateSolution); 1505 } 1506 1507 //-------------------------------------------------------- 1508 // apply_predictor: 1509 // apply the thermostat to the atoms in the first step 1510 // of the Verlet algorithm 1511 //-------------------------------------------------------- apply_pre_predictor(double dt)1512 void ThermostatIntegratorFixed::apply_pre_predictor(double dt) 1513 { 1514 // initialize values to be track change in finite element energy over the timestep 1515 initialize_delta_nodal_atomic_energy(dt); 1516 initialFeEnergy_ = -1.*((mdMassMatrix_.quantity())*(temperature_.quantity())); // initially stored as negative for efficiency 1517 1518 ThermostatGlcFs::apply_pre_predictor(dt); 1519 } 1520 1521 //-------------------------------------------------------- 1522 // apply_pre_corrector: 1523 // apply the thermostat to the atoms in the first part 1524 // of the corrector step of the Verlet algorithm 1525 //-------------------------------------------------------- apply_pre_corrector(double dt)1526 void ThermostatIntegratorFixed::apply_pre_corrector(double dt) 1527 { 1528 // do full prediction if we just redid the shape functions 1529 if (full_prediction()) { 1530 firstHalfAtomForces_ = atomThermostatForces_; // reset in case this time step needed special treatment 1531 _tempNodalAtomicEnergyFiltered_ = nodalAtomicEnergyFiltered_.quantity(); 1532 } 1533 1534 ThermostatGlcFs::apply_pre_corrector(dt); 1535 1536 if (full_prediction()) { 1537 // reset temporary variables 1538 nodalAtomicEnergyFiltered_ = _tempNodalAtomicEnergyFiltered_; 1539 } 1540 1541 if (halve_force()) { 1542 // save old velocities if we are doing halving calculation of lambda force 1543 // copy velocities over into temporary storage 1544 (*atomPredictedVelocities_) = atomVelocities_->quantity(); 1545 } 1546 } 1547 1548 //-------------------------------------------------------- 1549 // apply_post_corrector: 1550 // apply the thermostat to the atoms in the second part 1551 // of the corrector step of the Verlet algorithm 1552 //-------------------------------------------------------- apply_post_corrector(double dt)1553 void ThermostatIntegratorFixed::apply_post_corrector(double dt) 1554 { 1555 1556 bool halveForce = halve_force(); 1557 1558 ThermostatGlcFs::apply_post_corrector(dt); 1559 1560 // update filtered energy with lambda power 1561 DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity()); 1562 timeFilter_->apply_post_step2(nodalAtomicEnergyFiltered_.set_quantity(), 1563 myNodalAtomicLambdaPower,dt); 1564 1565 if (halveForce) { 1566 // Halve lambda force due to fixed temperature constraints 1567 // 1) makes up for poor initial condition 1568 // 2) accounts for possibly large value of lambda when atomic shape function values change 1569 // from eulerian mapping after more than 1 timestep 1570 // avoids unstable oscillations arising from 1571 // thermostat having to correct for error introduced in lambda changing the 1572 // shape function matrices 1573 lambdaSolver_->scale_lambda(0.5); 1574 firstHalfAtomForces_ = atomThermostatForcesPredVel_; 1575 atomThermostatForcesPredVel_->unfix_quantity(); 1576 } 1577 else { 1578 firstHalfAtomForces_ = atomThermostatForces_; 1579 } 1580 } 1581 1582 //-------------------------------------------------------- 1583 // add_to_temperature 1584 // add in contributions from lambda power and boundary 1585 // flux to the FE temperature 1586 //-------------------------------------------------------- add_to_energy(const DENS_MAT & nodalLambdaPower,DENS_MAT & deltaEnergy,double)1587 void ThermostatIntegratorFixed::add_to_energy(const DENS_MAT & nodalLambdaPower, 1588 DENS_MAT & deltaEnergy, 1589 double /* dt */) 1590 { 1591 deltaEnergy.resize(nNodes_,1); 1592 1593 SetDependencyManager<int> * myRegulatedNodes = 1594 (atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes"); 1595 const set<int> & regulatedNodes(myRegulatedNodes->quantity()); 1596 1597 for (int i = 0; i < nNodes_; i++) { 1598 if (regulatedNodes.find(i) != regulatedNodes.end()) { 1599 deltaEnergy(i,0) = 0.; 1600 } 1601 else { 1602 deltaEnergy(i,0) = nodalLambdaPower(i,0); 1603 } 1604 } 1605 } 1606 1607 //-------------------------------------------------------- 1608 // set_thermostat_rhs: 1609 // sets up the right-hand side including boundary 1610 // fluxes (coupling & prescribed), heat sources, and 1611 // fixed (uncoupled) nodes 1612 //-------------------------------------------------------- set_thermostat_rhs(DENS_MAT & rhs,double dt)1613 void ThermostatIntegratorFixed::set_thermostat_rhs(DENS_MAT & rhs, 1614 double dt) 1615 { 1616 // for essential bcs (fixed nodes) : 1617 // form rhs : (delThetaV - delTheta)/dt 1618 SetDependencyManager<int> * myRegulatedNodes = 1619 (atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes"); 1620 const set<int> & regulatedNodes(myRegulatedNodes->quantity()); 1621 double factor = (1./dt)/keMultiplier_; 1622 rhs.resize(nNodes_,1); 1623 1624 for (int i = 0; i < nNodes_; i++) { 1625 if (regulatedNodes.find(i) != regulatedNodes.end()) { 1626 rhs(i,0) = factor*(deltaNodalAtomicEnergy_(i,0) - deltaFeEnergy_(i,0)); 1627 } 1628 else { 1629 rhs(i,0) = 0.; 1630 } 1631 } 1632 } 1633 1634 //-------------------------------------------------------- 1635 //-------------------------------------------------------- 1636 // Class ThermostatIntegratorFluxFiltered 1637 //-------------------------------------------------------- 1638 //-------------------------------------------------------- 1639 1640 //-------------------------------------------------------- 1641 // Constructor 1642 // Grab references to ATC and thermostat data 1643 //-------------------------------------------------------- ThermostatIntegratorFluxFiltered(AtomicRegulator * thermostat,int lambdaMaxIterations,const string & regulatorPrefix)1644 ThermostatIntegratorFluxFiltered::ThermostatIntegratorFluxFiltered(AtomicRegulator * thermostat, 1645 int lambdaMaxIterations, 1646 const string & regulatorPrefix) : 1647 ThermostatIntegratorFlux(thermostat,lambdaMaxIterations,regulatorPrefix) 1648 { 1649 // do nothing 1650 } 1651 1652 //-------------------------------------------------------- 1653 // initialize 1654 // initializes all method data 1655 //-------------------------------------------------------- initialize()1656 void ThermostatIntegratorFluxFiltered::initialize() 1657 { 1658 ThermostatIntegratorFlux::initialize(); 1659 1660 TimeFilterManager * timeFilterManager = atc_->time_filter_manager(); 1661 if (!timeFilterManager->end_equilibrate()) { 1662 // always must start as zero because of filtering scheme 1663 heatSourceOld_.reset(nNodes_,1); 1664 instantHeatSource_.reset(nNodes_,1); 1665 timeStepSource_.reset(nNodes_,1); 1666 } 1667 } 1668 1669 //-------------------------------------------------------- 1670 // apply_post_corrector: 1671 // apply the thermostat to the atoms in the second part 1672 // of the corrector step of the Verlet algorithm 1673 //-------------------------------------------------------- apply_post_corrector(double dt)1674 void ThermostatIntegratorFluxFiltered::apply_post_corrector(double dt) 1675 { 1676 // compute lambda 1677 ThermostatIntegratorFlux::apply_post_corrector(dt); 1678 1679 // store data needed for filter inversion of heat flux for thermostat rhs 1680 instantHeatSource_ = rhs_; 1681 heatSourceOld_ = heatSource_.quantity(); 1682 } 1683 1684 //-------------------------------------------------------- 1685 // add_to_temperature 1686 // add in contributions from lambda power and boundary 1687 // flux to the FE temperature 1688 //-------------------------------------------------------- add_to_energy(const DENS_MAT & nodalLambdaPower,DENS_MAT & deltaEnergy,double dt)1689 void ThermostatIntegratorFluxFiltered::add_to_energy(const DENS_MAT & nodalLambdaPower, 1690 DENS_MAT & deltaEnergy, 1691 double dt) 1692 { 1693 deltaEnergy.reset(nNodes_,1); 1694 double coef = timeFilter_->unfiltered_coefficient_post_s1(2.*dt); 1695 const DENS_MAT & myBoundaryFlux(boundaryFlux_[TEMPERATURE].quantity()); 1696 for (int i = 0; i < nNodes_; i++) { 1697 deltaEnergy(i,0) = coef*nodalLambdaPower(i,0) + dt*myBoundaryFlux(i,0); 1698 } 1699 } 1700 1701 //-------------------------------------------------------- 1702 // set_thermostat_rhs: 1703 // sets up the right-hand side including boundary 1704 // fluxes (coupling & prescribed), heat sources, and 1705 // fixed (uncoupled) nodes 1706 //-------------------------------------------------------- set_thermostat_rhs(DENS_MAT & rhs,double dt)1707 void ThermostatIntegratorFluxFiltered::set_thermostat_rhs(DENS_MAT & rhs, 1708 double dt) 1709 { 1710 1711 // only tested with flux != 0 + ess bc = 0 1712 1713 // (a) for flux based : 1714 // form rhs : 2/3kB * W_I^-1 * \int N_I r dV 1715 // vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I 1716 // fluxes are set in ATC transfer 1717 1718 // invert heatSource_ to get unfiltered source 1719 // relevant coefficients from time filter 1720 1721 double coefF1 = timeFilter_->filtered_coefficient_pre_s1(2.*dt); 1722 double coefF2 = timeFilter_->filtered_coefficient_post_s1(2.*dt); 1723 double coefU1 = timeFilter_->unfiltered_coefficient_pre_s1(2.*dt); 1724 double coefU2 = timeFilter_->unfiltered_coefficient_post_s1(2.*dt); 1725 1726 const DENS_MAT & heatSource(heatSource_.quantity()); 1727 SetDependencyManager<int> * myApplicationNodes = 1728 (atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatApplicationNodes"); 1729 const set<int> & applicationNodes(myApplicationNodes->quantity()); 1730 rhs.resize(nNodes_,1); 1731 for (int i = 0; i < nNodes_; i++) { 1732 if (applicationNodes.find(i) != applicationNodes.end()) { 1733 rhs(i,0) = heatSource(i,0) - coefF1*coefF2*heatSourceOld_(i,0) - coefU1*coefF2*instantHeatSource_(i,0); 1734 rhs(i,0) /= coefU2; 1735 } 1736 else { 1737 rhs(i,0) = 0.; 1738 } 1739 } 1740 } 1741 1742 //-------------------------------------------------------- 1743 // output: 1744 // adds all relevant output to outputData 1745 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)1746 void ThermostatIntegratorFluxFiltered::output(OUTPUT_LIST & outputData) 1747 { 1748 _lambdaPowerOutput_ = lambdaPowerFiltered_->quantity(); 1749 // approximate value for lambda power 1750 double dt = LammpsInterface::instance()->dt(); 1751 _lambdaPowerOutput_ *= (2./dt); 1752 DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity()); 1753 if ((atc_->lammps_interface())->rank_zero()) { 1754 outputData[regulatorPrefix_+"Lambda"] = λ 1755 outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_); 1756 } 1757 } 1758 1759 //-------------------------------------------------------- 1760 //-------------------------------------------------------- 1761 // Class ThermostatIntegratorFixedFiltered 1762 //-------------------------------------------------------- 1763 //-------------------------------------------------------- 1764 1765 //-------------------------------------------------------- 1766 // Constructor 1767 // Grab references to ATC and thermostat data 1768 //-------------------------------------------------------- ThermostatIntegratorFixedFiltered(AtomicRegulator * thermostat,int lambdaMaxIterations,const string & regulatorPrefix)1769 ThermostatIntegratorFixedFiltered::ThermostatIntegratorFixedFiltered(AtomicRegulator * thermostat, 1770 int lambdaMaxIterations, 1771 const string & regulatorPrefix) : 1772 ThermostatIntegratorFixed(thermostat,lambdaMaxIterations,regulatorPrefix) 1773 { 1774 // do nothing 1775 } 1776 1777 //-------------------------------------------------------- 1778 // initialize_delta_nodal_atomic_energy: 1779 // initializes storage for the variable tracking 1780 // the change in the nodal atomic energy 1781 // that has occurred over the past timestep 1782 //-------------------------------------------------------- 1783 1784 initialize_delta_nodal_atomic_energy(double dt)1785 void ThermostatIntegratorFixedFiltered::initialize_delta_nodal_atomic_energy(double dt) 1786 { 1787 // initialize delta energy 1788 DENS_MAT & myNodalAtomicEnergyFiltered(nodalAtomicEnergyFiltered_.set_quantity()); 1789 initialNodalAtomicEnergy_ = myNodalAtomicEnergyFiltered; 1790 initialNodalAtomicEnergy_ *= -1.; // initially stored as negative for efficiency 1791 timeFilter_->apply_pre_step1(myNodalAtomicEnergyFiltered, 1792 nodalAtomicEnergy_->quantity(),dt); 1793 } 1794 1795 //-------------------------------------------------------- 1796 // compute_delta_nodal_atomic_energy: 1797 // computes the change in the nodal atomic energy 1798 // that has occurred over the past timestep 1799 //-------------------------------------------------------- compute_delta_nodal_atomic_energy(double dt)1800 void ThermostatIntegratorFixedFiltered::compute_delta_nodal_atomic_energy(double dt) 1801 { 1802 // set delta energy based on predicted atomic velocities 1803 DENS_MAT & myNodalAtomicEnergyFiltered(nodalAtomicEnergyFiltered_.set_quantity()); 1804 timeFilter_->apply_post_step1(myNodalAtomicEnergyFiltered, 1805 nodalAtomicEnergy_->quantity(),dt); 1806 deltaNodalAtomicEnergy_ = initialNodalAtomicEnergy_; 1807 deltaNodalAtomicEnergy_ += myNodalAtomicEnergyFiltered; 1808 } 1809 1810 //-------------------------------------------------------- 1811 // add_to_temperature 1812 // add in contributions from lambda power and boundary 1813 // flux to the FE temperature 1814 //-------------------------------------------------------- add_to_energy(const DENS_MAT & nodalLambdaPower,DENS_MAT & deltaEnergy,double dt)1815 void ThermostatIntegratorFixedFiltered::add_to_energy(const DENS_MAT & nodalLambdaPower, 1816 DENS_MAT & deltaEnergy, 1817 double dt) 1818 { 1819 deltaEnergy.resize(nNodes_,1); 1820 SetDependencyManager<int> * myRegulatedNodes = 1821 (atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes"); 1822 const set<int> & regulatedNodes(myRegulatedNodes->quantity()); 1823 double coef = timeFilter_->unfiltered_coefficient_post_s1(2.*dt); 1824 for (int i = 0; i < nNodes_; i++) { 1825 if (regulatedNodes.find(i) != regulatedNodes.end()) { 1826 deltaEnergy(i,0) = 0.; 1827 } 1828 else { 1829 deltaEnergy(i,0) = coef*nodalLambdaPower(i,0); 1830 } 1831 } 1832 } 1833 //-------------------------------------------------------- 1834 // set_thermostat_rhs: 1835 // sets up the right-hand side for fixed 1836 // (coupling & prescribed) temperature values 1837 //-------------------------------------------------------- set_thermostat_rhs(DENS_MAT & rhs,double dt)1838 void ThermostatIntegratorFixedFiltered::set_thermostat_rhs(DENS_MAT & rhs, 1839 double dt) 1840 { 1841 // (b) for essential bcs (fixed nodes): 1842 // form rhs : (delThetaV - delTheta)/dt 1843 SetDependencyManager<int> * myRegulatedNodes = 1844 (atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes"); 1845 const set<int> & regulatedNodes(myRegulatedNodes->quantity()); 1846 double factor = (1./dt)/keMultiplier_; 1847 factor /= timeFilter_->unfiltered_coefficient_post_s1(2.*dt); 1848 rhs.resize(nNodes_,1); 1849 1850 for (int i = 0; i < nNodes_; i++) { 1851 if (regulatedNodes.find(i) != regulatedNodes.end()) { 1852 rhs(i,0) = factor*(deltaNodalAtomicEnergy_(i,0) - deltaFeEnergy_(i,0)); 1853 } 1854 else { 1855 rhs(i,0) = 0.; 1856 } 1857 } 1858 } 1859 1860 //-------------------------------------------------------- 1861 // output: 1862 // adds all relevant output to outputData 1863 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)1864 void ThermostatIntegratorFixedFiltered::output(OUTPUT_LIST & outputData) 1865 { 1866 _lambdaPowerOutput_ = lambdaPowerFiltered_->quantity(); 1867 // approximate value for lambda power 1868 double dt = LammpsInterface::instance()->dt(); 1869 _lambdaPowerOutput_ *= (2./dt); 1870 DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity()); 1871 if ((atc_->lammps_interface())->rank_zero()) { 1872 outputData[regulatorPrefix_+"Lambda"] = λ 1873 outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_); 1874 } 1875 } 1876 1877 //-------------------------------------------------------- 1878 //-------------------------------------------------------- 1879 // Class ThermostatFluxFixed 1880 //-------------------------------------------------------- 1881 //-------------------------------------------------------- 1882 1883 //-------------------------------------------------------- 1884 // Constructor 1885 //-------------------------------------------------------- ThermostatFluxFixed(AtomicRegulator * thermostat,int lambdaMaxIterations,bool constructThermostats)1886 ThermostatFluxFixed::ThermostatFluxFixed(AtomicRegulator * thermostat, 1887 int lambdaMaxIterations, 1888 bool constructThermostats) : 1889 RegulatorMethod(thermostat), 1890 thermostatFlux_(nullptr), 1891 thermostatFixed_(nullptr), 1892 thermostatBcs_(nullptr) 1893 { 1894 if (constructThermostats) { 1895 thermostatFlux_ = new ThermostatIntegratorFlux(thermostat,lambdaMaxIterations,regulatorPrefix_+"Flux"); 1896 thermostatFixed_ = new ThermostatIntegratorFixed(thermostat,lambdaMaxIterations,regulatorPrefix_+"Fixed"); 1897 1898 // need to choose BC type based on coupling mode 1899 if (thermostat->coupling_mode() == AtomicRegulator::FLUX) { 1900 thermostatBcs_ = thermostatFlux_; 1901 } 1902 else if (thermostat->coupling_mode() == AtomicRegulator::FIXED) { 1903 thermostatBcs_ = thermostatFixed_; 1904 } 1905 else { 1906 throw ATC_Error("ThermostatFluxFixed:create_thermostats - invalid thermostat type provided"); 1907 } 1908 } 1909 } 1910 1911 //-------------------------------------------------------- 1912 // Destructor 1913 //-------------------------------------------------------- ~ThermostatFluxFixed()1914 ThermostatFluxFixed::~ThermostatFluxFixed() 1915 { 1916 if (thermostatFlux_) delete thermostatFlux_; 1917 if (thermostatFixed_) delete thermostatFixed_; 1918 } 1919 1920 //-------------------------------------------------------- 1921 // constructor_transfers 1922 // instantiates or obtains all dependency managed data 1923 //-------------------------------------------------------- construct_transfers()1924 void ThermostatFluxFixed::construct_transfers() 1925 { 1926 thermostatFlux_->construct_transfers(); 1927 thermostatFixed_->construct_transfers(); 1928 } 1929 1930 //-------------------------------------------------------- 1931 // initialize 1932 // initializes all method data 1933 //-------------------------------------------------------- initialize()1934 void ThermostatFluxFixed::initialize() 1935 { 1936 thermostatFixed_->initialize(); 1937 thermostatFlux_->initialize(); 1938 } 1939 1940 //-------------------------------------------------------- 1941 // apply_predictor: 1942 // apply the thermostat to the atoms in the first step 1943 // of the Verlet algorithm 1944 //-------------------------------------------------------- apply_pre_predictor(double dt)1945 void ThermostatFluxFixed::apply_pre_predictor(double dt) 1946 { 1947 thermostatFixed_->apply_pre_predictor(dt); 1948 thermostatFlux_->apply_pre_predictor(dt); 1949 } 1950 1951 //-------------------------------------------------------- 1952 // apply_pre_corrector: 1953 // apply the thermostat to the atoms in the first part 1954 // of the corrector step of the Verlet algorithm 1955 //-------------------------------------------------------- apply_pre_corrector(double dt)1956 void ThermostatFluxFixed::apply_pre_corrector(double dt) 1957 { 1958 thermostatFlux_->apply_pre_corrector(dt); 1959 if (thermostatFixed_->full_prediction()) { 1960 atc_->set_fixed_nodes(); 1961 } 1962 thermostatFixed_->apply_pre_corrector(dt); 1963 } 1964 1965 //-------------------------------------------------------- 1966 // apply_post_corrector: 1967 // apply the thermostat to the atoms in the second part 1968 // of the corrector step of the Verlet algorithm 1969 //-------------------------------------------------------- apply_post_corrector(double dt)1970 void ThermostatFluxFixed::apply_post_corrector(double dt) 1971 { 1972 thermostatFlux_->apply_post_corrector(dt); 1973 atc_->set_fixed_nodes(); 1974 thermostatFixed_->apply_post_corrector(dt); 1975 } 1976 1977 //-------------------------------------------------------- 1978 // output: 1979 // adds all relevant output to outputData 1980 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)1981 void ThermostatFluxFixed::output(OUTPUT_LIST & outputData) 1982 { 1983 thermostatFlux_->output(outputData); 1984 thermostatFixed_->output(outputData); 1985 } 1986 1987 //-------------------------------------------------------- 1988 //-------------------------------------------------------- 1989 // Class ThermostatFluxFixedFiltered 1990 //-------------------------------------------------------- 1991 //-------------------------------------------------------- 1992 1993 //-------------------------------------------------------- 1994 // Constructor 1995 //-------------------------------------------------------- ThermostatFluxFixedFiltered(AtomicRegulator * thermostat,int lambdaMaxIterations)1996 ThermostatFluxFixedFiltered::ThermostatFluxFixedFiltered(AtomicRegulator * thermostat, 1997 int lambdaMaxIterations) : 1998 ThermostatFluxFixed(thermostat,lambdaMaxIterations,false) 1999 { 2000 thermostatFlux_ = new ThermostatIntegratorFluxFiltered(thermostat,lambdaMaxIterations,regulatorPrefix_+"Flux"); 2001 thermostatFixed_ = new ThermostatIntegratorFixedFiltered(thermostat,lambdaMaxIterations,regulatorPrefix_+"Fixed"); 2002 2003 // need to choose BC type based on coupling mode 2004 if (thermostat->coupling_mode() == AtomicRegulator::FLUX) { 2005 thermostatBcs_ = thermostatFlux_; 2006 } 2007 else if (thermostat->coupling_mode() == AtomicRegulator::FIXED) { 2008 thermostatBcs_ = thermostatFixed_; 2009 } 2010 else { 2011 throw ATC_Error("ThermostatFluxFixed:create_thermostats - invalid thermostat type provided"); 2012 } 2013 } 2014 //-------------------------------------------------------- 2015 // Class ThermostatGlc 2016 //-------------------------------------------------------- 2017 //-------------------------------------------------------- 2018 2019 //-------------------------------------------------------- 2020 // Constructor 2021 //-------------------------------------------------------- ThermostatGlc(AtomicRegulator * thermostat)2022 ThermostatGlc::ThermostatGlc(AtomicRegulator * thermostat) : 2023 ThermostatShapeFunction(thermostat), 2024 timeFilter_(atomicRegulator_->time_filter()), 2025 lambdaPowerFiltered_(nullptr), 2026 atomThermostatForces_(nullptr), 2027 prescribedDataMgr_(atc_->prescribed_data_manager()), 2028 atomMasses_(nullptr) 2029 { 2030 // consistent with stage 3 of ATC_Method::initialize 2031 lambdaPowerFiltered_= atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaPowerFiltered",1); 2032 } 2033 2034 //-------------------------------------------------------- 2035 // constructor_transfers 2036 // instantiates or obtains all dependency managed data 2037 //-------------------------------------------------------- construct_transfers()2038 void ThermostatGlc::construct_transfers() 2039 { 2040 ThermostatShapeFunction::construct_transfers(); 2041 InterscaleManager & interscaleManager(atc_->interscale_manager()); 2042 2043 // get data from manager 2044 atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS); 2045 2046 // thermostat forces based on lambda and the atomic velocities 2047 AtomicThermostatForce * atomThermostatForces = new AtomicThermostatForce(atc_); 2048 interscaleManager.add_per_atom_quantity(atomThermostatForces, 2049 regulatorPrefix_+"AtomThermostatForce"); 2050 atomThermostatForces_ = atomThermostatForces; 2051 } 2052 2053 //-------------------------------------------------------- 2054 // apply_to_atoms: 2055 // determines what if any contributions to the 2056 // atomic moition is needed for 2057 // consistency with the thermostat 2058 //-------------------------------------------------------- apply_to_atoms(PerAtomQuantity<double> * atomVelocities,const DENS_MAT & lambdaForce,double dt)2059 void ThermostatGlc::apply_to_atoms(PerAtomQuantity<double> * atomVelocities, 2060 const DENS_MAT & lambdaForce, 2061 double dt) 2062 { 2063 _velocityDelta_ = lambdaForce; 2064 _velocityDelta_ /= atomMasses_->quantity(); 2065 _velocityDelta_ *= dt; 2066 (*atomVelocities) += _velocityDelta_; 2067 } 2068 2069 //-------------------------------------------------------- 2070 //-------------------------------------------------------- 2071 // Class ThermostatPowerVerlet 2072 //-------------------------------------------------------- 2073 //-------------------------------------------------------- 2074 2075 //-------------------------------------------------------- 2076 // Constructor 2077 // Grab references to ATC and thermostat data 2078 //-------------------------------------------------------- ThermostatPowerVerlet(AtomicRegulator * thermostat)2079 ThermostatPowerVerlet::ThermostatPowerVerlet(AtomicRegulator * thermostat) : 2080 ThermostatGlc(thermostat), 2081 nodalTemperatureRoc_(atc_->field_roc(TEMPERATURE)), 2082 heatSource_(atc_->atomic_source(TEMPERATURE)), 2083 nodalAtomicPower_(nullptr), 2084 nodalAtomicLambdaPower_(nullptr) 2085 { 2086 // do nothing 2087 } 2088 2089 //-------------------------------------------------------- 2090 // constructor_transfers 2091 // instantiates or obtains all dependency managed data 2092 //-------------------------------------------------------- construct_transfers()2093 void ThermostatPowerVerlet::construct_transfers() 2094 { 2095 InterscaleManager & interscaleManager(atc_->interscale_manager()); 2096 2097 // set up node mappings 2098 create_node_maps(); 2099 2100 // determine if mapping is needed and set up if so 2101 if (atomicRegulator_->use_localized_lambda()) { 2102 lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_); 2103 interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_, 2104 regulatorPrefix_+"LambdaAtomMap"); 2105 } 2106 2107 // set up linear solver 2108 if (atomicRegulator_->use_lumped_lambda_solve()) { 2109 shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_); 2110 linearSolverType_ = AtomicRegulator::RSL_SOLVE; 2111 } 2112 else { 2113 if (lambdaAtomMap_) { 2114 shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_, 2115 lambdaAtomMap_, 2116 nodeToOverlapMap_); 2117 } 2118 else { 2119 shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_); 2120 } 2121 linearSolverType_ = AtomicRegulator::CG_SOLVE; 2122 } 2123 interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_, 2124 regulatorPrefix_+"LambdaCouplingMatrixEnergy"); 2125 2126 // base class transfers, e.g. weights 2127 ThermostatGlc::construct_transfers(); 2128 2129 // get managed data 2130 nodalAtomicPower_ = interscaleManager.dense_matrix("NodalAtomicPower"); 2131 2132 // power induced by lambda 2133 DotTwiceKineticEnergy * atomicLambdaPower = 2134 new DotTwiceKineticEnergy(atc_,atomThermostatForces_); 2135 interscaleManager.add_per_atom_quantity(atomicLambdaPower, 2136 regulatorPrefix_+"AtomicLambdaPower"); 2137 2138 // restriction to nodes of power induced by lambda 2139 nodalAtomicLambdaPower_ = new AtfShapeFunctionRestriction(atc_, 2140 atomicLambdaPower, 2141 interscaleManager.per_atom_sparse_matrix("Interpolant")); 2142 interscaleManager.add_dense_matrix(nodalAtomicLambdaPower_, 2143 regulatorPrefix_+"NodalAtomicLambdaPower"); 2144 } 2145 2146 //-------------------------------------------------------- 2147 // initialize 2148 // initializes all method data 2149 //-------------------------------------------------------- initialize()2150 void ThermostatPowerVerlet::initialize() 2151 { 2152 ThermostatGlc::initialize(); 2153 2154 // sets up time filter for cases where variables temporally filtered 2155 TimeFilterManager * timeFilterManager = atc_->time_filter_manager(); 2156 if (!timeFilterManager->end_equilibrate()) { 2157 _nodalAtomicLambdaPowerOut_ = 0.; 2158 *lambdaPowerFiltered_ = 0.; 2159 timeFilter_->initialize(lambdaPowerFiltered_->quantity()); 2160 } 2161 } 2162 2163 //-------------------------------------------------------- 2164 // apply_predictor: 2165 // apply the thermostat to the atoms in the first step 2166 // of the Verlet algorithm 2167 //-------------------------------------------------------- apply_pre_predictor(double dt)2168 void ThermostatPowerVerlet::apply_pre_predictor(double dt) 2169 { 2170 atomThermostatForces_->unfix_quantity(); 2171 compute_thermostat(0.5*dt); 2172 2173 // apply lambda force to atoms 2174 const DENS_MAT & thermostatForces(atomThermostatForces_->quantity()); 2175 atomThermostatForces_->fix_quantity(); 2176 apply_to_atoms(atomVelocities_,thermostatForces,0.5*dt); 2177 } 2178 2179 //-------------------------------------------------------- 2180 // apply_pre_corrector: 2181 // apply the thermostat to the atoms in the first part 2182 // of the corrector step of the Verlet algorithm 2183 //-------------------------------------------------------- apply_pre_corrector(double dt)2184 void ThermostatPowerVerlet::apply_pre_corrector(double dt) 2185 { 2186 atomThermostatForces_->unfix_quantity(); 2187 compute_thermostat(0.5*dt); 2188 2189 // apply lambda force to atoms 2190 const DENS_MAT & thermostatForces(atomThermostatForces_->quantity()); 2191 atomThermostatForces_->fix_quantity(); 2192 apply_to_atoms(atomVelocities_,thermostatForces,0.5*dt); 2193 } 2194 2195 //-------------------------------------------------------- 2196 // add_to_rhs: 2197 // determines what if any contributions to the 2198 // finite element equations are needed for 2199 // consistency with the thermostat 2200 //-------------------------------------------------------- add_to_rhs(FIELDS & rhs)2201 void ThermostatPowerVerlet::add_to_rhs(FIELDS & rhs) 2202 { 2203 rhs[TEMPERATURE] += nodalAtomicLambdaPower_->quantity() + boundaryFlux_[TEMPERATURE].quantity(); 2204 } 2205 2206 //-------------------------------------------------------- 2207 // set_thermostat_rhs: 2208 // sets up the right-hand side including boundary 2209 // fluxes (coupling & prescribed), heat sources, and 2210 // fixed (uncoupled) nodes 2211 //-------------------------------------------------------- set_thermostat_rhs(DENS_MAT & rhs_nodes)2212 void ThermostatPowerVerlet::set_thermostat_rhs(DENS_MAT & rhs_nodes) 2213 { 2214 // (a) for flux based : 2215 // form rhs : \int N_I r dV 2216 // vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I 2217 // fluxes are set in ATC transfer 2218 rhs_nodes = heatSource_.quantity(); 2219 2220 // (b) for ess. bcs 2221 // form rhs : {sum_a (2 * N_Ia * v_ia * f_ia) - (dtheta/dt)_I} 2222 2223 // replace rhs for prescribed nodes 2224 const DENS_MAT & myNodalAtomicPower(nodalAtomicPower_->quantity()); 2225 const DIAG_MAT & myMdMassMatrix(mdMassMatrix_.quantity()); 2226 const DENS_MAT & myNodalTemperatureRoc(nodalTemperatureRoc_.quantity()); 2227 for (int i = 0; i < nNodes_; i++) { 2228 if (prescribedDataMgr_->is_fixed(i,TEMPERATURE,0)) { 2229 rhs_nodes(i,0) = 0.5*(myNodalAtomicPower(i,0) - myMdMassMatrix(i,i)*myNodalTemperatureRoc(i,0)); 2230 } 2231 } 2232 } 2233 2234 //-------------------------------------------------------- 2235 // compute_thermostat: 2236 // sets up and solves the thermostat equations since 2237 // they are the same at different parts of the time 2238 // step 2239 //-------------------------------------------------------- compute_thermostat(double dt)2240 void ThermostatPowerVerlet::compute_thermostat(double dt) 2241 { 2242 // set up rhs 2243 set_thermostat_rhs(_rhs_); 2244 2245 // solve linear system for lambda 2246 DENS_MAT & myLambda(lambda_->set_quantity()); 2247 solve_for_lambda(_rhs_,myLambda); 2248 2249 nodalAtomicLambdaPower_->unfix_quantity(); // enable computation of force applied by lambda 2250 timeFilter_->apply_pre_step1(lambdaPowerFiltered_->set_quantity(), 2251 nodalAtomicLambdaPower_->quantity(),dt); 2252 nodalAtomicLambdaPower_->fix_quantity(); 2253 } 2254 2255 //-------------------------------------------------------- 2256 // output: 2257 // adds all relevant output to outputData 2258 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)2259 void ThermostatPowerVerlet::output(OUTPUT_LIST & outputData) 2260 { 2261 _nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); 2262 DENS_MAT & lambda(lambda_->set_quantity()); 2263 if ((atc_->lammps_interface())->rank_zero()) { 2264 outputData["lambda"] = λ 2265 outputData["nodalLambdaPower"] = &(_nodalAtomicLambdaPowerOut_); 2266 } 2267 } 2268 2269 //-------------------------------------------------------- 2270 // finish: 2271 // final tasks after a run 2272 //-------------------------------------------------------- finish()2273 void ThermostatPowerVerlet::finish() 2274 { 2275 _nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); 2276 } 2277 2278 //-------------------------------------------------------- 2279 //-------------------------------------------------------- 2280 // Class ThermostatHooverVerlet 2281 //-------------------------------------------------------- 2282 //-------------------------------------------------------- 2283 2284 //-------------------------------------------------------- 2285 // Constructor 2286 // Grab references to ATC and thermostat data 2287 //-------------------------------------------------------- ThermostatHooverVerlet(AtomicRegulator * thermostat)2288 ThermostatHooverVerlet::ThermostatHooverVerlet(AtomicRegulator * thermostat) : 2289 ThermostatPowerVerlet(thermostat), 2290 lambdaHoover_(nullptr), 2291 nodalAtomicHooverLambdaPower_(nullptr) 2292 { 2293 // set up data consistent with stage 3 of ATC_Method::initialize 2294 lambdaHoover_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaHoover",1); 2295 2296 } 2297 2298 //-------------------------------------------------------- 2299 // constructor_transfers 2300 // instantiates or obtains all dependency managed data 2301 //-------------------------------------------------------- construct_transfers()2302 void ThermostatHooverVerlet::construct_transfers() 2303 { 2304 ThermostatPowerVerlet::construct_transfers(); 2305 InterscaleManager & interscaleManager(atc_->interscale_manager()); 2306 2307 2308 FtaShapeFunctionProlongation * atomHooverLambdas = new FtaShapeFunctionProlongation(atc_, 2309 lambdaHoover_, 2310 interscaleManager.per_atom_sparse_matrix("Interpolant")); 2311 interscaleManager.add_per_atom_quantity(atomHooverLambdas, 2312 regulatorPrefix_+"AtomHooverLambda"); 2313 AtomicThermostatForce * atomHooverThermostatForces = new AtomicThermostatForce(atc_,atomHooverLambdas); 2314 interscaleManager.add_per_atom_quantity(atomHooverThermostatForces, 2315 regulatorPrefix_+"AtomHooverThermostatForce"); 2316 SummedAtomicQuantity<double> * atomTotalThermostatForces = 2317 new SummedAtomicQuantity<double>(atc_,atomThermostatForces_,atomHooverThermostatForces); 2318 interscaleManager.add_per_atom_quantity(atomTotalThermostatForces, 2319 regulatorPrefix_+"AtomTotalThermostatForce"); 2320 atomThermostatForces_ = atomTotalThermostatForces; 2321 2322 // transfers dependent on time integration method 2323 DotTwiceKineticEnergy * atomicHooverLambdaPower = 2324 new DotTwiceKineticEnergy(atc_,atomHooverThermostatForces); 2325 interscaleManager.add_per_atom_quantity(atomicHooverLambdaPower, 2326 regulatorPrefix_+"AtomicHooverLambdaPower"); 2327 2328 nodalAtomicHooverLambdaPower_ = new AtfShapeFunctionRestriction(atc_, 2329 atomicHooverLambdaPower, 2330 interscaleManager.per_atom_sparse_matrix("Interpolant")); 2331 interscaleManager.add_dense_matrix(nodalAtomicHooverLambdaPower_, 2332 regulatorPrefix_+"NodalAtomicHooverLambdaPower"); 2333 } 2334 2335 //-------------------------------------------------------- 2336 // add_to_rhs: 2337 // determines what if any contributions to the 2338 // finite element equations are needed for 2339 // consistency with the thermostat 2340 //-------------------------------------------------------- add_to_rhs(FIELDS & rhs)2341 void ThermostatHooverVerlet::add_to_rhs(FIELDS & rhs) 2342 { 2343 rhs[TEMPERATURE] += _nodalAtomicLambdaPowerOut_; 2344 } 2345 2346 //-------------------------------------------------------- 2347 // compute_thermostat: 2348 // sets up and solves the thermostat equations since 2349 // they are the same at different parts of the time 2350 // step 2351 //-------------------------------------------------------- compute_thermostat(double dt)2352 void ThermostatHooverVerlet::compute_thermostat(double dt) 2353 { 2354 // apply prescribed/extrinsic sources and fixed nodes 2355 ThermostatPowerVerlet::compute_thermostat(0.5*dt); 2356 _nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); // save power from lambda in power-based thermostat 2357 2358 // set up Hoover rhs 2359 set_hoover_rhs(_rhs_); 2360 2361 // solve linear system for lambda 2362 DENS_MAT & myLambda(lambdaHoover_->set_quantity()); 2363 solve_for_lambda(_rhs_,myLambda); 2364 2365 // compute force applied by lambda 2366 // compute nodal atomic power from Hoover coupling 2367 // only add in contribution to uncoupled nodes 2368 if (atomicRegulator_->use_localized_lambda()) 2369 add_to_lambda_power(atomThermostatForces_->quantity(),0.5*dt); 2370 } 2371 2372 //-------------------------------------------------------- 2373 // set_hoover_rhs: 2374 // sets up the right-hand side for fixed value, 2375 // i.e. Hoover coupling 2376 //-------------------------------------------------------- set_hoover_rhs(DENS_MAT & rhs)2377 void ThermostatHooverVerlet::set_hoover_rhs(DENS_MAT & rhs) 2378 { 2379 // form rhs : sum_a ( N_Ia * v_ia * f_ia) - 0.5*M_MD*(dtheta/dt)_I 2380 rhs = nodalAtomicPower_->quantity(); 2381 rhs -= mdMassMatrix_.quantity()*nodalTemperatureRoc_.quantity(); 2382 rhs /= 2.; 2383 } 2384 2385 //-------------------------------------------------------- 2386 // add_to_nodal_lambda_power: 2387 // determines the power exerted by the Hoover 2388 // thermostat at each FE node 2389 //-------------------------------------------------------- add_to_lambda_power(const DENS_MAT &,double dt)2390 void ThermostatHooverVerlet::add_to_lambda_power(const DENS_MAT & /* myLambdaForce */, 2391 double dt) 2392 { 2393 _myNodalLambdaPower_ = nodalAtomicHooverLambdaPower_->quantity(); 2394 const INT_ARRAY & nodeToOverlapMap(nodeToOverlapMap_->quantity()); 2395 for (int i = 0; i < nNodes_; ++i) { 2396 if (nodeToOverlapMap(i,0)==-1) 2397 _nodalAtomicLambdaPowerOut_(i,0) += _myNodalLambdaPower_(i,0); 2398 else 2399 _myNodalLambdaPower_(i,0) = 0.; 2400 } 2401 timeFilter_->apply_post_step1(lambdaPowerFiltered_->set_quantity(),_myNodalLambdaPower_,dt); 2402 } 2403 2404 //-------------------------------------------------------- 2405 //-------------------------------------------------------- 2406 // Class ThermostatPowerVerletFiltered 2407 //-------------------------------------------------------- 2408 //-------------------------------------------------------- 2409 2410 //-------------------------------------------------------- 2411 // Constructor 2412 // Grab references to ATC and thermostat data 2413 //-------------------------------------------------------- ThermostatPowerVerletFiltered(AtomicRegulator * thermostat)2414 ThermostatPowerVerletFiltered::ThermostatPowerVerletFiltered(AtomicRegulator * thermostat) : 2415 ThermostatPowerVerlet(thermostat), 2416 nodalTemperature2Roc_(atc_->field_2roc(TEMPERATURE)), 2417 fieldsRoc_(atc_->fields_roc()), 2418 filterScale_((atc_->time_filter_manager())->filter_scale()) 2419 { 2420 heatSourceRoc_.reset(nNodes_,1); 2421 fluxRoc_[TEMPERATURE].reset(nNodes_,1); 2422 } 2423 2424 //-------------------------------------------------------- 2425 // compute_boundary_flux 2426 // also sets time derivatives of boundary flux and 2427 // heat sources 2428 //-------------------------------------------------------- compute_boundary_flux(FIELDS & fields)2429 void ThermostatPowerVerletFiltered::compute_boundary_flux(FIELDS & fields) 2430 { 2431 ThermostatPowerVerlet::compute_boundary_flux(fields); 2432 2433 // compute boundary flux rate of change 2434 fluxRoc_[TEMPERATURE] = 0.; 2435 atc_->compute_boundary_flux(fieldMask_, 2436 fieldsRoc_, 2437 fluxRoc_, 2438 atomMaterialGroups_, 2439 shpFcnDerivs_); 2440 2441 2442 2443 // compute extrinsic model rate of change 2444 (atc_->extrinsic_model_manager()).set_sources(fieldsRoc_,fluxRoc_); 2445 heatSourceRoc_ = fluxRoc_[TEMPERATURE].quantity(); 2446 } 2447 2448 //-------------------------------------------------------- 2449 // add_to_rhs: 2450 // determines what if any contributions to the 2451 // finite element equations are needed for 2452 // consistency with the thermostat 2453 //-------------------------------------------------------- add_to_rhs(FIELDS & rhs)2454 void ThermostatPowerVerletFiltered::add_to_rhs(FIELDS & rhs) 2455 { 2456 rhs[TEMPERATURE] += lambdaPowerFiltered_->quantity() + boundaryFlux_[TEMPERATURE].quantity(); 2457 } 2458 2459 //-------------------------------------------------------- 2460 // set_thermostat_rhs: 2461 // sets up the right-hand side including boundary 2462 // fluxes (coupling & prescribed), heat sources, and 2463 // fixed (uncoupled) nodes 2464 //-------------------------------------------------------- set_thermostat_rhs(DENS_MAT & rhs_nodes)2465 void ThermostatPowerVerletFiltered::set_thermostat_rhs(DENS_MAT & rhs_nodes) 2466 { 2467 // (a) for flux based : 2468 // form rhs : \int N_I r dV 2469 // vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I 2470 // fluxes are set in ATC transfer 2471 rhs_nodes = heatSource_.quantity() + filterScale_*heatSourceRoc_.quantity(); 2472 2473 // (b) for ess. bcs 2474 // form rhs : {sum_a (N_Ia * v_ia * f_ia) - 0.5*(dtheta/dt)_I} 2475 const DENS_MAT & myNodalAtomicPower(nodalAtomicPower_->quantity()); 2476 const DIAG_MAT & myMdMassMatrix(mdMassMatrix_.quantity()); 2477 const DENS_MAT & myNodalTemperatureRoc(nodalTemperatureRoc_.quantity()); 2478 const DENS_MAT & myNodalTemperature2Roc(nodalTemperature2Roc_.quantity()); 2479 for (int i = 0; i < nNodes_; i++) { 2480 if (prescribedDataMgr_->is_fixed(i,TEMPERATURE,0)) { 2481 rhs_nodes(i,0) = 0.5*(myNodalAtomicPower(i,0) - myMdMassMatrix(i,i)*(myNodalTemperatureRoc(i,0)+ filterScale_*myNodalTemperature2Roc(i,0))); 2482 } 2483 } 2484 } 2485 2486 //-------------------------------------------------------- 2487 // output: 2488 // adds all relevant output to outputData 2489 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)2490 void ThermostatPowerVerletFiltered::output(OUTPUT_LIST & outputData) 2491 { 2492 outputData["lambda"] = &(lambda_->set_quantity()); 2493 outputData["nodalLambdaPower"] = &(lambdaPowerFiltered_->set_quantity()); 2494 } 2495 2496 //-------------------------------------------------------- 2497 //-------------------------------------------------------- 2498 // Class ThermostatHooverVerletFiltered 2499 //-------------------------------------------------------- 2500 //-------------------------------------------------------- 2501 2502 //-------------------------------------------------------- 2503 // Constructor 2504 // Grab references to ATC and thermostat data 2505 //-------------------------------------------------------- ThermostatHooverVerletFiltered(AtomicRegulator * thermostat)2506 ThermostatHooverVerletFiltered::ThermostatHooverVerletFiltered(AtomicRegulator * thermostat) : 2507 ThermostatPowerVerletFiltered(thermostat), 2508 lambdaHoover_(nullptr), 2509 nodalAtomicHooverLambdaPower_(nullptr) 2510 { 2511 // consistent with stage 3 of ATC_Method::initialize 2512 lambdaHoover_ = atomicRegulator_->regulator_data("LambdaHoover",1); 2513 } 2514 2515 //-------------------------------------------------------- 2516 // constructor_transfers 2517 // instantiates or obtains all dependency managed data 2518 //-------------------------------------------------------- construct_transfers()2519 void ThermostatHooverVerletFiltered::construct_transfers() 2520 { 2521 ThermostatPowerVerletFiltered::construct_transfers(); 2522 InterscaleManager & interscaleManager(atc_->interscale_manager()); 2523 2524 FtaShapeFunctionProlongation * atomHooverLambdas = new FtaShapeFunctionProlongation(atc_, 2525 lambdaHoover_, 2526 interscaleManager.per_atom_sparse_matrix("Interpolant")); 2527 interscaleManager.add_per_atom_quantity(atomHooverLambdas, 2528 regulatorPrefix_+"AtomHooverLambda"); 2529 AtomicThermostatForce * atomHooverThermostatForces = new AtomicThermostatForce(atc_,atomHooverLambdas); 2530 interscaleManager.add_per_atom_quantity(atomHooverThermostatForces, 2531 regulatorPrefix_+"AtomHooverThermostatForce"); 2532 SummedAtomicQuantity<double> * atomTotalThermostatForces = 2533 new SummedAtomicQuantity<double>(atc_,atomThermostatForces_,atomHooverThermostatForces); 2534 interscaleManager.add_per_atom_quantity(atomTotalThermostatForces, 2535 regulatorPrefix_+"AtomTotalThermostatForce"); 2536 atomThermostatForces_ = atomTotalThermostatForces; 2537 2538 // transfers dependent on time integration method 2539 DotTwiceKineticEnergy * atomicHooverLambdaPower = 2540 new DotTwiceKineticEnergy(atc_,atomHooverThermostatForces); 2541 interscaleManager.add_per_atom_quantity(atomicHooverLambdaPower, 2542 regulatorPrefix_+"AtomicHooverLambdaPower"); 2543 2544 nodalAtomicHooverLambdaPower_ = new AtfShapeFunctionRestriction(atc_, 2545 atomicHooverLambdaPower, 2546 interscaleManager.per_atom_sparse_matrix("Interpolant")); 2547 interscaleManager.add_dense_matrix(nodalAtomicHooverLambdaPower_, 2548 regulatorPrefix_+"NodalAtomicHooverLambdaPower"); 2549 } 2550 2551 //-------------------------------------------------------- 2552 // compute_thermostat: 2553 // sets up and solves the thermostat equations since 2554 // they are the same at different parts of the time 2555 // step 2556 //-------------------------------------------------------- compute_thermostat(double dt)2557 void ThermostatHooverVerletFiltered::compute_thermostat(double dt) 2558 { 2559 // apply prescribed/extrinsic sources and fixed nodes 2560 ThermostatPowerVerletFiltered::compute_thermostat(0.5*dt); 2561 _nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); // save power from lambda in power-based thermostat 2562 2563 // set up Hoover rhs 2564 set_hoover_rhs(_rhs_); 2565 2566 // solve linear system for lambda 2567 DENS_MAT & myLambda(lambdaHoover_->set_quantity()); 2568 solve_for_lambda(_rhs_,myLambda); 2569 2570 2571 // compute force applied by lambda 2572 // compute nodal atomic power from Hoover coupling 2573 // only add in contribution to uncoupled nodes 2574 if (atomicRegulator_->use_localized_lambda()) 2575 add_to_lambda_power(atomThermostatForces_->quantity(),0.5*dt); 2576 } 2577 2578 //-------------------------------------------------------- 2579 // set_hoover_rhs: 2580 // sets up the right-hand side for fixed value, 2581 // i.e. Hoover coupling 2582 //-------------------------------------------------------- set_hoover_rhs(DENS_MAT & rhs)2583 void ThermostatHooverVerletFiltered::set_hoover_rhs(DENS_MAT & rhs) 2584 { 2585 // form rhs : sum_a (N_Ia * v_ia * f_ia) - 0.5*M_MD*(dtheta/dt)_I 2586 rhs = nodalAtomicPower_->quantity(); 2587 rhs -= mdMassMatrix_.quantity()*(nodalTemperatureRoc_.quantity() + filterScale_*nodalTemperature2Roc_.quantity()); 2588 rhs /= 2.; 2589 } 2590 2591 //-------------------------------------------------------- 2592 // add_to_nodal_lambda_power: 2593 // determines the power exerted by the Hoover 2594 // thermostat at each FE node 2595 //-------------------------------------------------------- add_to_lambda_power(const DENS_MAT &,double dt)2596 void ThermostatHooverVerletFiltered::add_to_lambda_power(const DENS_MAT & /* myLambdaForce */, 2597 double dt) 2598 { 2599 _myNodalLambdaPower_ = nodalAtomicHooverLambdaPower_->quantity(); 2600 const INT_ARRAY nodeToOverlapMap(nodeToOverlapMap_->quantity()); 2601 for (int i = 0; i < nNodes_; ++i) { 2602 if (nodeToOverlapMap(i,0)==-1) 2603 _nodalAtomicLambdaPowerOut_(i,0) += _myNodalLambdaPower_(i,0); 2604 else 2605 _myNodalLambdaPower_(i,0) = 0.; 2606 } 2607 timeFilter_->apply_post_step1(lambdaPowerFiltered_->set_quantity(),_myNodalLambdaPower_,dt); 2608 } 2609 2610 //-------------------------------------------------------- 2611 // add_to_rhs: 2612 // determines what if any contributions to the 2613 // finite element equations are needed for 2614 // consistency with the thermostat 2615 //-------------------------------------------------------- add_to_rhs(FIELDS & rhs)2616 void ThermostatHooverVerletFiltered::add_to_rhs(FIELDS & rhs) 2617 { 2618 rhs[TEMPERATURE] += lambdaPowerFiltered_->quantity(); 2619 } 2620 2621 }; 2622