1 // ATC transfer headers 2 #include "ElasticTimeIntegrator.h" 3 #include "ATC_Coupling.h" 4 #include "TimeFilter.h" 5 #include "ATC_Error.h" 6 #include "PerAtomQuantityLibrary.h" 7 8 namespace ATC { 9 10 //-------------------------------------------------------- 11 //-------------------------------------------------------- 12 // Class MomentumTimeIntegrator 13 //-------------------------------------------------------- 14 //-------------------------------------------------------- 15 16 //-------------------------------------------------------- 17 // Constructor 18 //-------------------------------------------------------- MomentumTimeIntegrator(ATC_Coupling * atc,TimeIntegrationType timeIntegrationType)19 MomentumTimeIntegrator::MomentumTimeIntegrator(ATC_Coupling * atc, 20 TimeIntegrationType timeIntegrationType) : 21 TimeIntegrator(atc, timeIntegrationType) 22 { 23 // do nothing 24 } 25 26 //-------------------------------------------------------- 27 // modify 28 // parses inputs and modifies state of the integrator 29 //-------------------------------------------------------- modify(int,char ** arg)30 bool MomentumTimeIntegrator::modify(int /* narg */, char **arg) 31 { 32 bool foundMatch = false; 33 int argIndex = 0; 34 // time integration scheme 35 /*! \page man_momentum_time_integration fix_modify AtC time_integration (momentum) 36 \section syntax 37 fix_modify AtC time_integration <descriptor> \n 38 - descriptor (string) = time integration type \n 39 40 various time integration methods for the finite elements\n 41 \section description 42 verlet - atomic velocity update with 2nd order Verlet, nodal temperature update with 2nd order Verlet, kinetostats based on controlling force \n 43 fractional_step - atomic velocity update with 2nd order Verlet, mixed nodal momentum update, 2nd order Verlet for continuum and exact 2nd order Verlet for atomic contributions, kinetostats based on controlling discrete momentum changes\n 44 gear - atomic velocity update with 2nd order Verlet, nodal temperature update with 3rd or 4th order Gear, kinetostats based on controlling power \n 45 \section examples 46 <TT> fix_modify atc time_integration verlet </TT> \n 47 <TT> fix_modify atc time_integration fractional_step </TT> \n 48 \section description 49 \section related 50 see \ref man_fix_atc 51 \section default 52 none 53 */ 54 if (strcmp(arg[argIndex],"verlet")==0) { 55 timeIntegrationType_ = VERLET; 56 needReset_ = true; 57 foundMatch = true; 58 } 59 else if (strcmp(arg[argIndex],"fractional_step")==0) { 60 timeIntegrationType_ = FRACTIONAL_STEP; 61 needReset_ = true; 62 foundMatch = true; 63 } 64 else if (strcmp(arg[argIndex],"gear")==0) { 65 timeIntegrationType_ = GEAR; 66 needReset_ = true; 67 foundMatch = true; 68 } 69 return foundMatch; 70 } 71 72 //-------------------------------------------------------- 73 // construct_methods 74 // creates algorithm objects 75 //-------------------------------------------------------- construct_methods()76 void MomentumTimeIntegrator::construct_methods() 77 { 78 if (atc_->reset_methods()) { 79 if (timeIntegrationMethod_) 80 delete timeIntegrationMethod_; 81 82 if (timeFilterManager_->need_reset()) { 83 switch (timeIntegrationType_) { 84 case VERLET: 85 timeFilter_ = timeFilterManager_->construct(TimeFilterManager::IMPLICIT); 86 atc_->set_mass_mat_time_filter(MOMENTUM,TimeFilterManager::IMPLICIT); 87 break; 88 case FRACTIONAL_STEP: 89 case GEAR: 90 timeFilter_ = timeFilterManager_->construct(TimeFilterManager::EXPLICIT_IMPLICIT); 91 atc_->set_mass_mat_time_filter(MOMENTUM,TimeFilterManager::EXPLICIT_IMPLICIT); 92 break; 93 default: 94 throw ATC_Error("Unknown time integration type in ThermalTimeIntegrator::Initialize()"); 95 } 96 } 97 98 if (timeFilterManager_->filter_dynamics()) { 99 switch (timeIntegrationType_) { 100 case VERLET: { 101 timeIntegrationMethod_ = new ElasticTimeIntegratorVerletFiltered(this); 102 break; 103 } 104 default: 105 throw ATC_Error("Unknown time integration type in MomentumTimeIntegrator::Initialize()"); 106 } 107 } 108 else { 109 switch (timeIntegrationType_) { 110 case VERLET: { 111 timeIntegrationMethod_ = new ElasticTimeIntegratorVerlet(this); 112 break; 113 } 114 case FRACTIONAL_STEP: { 115 timeIntegrationMethod_ = new ElasticTimeIntegratorFractionalStep(this); 116 break; 117 } 118 case GEAR: { 119 timeIntegrationMethod_ = new FluidsTimeIntegratorGear(this); 120 break; 121 } 122 default: 123 throw ATC_Error("Unknown time integration type in MomentumTimeIntegrator::Initialize()"); 124 } 125 } 126 } 127 } 128 129 //-------------------------------------------------------- 130 // pack_fields 131 // add persistent variables to data list 132 //-------------------------------------------------------- pack_fields(RESTART_LIST & data)133 void MomentumTimeIntegrator::pack_fields(RESTART_LIST & data) 134 { 135 data["NodalAtomicForceFiltered"] = & nodalAtomicForceFiltered_.set_quantity(); 136 data["NodalAtomicMomentumFiltered"] = & nodalAtomicMomentumFiltered_.set_quantity(); 137 TimeIntegrator::pack_fields(data); 138 } 139 140 //-------------------------------------------------------- 141 //-------------------------------------------------------- 142 // Class MomentumIntegrationMethod 143 //-------------------------------------------------------- 144 //-------------------------------------------------------- 145 146 //-------------------------------------------------------- 147 // Constructor 148 // Grab data from ATC 149 //-------------------------------------------------------- MomentumIntegrationMethod(MomentumTimeIntegrator * momentumTimeIntegrator)150 MomentumIntegrationMethod::MomentumIntegrationMethod(MomentumTimeIntegrator * momentumTimeIntegrator) : 151 TimeIntegrationMethod(momentumTimeIntegrator), 152 timeFilter_(timeIntegrator_->time_filter()), 153 velocity_(atc_->field(VELOCITY)), 154 acceleration_(atc_->field_roc(VELOCITY)), 155 nodalAtomicVelocityOut_(atc_->nodal_atomic_field(VELOCITY)), 156 velocityRhs_(atc_->field_rhs(VELOCITY)), 157 nodalAtomicForceOut_(atc_->nodal_atomic_field_roc(VELOCITY)) 158 { 159 // do nothing 160 } 161 162 //-------------------------------------------------------- 163 // construct_transfers 164 // Grab existing managed quantities, 165 // create the rest 166 //-------------------------------------------------------- construct_transfers()167 void MomentumIntegrationMethod::construct_transfers() 168 { 169 InterscaleManager & interscaleManager(atc_->interscale_manager()); 170 nodalAtomicVelocity_ = interscaleManager.dense_matrix("NodalAtomicVelocity"); 171 } 172 173 //-------------------------------------------------------- 174 //-------------------------------------------------------- 175 // Class ElasticTimeIntegratorVerlet 176 //-------------------------------------------------------- 177 //-------------------------------------------------------- 178 179 //-------------------------------------------------------- 180 // Constructor 181 //-------------------------------------------------------- ElasticTimeIntegratorVerlet(MomentumTimeIntegrator * momentumTimeIntegrator)182 ElasticTimeIntegratorVerlet::ElasticTimeIntegratorVerlet(MomentumTimeIntegrator * momentumTimeIntegrator) : 183 MomentumIntegrationMethod(momentumTimeIntegrator), 184 displacement_(atc_->field(DISPLACEMENT)), 185 nodalAtomicDisplacementOut_(atc_->nodal_atomic_field(DISPLACEMENT)), 186 nodalAtomicForceFiltered_(momentumTimeIntegrator->nodal_atomic_force_filtered()), 187 nodalAtomicDisplacement_(nullptr), 188 nodalAtomicForce_(nullptr) 189 { 190 // do nothing 191 } 192 193 //-------------------------------------------------------- 194 // construct_transfers 195 // Grab existing managed quantities, 196 // create the rest 197 //-------------------------------------------------------- construct_transfers()198 void ElasticTimeIntegratorVerlet::construct_transfers() 199 { 200 MomentumIntegrationMethod::construct_transfers(); 201 InterscaleManager & interscaleManager = atc_->interscale_manager(); 202 nodalAtomicDisplacement_ = interscaleManager.dense_matrix("NodalAtomicDisplacement"); 203 nodalAtomicForce_ = interscaleManager.dense_matrix("NodalAtomicForce"); 204 } 205 206 //-------------------------------------------------------- 207 // initialize 208 // initialize all data 209 //-------------------------------------------------------- initialize()210 void ElasticTimeIntegratorVerlet::initialize() 211 { 212 MomentumIntegrationMethod::initialize(); 213 214 // sets up time filter for cases where variables temporally filtered 215 // this time integrator should use an implicit filter 216 TimeFilterManager * timeFilterManager = (timeIntegrator_->atc())->time_filter_manager(); 217 if (timeFilterManager->need_reset()) { 218 timeFilter_->initialize(nodalAtomicForce_->quantity()); 219 } 220 221 if (!(timeFilterManager->end_equilibrate())) { 222 nodalAtomicForceFiltered_.reset(atc_->num_nodes(),atc_->nsd()); 223 } 224 225 if (!(timeFilterManager->filter_dynamics())){ 226 //post_process(); 227 //compute_nodal_forces(velocityRhs_.set_quantity()); 228 velocityRhs_ = nodalAtomicForce_->quantity(); 229 } 230 } 231 232 //-------------------------------------------------------- 233 // pre_initial_integrate1 234 // time integration before Verlet step 1 235 //-------------------------------------------------------- pre_initial_integrate1(double dt)236 void ElasticTimeIntegratorVerlet::pre_initial_integrate1(double dt) 237 { 238 explicit_1(velocity_.set_quantity(),acceleration_.quantity(),.5*dt); 239 } 240 241 242 //-------------------------------------------------------- 243 // post_initial_integrate1 244 // time integration after Verlet step 1 245 //-------------------------------------------------------- post_initial_integrate1(double dt)246 void ElasticTimeIntegratorVerlet::post_initial_integrate1(double dt) 247 { 248 249 // for improved accuracy, but this would be inconsistent with 250 // the atomic integration scheme 251 explicit_1(displacement_.set_quantity(),velocity_.quantity(),dt); 252 } 253 254 //-------------------------------------------------------- 255 // pre_final_integrate1 256 // first time integration computations 257 // before Verlet step 2 258 //-------------------------------------------------------- pre_final_integrate1(double dt)259 void ElasticTimeIntegratorVerlet::pre_final_integrate1(double dt) 260 { 261 // integrate filtered atomic force 262 timeFilter_->apply_post_step2(nodalAtomicForceFiltered_.set_quantity(), 263 nodalAtomicForce_->quantity(),dt); 264 } 265 266 //-------------------------------------------------------- 267 // post_final_integrate2 268 // second time integration computations 269 // after Verlet step 2 270 //-------------------------------------------------------- post_final_integrate2(double dt)271 void ElasticTimeIntegratorVerlet::post_final_integrate2(double dt) 272 { 273 atc_->apply_inverse_mass_matrix(velocityRhs_.quantity(), 274 acceleration_.set_quantity(), 275 VELOCITY); 276 explicit_1(velocity_.set_quantity(),acceleration_.quantity(),.5*dt); 277 } 278 279 //-------------------------------------------------------- 280 // add_to_rhs 281 // add integrated atomic force contributions 282 //-------------------------------------------------------- add_to_rhs()283 void ElasticTimeIntegratorVerlet::add_to_rhs() 284 { 285 // Compute MD contribution to FEM equation 286 velocityRhs_ += nodalAtomicForce_->quantity(); 287 } 288 289 //-------------------------------------------------------- 290 // post_process 291 // post processing of variables before output 292 //-------------------------------------------------------- post_process()293 void ElasticTimeIntegratorVerlet::post_process() 294 { 295 nodalAtomicDisplacementOut_ = nodalAtomicDisplacement_->quantity(); 296 nodalAtomicVelocityOut_ = nodalAtomicVelocity_->quantity(); 297 } 298 299 //-------------------------------------------------------- 300 // output 301 // add variables to output list 302 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)303 void ElasticTimeIntegratorVerlet::output(OUTPUT_LIST & outputData) 304 { 305 DENS_MAT & nodalAtomicForce(nodalAtomicForce_->set_quantity()); 306 if ((atc_->lammps_interface())->rank_zero()) { 307 outputData["NodalAtomicForce"] = &nodalAtomicForce; 308 } 309 } 310 311 //-------------------------------------------------------- 312 // finish 313 // finalize state of nodal atomic quantities 314 //-------------------------------------------------------- finish()315 void ElasticTimeIntegratorVerlet::finish() 316 { 317 post_process(); 318 } 319 320 //-------------------------------------------------------- 321 //-------------------------------------------------------- 322 // Class ElasticTimeIntegratorVerletFiltered 323 //-------------------------------------------------------- 324 //-------------------------------------------------------- 325 326 //-------------------------------------------------------- 327 // Constructor 328 //-------------------------------------------------------- ElasticTimeIntegratorVerletFiltered(MomentumTimeIntegrator * momentumTimeIntegrator)329 ElasticTimeIntegratorVerletFiltered::ElasticTimeIntegratorVerletFiltered(MomentumTimeIntegrator * momentumTimeIntegrator) : 330 ElasticTimeIntegratorVerlet(momentumTimeIntegrator), 331 nodalAtomicAcceleration_(atc_->nodal_atomic_field_roc(VELOCITY)) 332 { 333 // do nothing 334 } 335 336 //-------------------------------------------------------- 337 // pre_initial_integrate1 338 // time integration before Verlet step 1 339 //-------------------------------------------------------- pre_initial_integrate1(double dt)340 void ElasticTimeIntegratorVerletFiltered::pre_initial_integrate1(double dt) 341 { 342 explicit_1(velocity_.set_quantity(),acceleration_.quantity(),.5*dt); 343 explicit_1(nodalAtomicVelocityOut_.set_quantity(),nodalAtomicAcceleration_.quantity(),.5*dt); 344 } 345 346 //-------------------------------------------------------- 347 // post_initial_integrate1 348 // time integration after Verlet step 1 349 //-------------------------------------------------------- post_initial_integrate1(double dt)350 void ElasticTimeIntegratorVerletFiltered::post_initial_integrate1(double dt) 351 { 352 353 // for improved accuracy, but this would be inconsistent with 354 // the atomic integration scheme 355 explicit_1(displacement_.set_quantity(),velocity_.quantity(),dt); 356 explicit_1(nodalAtomicDisplacementOut_.set_quantity(),nodalAtomicVelocityOut_.quantity(),dt); 357 } 358 359 //-------------------------------------------------------- 360 // post_final_integrate2 361 // second time integration after Verlet step 2 362 //-------------------------------------------------------- post_final_integrate2(double dt)363 void ElasticTimeIntegratorVerletFiltered::post_final_integrate2(double dt) 364 { 365 DENS_MAT velocityRoc(velocityRhs_.nRows(),velocityRhs_.nCols()); 366 atc_->apply_inverse_mass_matrix(velocityRhs_.quantity(), 367 acceleration_.set_quantity(), 368 VELOCITY); 369 explicit_1(velocity_.set_quantity(),acceleration_.quantity(),.5*dt); 370 371 atc_->apply_inverse_md_mass_matrix(nodalAtomicForceFiltered_.quantity(), 372 nodalAtomicAcceleration_.set_quantity(), 373 VELOCITY); 374 explicit_1(nodalAtomicVelocityOut_.set_quantity(),nodalAtomicAcceleration_.quantity(),.5*dt); 375 } 376 377 //-------------------------------------------------------- 378 // add_to_rhs 379 // add integrated atomic force contributions 380 //-------------------------------------------------------- add_to_rhs()381 void ElasticTimeIntegratorVerletFiltered::add_to_rhs() 382 { 383 // MD contributions to FE equations 384 velocityRhs_ += nodalAtomicForceFiltered_.set_quantity(); 385 } 386 387 //-------------------------------------------------------- 388 // output 389 // add variables to output list 390 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)391 void ElasticTimeIntegratorVerletFiltered::output(OUTPUT_LIST & outputData) 392 { 393 DENS_MAT & nodalAtomicForce(nodalAtomicForceFiltered_.set_quantity()); 394 if ((atc_->lammps_interface())->rank_zero()) { 395 outputData["NodalAtomicForce"] = &nodalAtomicForce; 396 } 397 } 398 399 //-------------------------------------------------------- 400 //-------------------------------------------------------- 401 // Class ElasticTimeIntegratorFractionalStep 402 //-------------------------------------------------------- 403 //-------------------------------------------------------- 404 405 //-------------------------------------------------------- 406 // Constructor 407 //-------------------------------------------------------- ElasticTimeIntegratorFractionalStep(MomentumTimeIntegrator * momentumTimeIntegrator)408 ElasticTimeIntegratorFractionalStep::ElasticTimeIntegratorFractionalStep(MomentumTimeIntegrator * momentumTimeIntegrator) : 409 MomentumIntegrationMethod(momentumTimeIntegrator), 410 displacement_(atc_->field(DISPLACEMENT)), 411 nodalAtomicDisplacementOut_(atc_->nodal_atomic_field(DISPLACEMENT)), 412 nodalAtomicForceFiltered_(momentumTimeIntegrator->nodal_atomic_force_filtered()), 413 nodalAtomicMomentum_(nullptr), 414 nodalAtomicMomentumFiltered_(momentumTimeIntegrator->nodal_atomic_momentum_filtered()), 415 nodalAtomicDisplacement_(nullptr), 416 nodalAtomicMomentumOld_(atc_->num_nodes(),atc_->nsd()), 417 nodalAtomicVelocityOld_(atc_->num_nodes(),atc_->nsd()) 418 { 419 // do nothing 420 } 421 422 //-------------------------------------------------------- 423 // construct_transfers 424 // Grab existing managed quantities, 425 // create the rest 426 //-------------------------------------------------------- construct_transfers()427 void ElasticTimeIntegratorFractionalStep::construct_transfers() 428 { 429 MomentumIntegrationMethod::construct_transfers(); 430 InterscaleManager & interscaleManager = atc_->interscale_manager(); 431 nodalAtomicMomentum_ = interscaleManager.dense_matrix("NodalAtomicMomentum"); 432 nodalAtomicDisplacement_ = interscaleManager.dense_matrix("NodalAtomicDisplacement"); 433 } 434 435 //-------------------------------------------------------- 436 // initialize 437 // initialize all data 438 //-------------------------------------------------------- initialize()439 void ElasticTimeIntegratorFractionalStep::initialize() 440 { 441 MomentumIntegrationMethod::initialize(); 442 443 // initial force to zero 444 nodalAtomicForce_.reset(atc_->num_nodes(),atc_->nsd()); 445 446 // sets up time filter for cases where variables temporally filtered 447 // this time integrator should use Crank-Nicholson filter for 2nd order accuracy 448 TimeFilterManager * timeFilterManager = (timeIntegrator_->atc())->time_filter_manager(); 449 if (timeFilterManager->need_reset()) { 450 // the form of this integrator implies no time filters that require history data can be used 451 timeFilter_->initialize(); 452 } 453 454 // sets up time filter for post-processing the filtered power 455 // this time integrator should use an explicit-implicit filter 456 // to mirror the 2nd order Verlet integration scheme 457 // It requires no history information so initial value just sizes arrays 458 if (!(timeFilterManager->end_equilibrate())) { 459 // implies an initial condition of the instantaneous atomic energy 460 nodalAtomicMomentumFiltered_ = nodalAtomicMomentum_->quantity(); 461 nodalAtomicForceFiltered_.reset(atc_->num_nodes(),atc_->nsd()); 462 } 463 } 464 465 //-------------------------------------------------------- 466 // pre_initial_integrate1 467 // time integration before Verlet step 1, used to 468 // provide the baseline momentum and displacement to 469 // quantify the change over the timestep 470 //-------------------------------------------------------- pre_initial_integrate1(double dt)471 void ElasticTimeIntegratorFractionalStep::pre_initial_integrate1(double dt) 472 { 473 // initialize changes in momentum 474 const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity()); 475 // updated filtered energy using explicit-implicit scheme 476 timeFilter_->apply_pre_step1(nodalAtomicMomentumFiltered_.set_quantity(), 477 myNodalAtomicMomentum,dt); 478 } 479 480 //-------------------------------------------------------- 481 // pre_initial_integrate2 482 // second time integration after kinetostat application 483 // to compute MD contributions to momentum change 484 //-------------------------------------------------------- pre_initial_integrate2(double dt)485 void ElasticTimeIntegratorFractionalStep::pre_initial_integrate2(double dt) 486 { 487 // used for updating change in velocity from mass matrix change 488 this->compute_old_time_data(); 489 490 // update filtered nodal atomic force 491 timeFilter_->apply_pre_step1(nodalAtomicForceFiltered_.set_quantity(), 492 nodalAtomicForce_,dt); 493 494 // store current force for use later 495 nodalAtomicForce_ = nodalAtomicMomentum_->quantity(); 496 nodalAtomicForce_ *= -1.; 497 } 498 499 //-------------------------------------------------------- 500 // post_initial_integrate1 501 // time integration after Verlet step 1 502 //-------------------------------------------------------- post_initial_integrate1(double dt)503 void ElasticTimeIntegratorFractionalStep::post_initial_integrate1(double dt) 504 { 505 // atomic contributions to change in momentum 506 // compute change in restricted atomic momentum 507 const DENS_MAT & nodalAtomicMomentum(nodalAtomicMomentum_->quantity()); 508 nodalAtomicForce_ += nodalAtomicMomentum; 509 510 // update FE velocity with change in velocity from MD 511 DENS_MAT & atomicVelocityDelta(atomicVelocityDelta_.set_quantity()); 512 atc_->apply_inverse_mass_matrix(nodalAtomicForce_, 513 atomicVelocityDelta, 514 VELOCITY); 515 velocity_ += atomicVelocityDelta; 516 517 // approximation to force for output 518 nodalAtomicForce_ /= 0.5*dt; 519 timeFilter_->apply_post_step1(nodalAtomicForceFiltered_.set_quantity(), 520 nodalAtomicForce_,dt); 521 522 // change to velocity from FE dynamics 523 atc_->apply_inverse_mass_matrix(velocityRhs_.quantity(), 524 acceleration_.set_quantity(), 525 VELOCITY); 526 explicit_1(velocity_.set_quantity(),acceleration_.quantity(),0.5*dt); 527 528 // used for updating change in momentum from mass matrix change 529 atc_->apply_inverse_mass_matrix(nodalAtomicMomentum, 530 nodalAtomicVelocityOld_, 531 VELOCITY); 532 nodalAtomicMomentumOld_ = nodalAtomicMomentum; 533 534 // get nodal momentum for second part of force update 535 nodalAtomicForce_ = nodalAtomicMomentum; 536 nodalAtomicForce_ *= -1.; 537 538 // update nodal displacements 539 explicit_1(displacement_.set_quantity(),velocity_.quantity(),dt); 540 } 541 542 //-------------------------------------------------------- 543 // post_final_integrate2 544 // second time integration computations 545 // after Verlet step 2 546 //-------------------------------------------------------- post_final_integrate2(double dt)547 void ElasticTimeIntegratorFractionalStep::post_final_integrate2(double dt) 548 { 549 // atomic contributions to change in momentum 550 // compute change in restricted atomic momentum 551 nodalAtomicForce_ += nodalAtomicMomentum_->quantity(); 552 553 // update FE temperature with change in temperature from MD 554 compute_velocity_delta(nodalAtomicForce_,dt); 555 velocity_ += atomicVelocityDelta_.quantity(); 556 557 // approximation to power for output 558 nodalAtomicForce_ /= 0.5*dt; 559 timeFilter_->apply_post_step1(nodalAtomicForceFiltered_.set_quantity(), 560 nodalAtomicForce_,dt); 561 562 // change to velocity from FE dynamics 563 atc_->apply_inverse_mass_matrix(velocityRhs_.quantity(), 564 acceleration_.set_quantity(), 565 VELOCITY); 566 explicit_1(velocity_.set_quantity(),acceleration_.quantity(),0.5*dt); 567 } 568 //-------------------------------------------------------- 569 // post_process 570 // post processing of variables before output 571 //-------------------------------------------------------- post_process()572 void ElasticTimeIntegratorFractionalStep::post_process() 573 { 574 nodalAtomicDisplacementOut_ = nodalAtomicDisplacement_->quantity(); 575 nodalAtomicVelocityOut_ = nodalAtomicVelocity_->quantity(); 576 } 577 578 //-------------------------------------------------------- 579 // output 580 // add variables to output list 581 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)582 void ElasticTimeIntegratorFractionalStep::output(OUTPUT_LIST & outputData) 583 { 584 if ((atc_->lammps_interface())->rank_zero()) { 585 outputData["NodalAtomicForce"] = & nodalAtomicForce_; 586 } 587 } 588 589 //-------------------------------------------------------- 590 // finish 591 // finalize state of nodal atomic quantities 592 //-------------------------------------------------------- finish()593 void ElasticTimeIntegratorFractionalStep::finish() 594 { 595 post_process(); 596 } 597 598 //-------------------------------------------------------- 599 // compute_old_time_data 600 //-------------------------------------------------------- compute_old_time_data()601 void ElasticTimeIntegratorFractionalStep::compute_old_time_data() 602 { 603 const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity()); 604 atc_->apply_inverse_mass_matrix(myNodalAtomicMomentum, 605 nodalAtomicVelocityOld_, 606 VELOCITY); 607 nodalAtomicMomentumOld_ = myNodalAtomicMomentum; 608 } 609 610 //-------------------------------------------------------- 611 // compute_velocity_delta 612 //-------------------------------------------------------- compute_velocity_delta(const DENS_MAT & atomicMomentumDelta,double)613 void ElasticTimeIntegratorFractionalStep::compute_velocity_delta(const DENS_MAT & atomicMomentumDelta, 614 double /* dt */) 615 { 616 DENS_MAT & myAtomicVelocityDelta(atomicVelocityDelta_.set_quantity()); 617 myAtomicVelocityDelta = nodalAtomicMomentumOld_ + atomicMomentumDelta; 618 atc_->apply_inverse_mass_matrix(myAtomicVelocityDelta, 619 VELOCITY); 620 myAtomicVelocityDelta += -1.*nodalAtomicVelocityOld_; 621 } 622 //-------------------------------------------------------- 623 //-------------------------------------------------------- 624 // Class FluidsTimeIntegratorGear 625 //-------------------------------------------------------- 626 //-------------------------------------------------------- 627 628 //-------------------------------------------------------- 629 // Constructor 630 // Grab data from ATC 631 //-------------------------------------------------------- 632 FluidsTimeIntegratorGear(MomentumTimeIntegrator * momentumTimeIntegrator)633 FluidsTimeIntegratorGear::FluidsTimeIntegratorGear(MomentumTimeIntegrator * momentumTimeIntegrator) : 634 MomentumIntegrationMethod(momentumTimeIntegrator), 635 nodalAtomicForceFiltered_(momentumTimeIntegrator->nodal_atomic_force_filtered()), 636 nodalAtomicMomentum_(nullptr), 637 nodalAtomicMomentumFiltered_(momentumTimeIntegrator->nodal_atomic_momentum_filtered()), 638 atomicVelocityDelta_(atc_->num_nodes(),atc_->nsd()), 639 nodalAtomicMomentumOld_(atc_->num_nodes(),atc_->nsd()), 640 nodalAtomicVelocityOld_(atc_->num_nodes(),atc_->nsd()), 641 velocity2Roc_(atc_->field_2roc(VELOCITY)) 642 { 643 // do nothing 644 } 645 646 //-------------------------------------------------------- 647 // construct_transfers 648 // Grab existing managed quantities, 649 // create the rest 650 //-------------------------------------------------------- construct_transfers()651 void FluidsTimeIntegratorGear::construct_transfers() 652 { 653 MomentumIntegrationMethod::construct_transfers(); 654 InterscaleManager & interscaleManager(atc_->interscale_manager()); 655 nodalAtomicMomentum_ = interscaleManager.dense_matrix("NodalAtomicMomentum"); 656 } 657 658 //-------------------------------------------------------- 659 // initialize 660 // initialize all data 661 //-------------------------------------------------------- initialize()662 void FluidsTimeIntegratorGear::initialize() 663 { 664 MomentumIntegrationMethod::initialize(); 665 666 // initial power to zero 667 nodalAtomicForce_.reset(atc_->num_nodes(),atc_->nsd()); 668 669 // sets up time filter for cases where variables temporally filtered 670 // this time integrator should use Crank-Nicholson filter for 2nd order accuracy 671 TimeFilterManager * timeFilterManager = atc_->time_filter_manager(); 672 if (timeFilterManager->need_reset()) { 673 // the form of this integrator implies no time filters that require history data can be used 674 timeFilter_->initialize(); 675 } 676 677 // sets up time filter for post-processing the filtered power 678 // this time integrator should use an explicit-implicit filter 679 // to mirror the 2nd order Verlet integration scheme 680 // It requires no history information so initial value just sizes arrays 681 if (!timeFilterManager->end_equilibrate()) { 682 // implies an initial condition of the instantaneous atomic energy 683 // for the corresponding filtered variable, consistent with the temperature 684 nodalAtomicMomentumFiltered_ = nodalAtomicMomentum_->quantity(); 685 nodalAtomicForceFiltered_.reset(atc_->num_nodes(),atc_->nsd()); 686 } 687 } 688 689 //-------------------------------------------------------- 690 // pre_initial_integrate1 691 //-------------------------------------------------------- pre_initial_integrate1(double dt)692 void FluidsTimeIntegratorGear::pre_initial_integrate1(double dt) 693 { 694 const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity()); 695 // updated filtered momentum using explicit-implicit scheme 696 timeFilter_->apply_pre_step1(nodalAtomicMomentumFiltered_.set_quantity(), 697 myNodalAtomicMomentum,dt); 698 } 699 700 //-------------------------------------------------------- 701 // pre_initial_integrate2 702 //-------------------------------------------------------- pre_initial_integrate2(double dt)703 void FluidsTimeIntegratorGear::pre_initial_integrate2(double dt) 704 { 705 // used for updating change in velocity from mass matrix change 706 this->compute_old_time_data(); 707 708 // update FE contributions 709 apply_gear_predictor(dt); 710 711 // update filtered nodal atomic force 712 713 // that way kinetostat and integrator can be consistent 714 timeFilter_->apply_pre_step1(nodalAtomicForceFiltered_.set_quantity(), 715 nodalAtomicForce_,dt); 716 717 // store current momentum for use later 718 nodalAtomicForce_ = nodalAtomicMomentum_->quantity(); 719 nodalAtomicForce_ *= -1.; 720 } 721 722 //-------------------------------------------------------- 723 // pre_final_integrate1 724 //-------------------------------------------------------- pre_final_integrate1(double dt)725 void FluidsTimeIntegratorGear::pre_final_integrate1(double dt) 726 { 727 728 // before the new rhs is computed but after atomic velocity is updated. 729 // compute change in restricted atomic momentum 730 nodalAtomicForce_ += nodalAtomicMomentum_->quantity(); 731 732 // update FE velocity with change in velocity from MD 733 compute_velocity_delta(nodalAtomicForce_,dt); 734 velocity_ += atomicVelocityDelta_.quantity(); 735 736 // approximation to force for output 737 nodalAtomicForce_ /= dt; 738 timeFilter_->apply_post_step1(nodalAtomicForceFiltered_.set_quantity(), 739 nodalAtomicForce_,dt); 740 741 // make sure nodes are fixed 742 atc_->set_fixed_nodes(); 743 } 744 745 //-------------------------------------------------------- 746 // post_final_integrate1 747 //-------------------------------------------------------- post_final_integrate1(double dt)748 void FluidsTimeIntegratorGear::post_final_integrate1(double dt) 749 { 750 // Finish updating temperature with FE contributions 751 atc_->apply_inverse_mass_matrix(velocityRhs_.quantity(), 752 _velocityResidual_,VELOCITY); 753 _velocityResidual_ -= acceleration_.quantity(); 754 _velocityResidual_ *= dt; 755 apply_gear_corrector(_velocityResidual_,dt); 756 } 757 758 //-------------------------------------------------------- 759 // post_process 760 //-------------------------------------------------------- post_final_integrate2(double dt)761 void FluidsTimeIntegratorGear::post_final_integrate2(double dt) 762 { 763 // update filtered atomic energy 764 timeFilter_->apply_post_step1(nodalAtomicMomentumFiltered_.set_quantity(), 765 nodalAtomicMomentum_->quantity(),dt); 766 } 767 768 //-------------------------------------------------------- 769 // output 770 // add variables to output list 771 //-------------------------------------------------------- post_process()772 void FluidsTimeIntegratorGear::post_process() 773 { 774 nodalAtomicForceOut_ = nodalAtomicForce_; 775 nodalAtomicVelocityOut_ = nodalAtomicVelocity_->quantity(); 776 } 777 778 //-------------------------------------------------------- 779 // output 780 // add variables to output list 781 //-------------------------------------------------------- output(OUTPUT_LIST & outputData)782 void FluidsTimeIntegratorGear::output(OUTPUT_LIST & outputData) 783 { 784 if ((atc_->lammps_interface())->rank_zero()) { 785 outputData["NodalAtomicForce"] = & nodalAtomicForce_; 786 } 787 } 788 789 //-------------------------------------------------------- 790 // finish 791 // finalize state of nodal atomic quantities 792 //-------------------------------------------------------- finish()793 void FluidsTimeIntegratorGear::finish() 794 { 795 post_process(); 796 } 797 798 //-------------------------------------------------------- 799 // apply_gear_predictor 800 //-------------------------------------------------------- apply_gear_predictor(double dt)801 void FluidsTimeIntegratorGear::apply_gear_predictor(double dt) 802 { 803 gear1_3_predict(velocity_.set_quantity(), 804 acceleration_.set_quantity(), 805 velocity2Roc_.quantity(),dt); 806 } 807 808 //-------------------------------------------------------- 809 // apply_gear_corrector 810 //-------------------------------------------------------- apply_gear_corrector(const DENS_MAT & residual,double dt)811 void FluidsTimeIntegratorGear::apply_gear_corrector(const DENS_MAT & residual, double dt) 812 { 813 gear1_3_correct(velocity_.set_quantity(), 814 acceleration_.set_quantity(), 815 velocity2Roc_.set_quantity(), 816 residual,dt); 817 } 818 819 //-------------------------------------------------------- 820 // compute_old_time_data 821 //-------------------------------------------------------- compute_old_time_data()822 void FluidsTimeIntegratorGear::compute_old_time_data() 823 { 824 const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity()); 825 atc_->apply_inverse_mass_matrix(myNodalAtomicMomentum, 826 nodalAtomicVelocityOld_, 827 VELOCITY); 828 nodalAtomicMomentumOld_ = myNodalAtomicMomentum; 829 } 830 831 //-------------------------------------------------------- 832 // compute_velocity_delta 833 //-------------------------------------------------------- compute_velocity_delta(const DENS_MAT & atomicMomentumDelta,double)834 void FluidsTimeIntegratorGear::compute_velocity_delta(const DENS_MAT & atomicMomentumDelta, 835 double /* dt */) 836 { 837 DENS_MAT & myAtomicVelocityDelta(atomicVelocityDelta_.set_quantity()); 838 myAtomicVelocityDelta = nodalAtomicMomentumOld_ + atomicMomentumDelta; 839 atc_->apply_inverse_mass_matrix(myAtomicVelocityDelta, 840 VELOCITY); 841 myAtomicVelocityDelta -= nodalAtomicVelocityOld_; 842 } 843 }; 844