1 // ATC_Transfer headers 2 #include "ATC_Transfer.h" 3 #include "ATC_Error.h" 4 #include "FE_Engine.h" 5 #include "LammpsInterface.h" 6 #include "Quadrature.h" 7 #include "VoigtOperations.h" 8 #include "TransferLibrary.h" 9 #include "Stress.h" 10 #include "KernelFunction.h" 11 #include "PerPairQuantity.h" 12 #include "FieldManager.h" 13 #define ESHELBY_VIRIAL 14 #include "LinearSolver.h" 15 16 17 // Other Headers 18 #include <vector> 19 #include <map> 20 #include <set> 21 #include <utility> 22 #include <fstream> 23 #include <sstream> 24 #include <exception> 25 26 27 28 // PLAN: 29 //* energies 30 //* filters - make filterFields class 31 //* output directly 32 //* enum, tagged, computes, mat(field to field) functions 33 //* grads & rates 34 //* on-the-fly 35 // * remove derived classes 36 37 using namespace std; 38 using namespace ATC_Utility; 39 using namespace voigt3; 40 41 namespace ATC { 42 43 const int numFields_ = 17; 44 FieldName indices_[numFields_] = { 45 CHARGE_DENSITY, 46 MASS_DENSITY, 47 SPECIES_CONCENTRATION, 48 NUMBER_DENSITY, 49 MOMENTUM, 50 VELOCITY, 51 PROJECTED_VELOCITY, 52 DISPLACEMENT, 53 POTENTIAL_ENERGY, 54 KINETIC_ENERGY, 55 KINETIC_TEMPERATURE, 56 TEMPERATURE, 57 CHARGE_FLUX, 58 SPECIES_FLUX, 59 THERMAL_ENERGY, 60 ENERGY, 61 INTERNAL_ENERGY 62 }; 63 //KINETIC_STRESS; 64 //ELECTRIC_POTENTIAL}; 65 ATC_Transfer(string groupName,double ** & perAtomArray,LAMMPS_NS::Fix * thisFix,string matParamFile)66 ATC_Transfer::ATC_Transfer(string groupName, 67 double ** & perAtomArray, 68 LAMMPS_NS::Fix * thisFix, 69 string matParamFile) 70 : ATC_Method(groupName,perAtomArray,thisFix), 71 xPointer_(NULL), 72 outputStepZero_(true), 73 neighborReset_(false), 74 pairMap_(NULL), 75 bondMatrix_(NULL), 76 pairVirial_(NULL), 77 pairHeatFlux_(NULL), 78 nComputes_(0), 79 hasPairs_(true), 80 hasBonds_(false), 81 resetKernelFunction_(false), 82 dxaExactMode_(true), 83 cauchyBornStress_(NULL) 84 { 85 nTypes_ = lammpsInterface_->ntypes(); 86 87 peScale_=1.; 88 keScale_= lammpsInterface_->mvv2e(); 89 // if surrogate model of md (no physics model created) 90 if (matParamFile != "none") { 91 fstream fileId(matParamFile.c_str(), std::ios::in); 92 if (!fileId.is_open()) throw ATC_Error("cannot open material file"); 93 CbData cb; 94 LammpsInterface *lmp = LammpsInterface::instance(); 95 lmp->lattice(cb.cell_vectors, cb.basis_vectors); 96 cb.inv_atom_volume = 1.0 / lmp->volume_per_atom(); 97 cb.e2mvv = 1.0 / lmp->mvv2e(); 98 cb.atom_mass = lmp->atom_mass(1); 99 cb.boltzmann = lmp->boltz(); 100 cb.hbar = lmp->hbar(); 101 cauchyBornStress_ = new StressCauchyBorn(fileId, cb); 102 } 103 104 // Defaults 105 set_time(); 106 107 outputFlags_.reset(NUM_TOTAL_FIELDS); 108 outputFlags_ = false; 109 fieldFlags_.reset(NUM_TOTAL_FIELDS); 110 fieldFlags_ = false; 111 gradFlags_.reset(NUM_TOTAL_FIELDS); 112 gradFlags_ = false; 113 rateFlags_.reset(NUM_TOTAL_FIELDS); 114 rateFlags_ = false; 115 116 outputFields_.resize(NUM_TOTAL_FIELDS); 117 for (int i = 0; i < NUM_TOTAL_FIELDS; i++) { outputFields_[i] = NULL; } 118 119 // Hardy requires ref positions for processor ghosts for bond list 120 121 //needXrefProcessorGhosts_ = true; 122 } 123 124 //------------------------------------------------------------------- ~ATC_Transfer()125 ATC_Transfer::~ATC_Transfer() 126 { 127 interscaleManager_.clear(); 128 if (cauchyBornStress_) delete cauchyBornStress_; 129 } 130 131 //------------------------------------------------------------------- 132 // called before the beginning of a "run" initialize()133 void ATC_Transfer::initialize() 134 { 135 if (kernelOnTheFly_ && !readRefPE_ && !setRefPEvalue_) { 136 if (setRefPE_) { 137 stringstream ss; 138 ss << "WARNING: Reference PE requested from atoms, but not yet implemented for on-the-fly, ignoring"; 139 lammpsInterface_->print_msg_once(ss.str()); 140 setRefPE_ = false; 141 } 142 } 143 144 ATC_Method::initialize(); 145 146 if (!initialized_) { 147 if (cauchyBornStress_) cauchyBornStress_->initialize(); 148 } 149 150 if (!initialized_ || ATC::LammpsInterface::instance()->atoms_sorted() || resetKernelFunction_) { 151 // initialize kernel funciton matrix N_Ia 152 if (! kernelOnTheFly_) { 153 try{ 154 if (!moleculeIds_.empty()) compute_kernel_matrix_molecule(); //KKM add 155 } 156 catch(bad_alloc&) { 157 ATC::LammpsInterface::instance()->print_msg("kernel will be computed on-the-fly"); 158 kernelOnTheFly_ = true; 159 } 160 } 161 resetKernelFunction_ = false; 162 } 163 // clears need for reset 164 ghostManager_.initialize(); 165 166 // initialize bond matrix B_Iab 167 if ((! bondOnTheFly_) 168 && ( ( fieldFlags_(STRESS) 169 || fieldFlags_(ESHELBY_STRESS) 170 || fieldFlags_(HEAT_FLUX) ) ) ) { 171 try { 172 compute_bond_matrix(); 173 } 174 catch(bad_alloc&) { 175 ATC::LammpsInterface::instance()->print_msg("stress/heat_flux will be computed on-the-fly"); 176 177 bondOnTheFly_ = true; 178 } 179 } 180 181 // set sample frequency to output if sample has not be specified 182 if (sampleFrequency_ == 0) sampleFrequency_ = outputFrequency_; 183 184 // output for step 0 185 if (!initialized_) { 186 if (outputFrequency_ > 0) { 187 // initialize filtered data 188 compute_fields(); 189 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 190 if(fieldFlags_(index)) { 191 string name = field_to_string((FieldName) index); 192 filteredData_[name] = hardyData_[name]; 193 timeFilters_(index)->initialize(filteredData_[name].quantity()); 194 } 195 if (rateFlags_(index)) { 196 string name = field_to_string((FieldName) index); 197 string rate_field = name + "_rate"; 198 filteredData_[rate_field] = hardyData_[rate_field]; 199 } 200 if (gradFlags_(index)) { 201 string name = field_to_string((FieldName) index); 202 string grad_field = name + "_gradient"; 203 filteredData_[grad_field] = hardyData_[grad_field]; 204 } 205 } 206 int index = NUM_TOTAL_FIELDS; 207 map <string,int>::const_iterator iter; 208 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 209 string tag = iter->first; 210 filteredData_[tag] = hardyData_[tag]; 211 timeFilters_(index)->initialize(filteredData_[tag].quantity()); 212 #ifdef ESHELBY_VIRIAL 213 if (tag == "virial" && fieldFlags_(ESHELBY_STRESS)) { 214 filteredData_["eshelby_virial"] = hardyData_["eshelby_virial"]; 215 } 216 #endif 217 index++; 218 } 219 output(); 220 } 221 } 222 223 initialized_ = true; 224 225 lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+sampleFrequency_); 226 227 228 //remap_ghost_ref_positions(); 229 update_peratom_output(); 230 } 231 232 //------------------------------------------------------------------- set_continuum_data()233 void ATC_Transfer::set_continuum_data() 234 { 235 ATC_Method::set_continuum_data(); 236 if (!initialized_) { 237 nNodesGlobal_ = feEngine_->fe_mesh()->num_nodes(); 238 } 239 } 240 241 //------------------------------------------------------------------- construct_time_integration_data()242 void ATC_Transfer::construct_time_integration_data() 243 { 244 if (!initialized_) { 245 246 // size arrays for requested/required fields 247 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 248 if (fieldFlags_(index)) { 249 250 int size = FieldSizes[index]; 251 if (atomToElementMapType_ == EULERIAN) { 252 if (index == STRESS) size=6; 253 if (index == CAUCHY_BORN_STRESS) size=6; 254 } 255 if (size == 0) { 256 if (index == SPECIES_CONCENTRATION) size=typeList_.size()+groupList_.size(); 257 } 258 string name = field_to_string((FieldName) index); 259 hardyData_ [name].reset(nNodes_,size); 260 filteredData_[name].reset(nNodes_,size); 261 } 262 } 263 264 // size arrays for projected compute fields 265 map <string,int>::const_iterator iter; 266 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 267 string tag = iter->first; 268 COMPUTE_POINTER cmpt = lammpsInterface_->compute_pointer(tag); 269 int ncols = lammpsInterface_->compute_ncols_peratom(cmpt); 270 hardyData_ [tag].reset(nNodes_,ncols); 271 filteredData_[tag].reset(nNodes_,ncols); 272 #ifdef ESHELBY_VIRIAL 273 if (tag == "virial" && fieldFlags_(ESHELBY_STRESS)) { 274 string esh = "eshelby_virial"; 275 int size = FieldSizes[ESHELBY_STRESS]; 276 hardyData_ [esh].reset(nNodes_,size); 277 filteredData_[esh].reset(nNodes_,size); 278 } 279 #endif 280 } 281 } 282 } 283 //-------------------------------------------------------- 284 // set_computational_geometry 285 // constructs needed transfer operators which define 286 // hybrid atom/FE computational geometry 287 //-------------------------------------------------------- set_computational_geometry()288 void ATC_Transfer::set_computational_geometry() 289 { 290 ATC_Method::set_computational_geometry(); 291 } 292 293 //------------------------------------------------------------------- 294 // construct_interpolant 295 // constructs: interpolatn, accumulant, weights, and spatial derivatives 296 //-------------------------------------------------------- construct_interpolant()297 void ATC_Transfer::construct_interpolant() 298 { 299 // interpolant 300 if (!(kernelOnTheFly_)) { 301 // finite element shape functions for interpolants 302 PerAtomShapeFunction * atomShapeFunctions = new PerAtomShapeFunction(this); 303 interscaleManager_.add_per_atom_sparse_matrix(atomShapeFunctions,"Interpolant"); 304 shpFcn_ = atomShapeFunctions; 305 } 306 // accummulant and weights 307 308 this->create_atom_volume(); 309 // accumulants 310 if (kernelFunction_) { 311 // kernel-based accumulants 312 if (kernelOnTheFly_) { 313 ConstantQuantity<double> * atomCount = new ConstantQuantity<double>(this,1.); 314 interscaleManager_.add_per_atom_quantity(atomCount,"AtomCount"); 315 OnTheFlyKernelAccumulation * myWeights 316 = new OnTheFlyKernelAccumulation(this, 317 atomCount, kernelFunction_, atomCoarseGrainingPositions_); 318 interscaleManager_.add_dense_matrix(myWeights, 319 "KernelInverseWeights"); 320 accumulantWeights_ = new OnTheFlyKernelWeights(myWeights); 321 } 322 else { 323 PerAtomKernelFunction * atomKernelFunctions = new PerAtomKernelFunction(this); 324 interscaleManager_.add_per_atom_sparse_matrix(atomKernelFunctions, 325 "Accumulant"); 326 accumulant_ = atomKernelFunctions; 327 accumulantWeights_ = new AccumulantWeights(accumulant_); 328 } 329 accumulantInverseVolumes_ = new KernelInverseVolumes(this,kernelFunction_); 330 interscaleManager_.add_diagonal_matrix(accumulantInverseVolumes_, 331 "AccumulantInverseVolumes"); 332 interscaleManager_.add_diagonal_matrix(accumulantWeights_, 333 "AccumulantWeights"); 334 } 335 else { 336 // mesh-based accumulants 337 if (kernelOnTheFly_) { 338 ConstantQuantity<double> * atomCount = new ConstantQuantity<double>(this,1.); 339 interscaleManager_.add_per_atom_quantity(atomCount,"AtomCount"); 340 OnTheFlyMeshAccumulation * myWeights 341 = new OnTheFlyMeshAccumulation(this, 342 atomCount, atomCoarseGrainingPositions_); 343 interscaleManager_.add_dense_matrix(myWeights, 344 "KernelInverseWeights"); 345 accumulantWeights_ = new OnTheFlyKernelWeights(myWeights); 346 } else { 347 accumulant_ = shpFcn_; 348 accumulantWeights_ = new AccumulantWeights(accumulant_); 349 interscaleManager_.add_diagonal_matrix(accumulantWeights_, 350 "AccumulantWeights"); 351 } 352 } 353 // gradient matrix 354 if (gradFlags_.has_member(true)) { 355 NativeShapeFunctionGradient * gradientMatrix = new NativeShapeFunctionGradient(this); 356 interscaleManager_.add_vector_sparse_matrix(gradientMatrix,"GradientMatrix"); 357 gradientMatrix_ = gradientMatrix; 358 } 359 } 360 //------------------------------------------------------------------- construct_molecule_transfers()361 void ATC_Transfer::construct_molecule_transfers() 362 { 363 // molecule centroid, molecule charge, dipole moment and quadrupole moment calculations KKM add 364 if (!moleculeIds_.empty()) { 365 map<string,pair<MolSize,int> >::const_iterator molecule; 366 InterscaleManager & interscaleManager = this->interscale_manager(); // KKM add, may be we do not need this as interscaleManager_ already exists. 367 PerAtomQuantity<double> * atomProcGhostCoarseGrainingPositions_ = interscaleManager.per_atom_quantity("AtomicProcGhostCoarseGrainingPositions"); 368 FundamentalAtomQuantity * mass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,PROC_GHOST); 369 molecule = moleculeIds_.begin(); 370 int groupbit = (molecule->second).second; 371 smallMoleculeSet_ = new SmallMoleculeSet(this,groupbit); 372 smallMoleculeSet_->initialize(); // KKM add, why should we? 373 interscaleManager_.add_small_molecule_set(smallMoleculeSet_,"MoleculeSet"); 374 moleculeCentroid_ = new SmallMoleculeCentroid(this,mass,smallMoleculeSet_,atomProcGhostCoarseGrainingPositions_); 375 interscaleManager_.add_dense_matrix(moleculeCentroid_,"MoleculeCentroid"); 376 AtomToSmallMoleculeTransfer<double> * moleculeMass = 377 new AtomToSmallMoleculeTransfer<double>(this,mass,smallMoleculeSet_); 378 interscaleManager_.add_dense_matrix(moleculeMass,"MoleculeMass"); 379 FundamentalAtomQuantity * atomicCharge = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_CHARGE,PROC_GHOST); 380 AtomToSmallMoleculeTransfer<double> * moleculeCharge = 381 new AtomToSmallMoleculeTransfer<double>(this,atomicCharge,smallMoleculeSet_); 382 interscaleManager_.add_dense_matrix(moleculeCharge,"MoleculeCharge"); 383 dipoleMoment_ = new SmallMoleculeDipoleMoment(this,atomicCharge,smallMoleculeSet_,atomProcGhostCoarseGrainingPositions_,moleculeCentroid_); 384 interscaleManager_.add_dense_matrix(dipoleMoment_,"DipoleMoment"); 385 quadrupoleMoment_ = new SmallMoleculeQuadrupoleMoment(this,atomicCharge,smallMoleculeSet_,atomProcGhostCoarseGrainingPositions_,moleculeCentroid_); 386 interscaleManager_.add_dense_matrix(quadrupoleMoment_,"QuadrupoleMoment"); 387 } 388 } 389 //---------------------------------------------------------------------- 390 // constructs quantities construct_transfers()391 void ATC_Transfer::construct_transfers() 392 { 393 394 // set pointer to positions 395 // REFACTOR use method's handling of xref/xpointer 396 set_xPointer(); 397 398 ATC_Method::construct_transfers(); 399 400 // reference potential energy 401 if (setRefPE_) { 402 if (!setRefPEvalue_ && !readRefPE_) { 403 FieldManager fmgr(this); 404 nodalRefPotentialEnergy_ = fmgr.nodal_atomic_field(REFERENCE_POTENTIAL_ENERGY); 405 } 406 else { 407 nodalRefPotentialEnergy_ = new DENS_MAN(nNodes_,1); 408 nodalRefPotentialEnergy_->set_memory_type(PERSISTENT); 409 interscaleManager_.add_dense_matrix(nodalRefPotentialEnergy_, 410 field_to_string(REFERENCE_POTENTIAL_ENERGY)); 411 } 412 } 413 414 // for hardy-based fluxes 415 416 bool needsBondMatrix = (! bondOnTheFly_ ) && 417 (fieldFlags_(STRESS) 418 || fieldFlags_(ESHELBY_STRESS) 419 || fieldFlags_(HEAT_FLUX)); 420 if (needsBondMatrix) { 421 if (hasPairs_ && hasBonds_) { 422 pairMap_ = new PairMapBoth(lammpsInterface_,groupbit_); 423 } 424 else if (hasBonds_) { 425 pairMap_ = new PairMapBond(lammpsInterface_,groupbit_); 426 } 427 else if (hasPairs_) { 428 pairMap_ = new PairMapNeighbor(lammpsInterface_,groupbit_); 429 } 430 } 431 if (pairMap_) interscaleManager_.add_pair_map(pairMap_,"PairMap"); 432 433 if ( fieldFlags_(STRESS) || fieldFlags_(ESHELBY_STRESS) || fieldFlags_(HEAT_FLUX) ) { 434 435 const FE_Mesh * fe_mesh = feEngine_->fe_mesh(); 436 if (!kernelBased_) { 437 bondMatrix_ = new BondMatrixPartitionOfUnity(lammpsInterface_,*pairMap_,xPointer_,fe_mesh,accumulantInverseVolumes_); 438 } 439 else { 440 bondMatrix_ = new BondMatrixKernel(lammpsInterface_,*pairMap_,xPointer_,fe_mesh,kernelFunction_); 441 } 442 } 443 if (bondMatrix_) interscaleManager_.add_sparse_matrix(bondMatrix_,"BondMatrix"); 444 445 if ( fieldFlags_(STRESS) || fieldFlags_(ESHELBY_STRESS) ) { 446 if (atomToElementMapType_ == LAGRANGIAN) { 447 pairVirial_ = new PairVirialLagrangian(lammpsInterface_,*pairMap_,xref_); 448 } 449 else if (atomToElementMapType_ == EULERIAN) { 450 pairVirial_ = new PairVirialEulerian(lammpsInterface_,*pairMap_); 451 } 452 else { 453 throw ATC_Error("no atom to element map specified"); 454 } 455 } 456 if (pairVirial_) interscaleManager_.add_dense_matrix(pairVirial_,"PairVirial"); 457 458 if ( fieldFlags_(HEAT_FLUX) ) { 459 if (atomToElementMapType_ == LAGRANGIAN) { 460 pairHeatFlux_ = new PairPotentialHeatFluxLagrangian(lammpsInterface_,*pairMap_,xref_); 461 } 462 else if (atomToElementMapType_ == EULERIAN) { 463 pairHeatFlux_ = new PairPotentialHeatFluxEulerian(lammpsInterface_,*pairMap_); 464 } 465 else { 466 throw ATC_Error("no atom to element map specified"); 467 } 468 } 469 if (pairHeatFlux_) interscaleManager_.add_dense_matrix(pairHeatFlux_,"PairHeatFlux"); 470 471 FieldManager fmgr(this); 472 473 // for(int index=0; index < NUM_TOTAL_FIELDS; ++index) 474 for(int i=0; i < numFields_; ++i) { 475 FieldName index = indices_[i]; 476 if (fieldFlags_(index)) { 477 outputFields_[index] = fmgr.nodal_atomic_field(index); 478 } 479 } 480 481 // WIP REJ - move to fmgr 482 if (fieldFlags_(ELECTRIC_POTENTIAL)) { 483 PerAtomQuantity<double> * atomCharge = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_CHARGE); 484 restrictedCharge_ = fmgr.restricted_atom_quantity(CHARGE_DENSITY,"default",atomCharge); 485 } 486 487 // computes 488 map <string,int>::const_iterator iter; 489 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 490 string tag = iter->first; 491 ComputedAtomQuantity * c = new ComputedAtomQuantity(this, tag); 492 interscaleManager_.add_per_atom_quantity(c,tag); 493 int projection = iter->second; 494 DIAG_MAN * w = NULL; 495 if (projection == VOLUME_NORMALIZATION ) 496 { w = accumulantInverseVolumes_; } 497 else if (projection == NUMBER_NORMALIZATION ) 498 { w = accumulantWeights_; } 499 if (kernelFunction_ && kernelOnTheFly_) { 500 OnTheFlyKernelAccumulationNormalized * C = new OnTheFlyKernelAccumulationNormalized(this, c, kernelFunction_, atomCoarseGrainingPositions_, w); 501 interscaleManager_.add_dense_matrix(C,tag); 502 outputFieldsTagged_[tag] = C; 503 } 504 else { 505 AtfProjection * C = new AtfProjection(this, c, accumulant_, w); 506 interscaleManager_.add_dense_matrix(C,tag); 507 outputFieldsTagged_[tag] = C; 508 } 509 } 510 511 } 512 513 //------------------------------------------------------------------- 514 // sets initial values of filtered quantities construct_methods()515 void ATC_Transfer::construct_methods() 516 { 517 ATC_Method::construct_methods(); 518 519 if ((!initialized_) || timeFilterManager_.need_reset()) { 520 timeFilters_.reset(NUM_TOTAL_FIELDS+nComputes_); 521 sampleCounter_ = 0; 522 523 // for filtered fields 524 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 525 if (fieldFlags_(index)) { 526 string name = field_to_string((FieldName) index); 527 filteredData_[name] = 0.0; 528 timeFilters_(index) = timeFilterManager_.construct(); 529 } 530 } 531 532 // for filtered projected computes 533 534 // lists/accessing of fields ( & computes) 535 map <string,int>::const_iterator iter; 536 int index = NUM_TOTAL_FIELDS; 537 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 538 string tag = iter->first; 539 filteredData_[tag] = 0.0; 540 timeFilters_(index) = timeFilterManager_.construct(); 541 index++; 542 } 543 } 544 } 545 546 547 //------------------------------------------------------------------- 548 // called after the end of a "run" finish()549 void ATC_Transfer::finish() 550 { 551 // base class 552 ATC_Method::finish(); 553 } 554 555 //------------------------------------------------------------------- 556 // this is the parser modify(int narg,char ** arg)557 bool ATC_Transfer::modify(int narg, char **arg) 558 { 559 bool match = false; 560 561 int argIdx = 0; 562 // check to see if it is a transfer class command 563 /*! \page man_hardy_fields fix_modify AtC fields 564 \section syntax 565 fix_modify AtC fields <all | none> \n 566 fix_modify AtC fields <add | delete> <list_of_fields> \n 567 - all | none (keyword) = output all or no fields \n 568 - add | delete (keyword) = add or delete the listed output fields \n 569 - fields (keyword) = \n 570 density : mass per unit volume \n 571 displacement : displacement vector \n 572 momentum : momentum per unit volume \n 573 velocity : defined by momentum divided by density \n 574 projected_velocity : simple kernel estimation of atomic velocities \n 575 temperature : temperature derived from the relative atomic kinetic energy (as done by ) \n 576 kinetic_temperature : temperature derived from the full kinetic energy \n 577 number_density : simple kernel estimation of number of atoms per unit volume \n 578 stress : 579 Cauchy stress tensor for eulerian analysis (atom_element_map), or 580 1st Piola-Kirchhoff stress tensor for lagrangian analysis \n 581 transformed_stress : 582 1st Piola-Kirchhoff stress tensor for eulerian analysis (atom_element_map), or 583 Cauchy stress tensor for lagrangian analysis \n 584 heat_flux : spatial heat flux vector for eulerian, 585 or referential heat flux vector for lagrangian \n 586 potential_energy : potential energy per unit volume \n 587 kinetic_energy : kinetic energy per unit volume \n 588 thermal_energy : thermal energy (kinetic energy - continuum kinetic energy) per unit volume \n 589 internal_energy : total internal energy (potential + thermal) per unit volume \n 590 energy : total energy (potential + kinetic) per unit volume \n 591 number_density : number of atoms per unit volume \n 592 eshelby_stress: configurational stress (energy-momentum) tensor defined by Eshelby 593 [References: Philos. Trans. Royal Soc. London A, Math. Phys. Sci., Vol. 244, 594 No. 877 (1951) pp. 87-112; J. Elasticity, Vol. 5, Nos. 3-4 (1975) pp. 321-335] \n 595 vacancy_concentration: volume fraction of vacancy content \n 596 type_concentration: volume fraction of a specific atom type \n 597 \section examples 598 <TT> fix_modify AtC fields add velocity temperature </TT> 599 \section description 600 Allows modification of the fields calculated and output by the 601 transfer class. The commands are cumulative, e.g.\n 602 <TT> fix_modify AtC fields none </TT> \n 603 followed by \n 604 <TT> fix_modify AtC fields add velocity temperature </TT> \n 605 will only output the velocity and temperature fields. 606 \section restrictions 607 Must be used with the hardy/field type of AtC fix, see \ref man_fix_atc. 608 Currently, the stress and heat flux formulas are only correct for 609 central force potentials, e.g. Lennard-Jones and EAM 610 but not Stillinger-Weber. 611 \section related 612 See \ref man_hardy_gradients , \ref man_hardy_rates and \ref man_hardy_computes 613 \section default 614 By default, no fields are output 615 */ 616 if (strcmp(arg[argIdx],"fields")==0) { 617 argIdx++; 618 if (strcmp(arg[argIdx],"all")==0) { 619 outputFlags_ = true; 620 match = true; 621 } 622 else if (strcmp(arg[argIdx],"none")==0) { 623 outputFlags_ = false; 624 match = true; 625 } 626 else if (strcmp(arg[argIdx],"add")==0) { 627 argIdx++; 628 for (int i = argIdx; i < narg; ++i) { 629 FieldName field_name = string_to_field(arg[i]); 630 outputFlags_(field_name) = true; 631 } 632 match = true; 633 } 634 else if (strcmp(arg[argIdx],"delete")==0) { 635 argIdx++; 636 for (int i = argIdx; i < narg; ++i) { 637 FieldName field_name = string_to_field(arg[i]); 638 outputFlags_(field_name) = false; 639 } 640 match = true; 641 } 642 check_field_dependencies(); 643 if (fieldFlags_(DISPLACEMENT)) { trackDisplacement_ = true; } 644 } 645 646 /*! \page man_hardy_gradients fix_modify AtC gradients 647 \section syntax 648 fix_modify AtC gradients <add | delete> <list_of_fields> \n 649 - add | delete (keyword) = add or delete the calculation of gradients for the listed output fields \n 650 - fields (keyword) = \n 651 gradients can be calculated for all fields listed in \ref man_hardy_fields 652 653 \section examples 654 <TT> fix_modify AtC gradients add temperature velocity stress </TT> \n 655 <TT> fix_modify AtC gradients delete velocity </TT> \n 656 \section description 657 Requests calculation and ouput of gradients of the fields from the 658 transfer class. These gradients will be with regard to spatial or material 659 coordinate for eulerian or lagrangian analysis, respectively, as specified by 660 atom_element_map (see \ref man_atom_element_map ) 661 \section restrictions 662 Must be used with the hardy/field type of AtC fix 663 ( see \ref man_fix_atc ) 664 \section related 665 \section default 666 No gradients are calculated by default 667 */ 668 else if (strcmp(arg[argIdx],"gradients")==0) { 669 argIdx++; 670 if (strcmp(arg[argIdx],"add")==0) { 671 argIdx++; 672 FieldName field_name; 673 for (int i = argIdx; i < narg; ++i) { 674 field_name = string_to_field(arg[i]); 675 gradFlags_(field_name) = true; 676 } 677 match = true; 678 } 679 else if (strcmp(arg[argIdx],"delete")==0) { 680 argIdx++; 681 FieldName field_name; 682 for (int i = argIdx; i < narg; ++i) { 683 field_name = string_to_field(arg[i]); 684 gradFlags_(field_name) = false; 685 } 686 match = true; 687 } 688 } 689 690 /*! \page man_hardy_rates fix_modify AtC rates 691 \section syntax 692 fix_modify AtC rates <add | delete> <list_of_fields> \n 693 - add | delete (keyword) = add or delete the calculation of rates (time derivatives) for the listed output fields \n 694 - fields (keyword) = \n 695 rates can be calculated for all fields listed in \ref man_hardy_fields 696 697 \section examples 698 <TT> fix_modify AtC rates add temperature velocity stress </TT> \n 699 <TT> fix_modify AtC rates delete stress </TT> \n 700 \section description 701 Requests calculation and ouput of rates (time derivatives) of the fields from the 702 transfer class. For eulerian analysis (see \ref man_atom_element_map ), these rates 703 are the partial time derivatives of the nodal fields, not the full (material) time 704 derivatives. \n 705 \section restrictions 706 Must be used with the hardy/field type of AtC fix 707 ( see \ref man_fix_atc ) 708 \section related 709 \section default 710 No rates are calculated by default 711 */ 712 else if (strcmp(arg[argIdx],"rates")==0) { 713 argIdx++; 714 if (strcmp(arg[argIdx],"add")==0) { 715 argIdx++; 716 FieldName field_name; 717 for (int i = argIdx; i < narg; ++i) { 718 field_name = string_to_field(arg[i]); 719 rateFlags_(field_name) = true; 720 } 721 match = true; 722 } 723 else if (strcmp(arg[argIdx],"delete")==0) { 724 argIdx++; 725 FieldName field_name; 726 for (int i = argIdx; i < narg; ++i) { 727 field_name = string_to_field(arg[i]); 728 rateFlags_(field_name) = false; 729 } 730 match = true; 731 } 732 } 733 734 735 /*! \page man_pair_interactions fix_modify AtC pair_interactions/bond_interactions 736 \section syntax 737 fix_modify AtC pair_interactions <on|off> \n 738 fix_modify AtC bond_interactions <on|off> \n 739 740 \section examples 741 <TT> fix_modify AtC bond_interactions on </TT> \n 742 \section description 743 include bonds and/or pairs in the stress and heat flux computations 744 \section restrictions 745 \section related 746 \section default 747 pair interactions: on, bond interactions: off 748 */ 749 if (strcmp(arg[argIdx],"pair_interactions")==0) { // default true 750 argIdx++; 751 if (strcmp(arg[argIdx],"on")==0) { hasPairs_ = true; } 752 else { hasPairs_ = false;} 753 match = true; 754 } 755 if (strcmp(arg[argIdx],"bond_interactions")==0) { // default false 756 argIdx++; 757 if (strcmp(arg[argIdx],"on")==0) { hasBonds_ = true; } 758 else { hasBonds_ = false;} 759 match = true; 760 } 761 762 /*! \page man_hardy_computes fix_modify AtC computes 763 \section syntax 764 fix_modify AtC computes <add | delete> [per-atom compute id] <volume | number> \n 765 - add | delete (keyword) = add or delete the calculation of an equivalent continuum field 766 for the specified per-atom compute as volume or number density quantity \n 767 - per-atom compute id = name/id for per-atom compute, 768 fields can be calculated for all per-atom computes available from LAMMPS \n 769 - volume | number (keyword) = field created is a per-unit-volume quantity 770 or a per-atom quantity as weighted by kernel functions \n 771 772 \section examples 773 <TT> compute virial all stress/atom </TT> \n 774 <TT> fix_modify AtC computes add virial volume </TT> \n 775 <TT> fix_modify AtC computes delete virial </TT> \n 776 \n 777 <TT> compute centrosymmetry all centro/atom </TT> \n 778 <TT> fix_modify AtC computes add centrosymmetry number </TT> \n 779 \section description 780 Calculates continuum fields corresponding to specified per-atom computes created by LAMMPS \n 781 \section restrictions 782 Must be used with the hardy/field type of AtC fix ( see \ref man_fix_atc ) \n 783 Per-atom compute must be specified before corresponding continuum field can be requested \n 784 \section related 785 See manual page for compute 786 \section default 787 No defaults exist for this command 788 */ 789 else if (strcmp(arg[argIdx],"computes")==0) { 790 argIdx++; 791 if (strcmp(arg[argIdx],"add")==0) { 792 argIdx++; 793 string tag(arg[argIdx++]); 794 int normalization = NO_NORMALIZATION; 795 if (narg > argIdx) { 796 if (strcmp(arg[argIdx],"volume")==0) { 797 normalization = VOLUME_NORMALIZATION; 798 } 799 else if (strcmp(arg[argIdx],"number")==0) { 800 normalization = NUMBER_NORMALIZATION; 801 } 802 else if (strcmp(arg[argIdx],"mass")==0) { 803 normalization = MASS_NORMALIZATION; 804 throw ATC_Error("mass normalized not implemented"); 805 } 806 } 807 computes_[tag] = normalization; 808 nComputes_++; 809 match = true; 810 } 811 else if (strcmp(arg[argIdx],"delete")==0) { 812 argIdx++; 813 string tag(arg[argIdx]); 814 if (computes_.find(tag) != computes_.end()) { 815 computes_.erase(tag); 816 nComputes_--; 817 } 818 else { 819 throw ATC_Error(tag+" compute is not in list"); 820 } 821 match = true; 822 } 823 } 824 825 826 /*! \page man_sample_frequency fix_modify AtC sample_frequency 827 \section syntax 828 fix_modify AtC sample_frequency [freq] 829 - freq (int) : frequency to sample field in number of steps 830 \section examples 831 <TT> fix_modify AtC sample_frequency 10 832 \section description 833 Specifies a frequency at which fields are computed for the case 834 where time filters are being applied. 835 \section restrictions 836 Must be used with the hardy/field AtC fix ( see \ref man_fix_atc ) 837 and is only relevant when time filters are being used. 838 \section related 839 \section default 840 none 841 */ 842 else if (strcmp(arg[argIdx],"sample_frequency")==0) { 843 argIdx++; 844 int value = outputFrequency_; // default to output frequency 845 if (narg > 1) { 846 if (atoi(arg[argIdx]) > 0) value = atoi(arg[argIdx]); 847 } 848 sampleFrequency_ = value; 849 match = true; 850 } // end "sample_frequency" 851 852 // no match, call base class parser 853 if (!match) { 854 match = ATC_Method::modify(narg, arg); 855 } 856 857 return match; 858 } 859 860 //------------------------------------------------------------------- 861 // called at the beginning of a timestep pre_init_integrate()862 void ATC_Transfer::pre_init_integrate() 863 { 864 ATC_Method::pre_init_integrate(); 865 } 866 867 //------------------------------------------------------------------- 868 // called at the begining of second half timestep 869 // REFACTOR move this to post_neighbor pre_final_integrate()870 void ATC_Transfer::pre_final_integrate() 871 { 872 // update time 873 update_time(); // time uses step if dt = 0 874 875 876 877 if ( neighborReset_ && sample_now() ) { 878 if (! kernelOnTheFly_ ) { 879 if (!moleculeIds_.empty()) compute_kernel_matrix_molecule(); //KKM add 880 } 881 neighborReset_ = false; 882 } 883 } 884 885 //------------------------------------------------------------------- 886 // called at the end of second half timestep post_final_integrate()887 void ATC_Transfer::post_final_integrate() 888 { 889 // compute spatially smoothed quantities 890 double dt = lammpsInterface_->dt(); 891 if ( sample_now() ) { 892 893 bool needsBond = (! bondOnTheFly_ ) && 894 (fieldFlags_(STRESS) 895 || fieldFlags_(ESHELBY_STRESS) 896 || fieldFlags_(HEAT_FLUX)); 897 898 if ( needsBond ) { 899 if (pairMap_->need_reset()) { 900 // ATC::LammpsInterface::instance()->print_msg("Recomputing bond matrix due to atomReset_ value"); 901 compute_bond_matrix(); 902 } 903 } 904 time_filter_pre (dt); 905 compute_fields(); 906 time_filter_post(dt); 907 lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+sampleFrequency_); 908 } 909 910 // output 911 if ( output_now() && !outputStepZero_ ) output(); 912 outputStepZero_ = false; 913 914 //ATC_Method::post_final_integrate(); 915 916 } 917 918 //------------------------------------------------------------------- compute_bond_matrix(void)919 void ATC_Transfer::compute_bond_matrix(void) 920 { 921 bondMatrix_->reset(); 922 } 923 //------------------------------------------------------------------- compute_fields(void)924 void ATC_Transfer::compute_fields(void) 925 { 926 927 // keep per-atom computes fresh. JAZ and REJ not sure why; 928 // need to confer with JAT. (JAZ, 4/5/12) 929 interscaleManager_.lammps_force_reset(); 930 931 // (1) direct quantities 932 for(int i=0; i < numFields_; ++i) { 933 FieldName index = indices_[i]; 934 if (fieldFlags_(index)) { 935 DENS_MAT & data(hardyData_[field_to_string(index)].set_quantity()); 936 data = (outputFields_[index])->quantity(); 937 } 938 } 939 940 if (fieldFlags_(STRESS)) 941 compute_stress(hardyData_["stress"].set_quantity()); 942 if (fieldFlags_(HEAT_FLUX)) 943 compute_heatflux(hardyData_["heat_flux"].set_quantity()); 944 // molecule data 945 if (fieldFlags_(DIPOLE_MOMENT)) 946 compute_dipole_moment(hardyData_["dipole_moment"].set_quantity()); 947 if (fieldFlags_(QUADRUPOLE_MOMENT)) 948 compute_quadrupole_moment(hardyData_["quadrupole_moment"].set_quantity()); 949 if (fieldFlags_(DISLOCATION_DENSITY)) 950 compute_dislocation_density(hardyData_["dislocation_density"].set_quantity()); 951 952 // (2) derived quantities 953 // compute: gradients 954 if (gradFlags_.has_member(true)) { 955 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 956 if (gradFlags_(index)) { 957 string field= field_to_string((FieldName) index); 958 string grad_field = field + "_gradient"; 959 if (hardyData_.find(field) == hardyData_.end() ) { 960 throw ATC_Error("field " + field + " needs to be defined for gradient"); 961 } 962 gradient_compute(hardyData_[field].quantity(), hardyData_[grad_field].set_quantity()); 963 } 964 } 965 } 966 // compute: eshelby stress 967 if (fieldFlags_(ESHELBY_STRESS)) { 968 { 969 compute_eshelby_stress(hardyData_["eshelby_stress"].set_quantity(), 970 hardyData_["internal_energy"].quantity(), 971 hardyData_["stress"].quantity(), 972 hardyData_["displacement_gradient"].quantity()); 973 } 974 } 975 if (fieldFlags_(CAUCHY_BORN_ESHELBY_STRESS)) { 976 DENS_MAT & H = hardyData_["displacement_gradient"].set_quantity(); 977 DENS_MAT E(H.nRows(),1); 978 ATOMIC_DATA::const_iterator tfield = hardyData_.find("temperature"); 979 const DENS_MAT *temp = tfield==hardyData_.end() ? NULL : &((tfield->second).quantity()); 980 //DENS_MAT & T = hardyData_["temperature"]; 981 //cauchy_born_entropic_energy(H,E,T); E += hardyData_["internal_energy"]; 982 cauchy_born_energy(H, E, temp); 983 984 compute_eshelby_stress(hardyData_["cauchy_born_eshelby_stress"].set_quantity(), 985 E,hardyData_["stress"].quantity(), 986 hardyData_["displacement_gradient"].quantity()); 987 } 988 // compute: cauchy born stress 989 if (fieldFlags_(CAUCHY_BORN_STRESS)) { 990 ATOMIC_DATA::const_iterator tfield = hardyData_.find("temperature"); 991 const DENS_MAT *temp = tfield==hardyData_.end() ? NULL : &((tfield->second).quantity()); 992 cauchy_born_stress(hardyData_["displacement_gradient"].quantity(), 993 hardyData_["cauchy_born_stress"].set_quantity(), temp); 994 } 995 // compute: cauchy born energy 996 if (fieldFlags_(CAUCHY_BORN_ENERGY)) { 997 ATOMIC_DATA::const_iterator tfield = hardyData_.find("temperature"); 998 const DENS_MAT *temp = tfield==hardyData_.end() ? NULL : &((tfield->second).quantity()); 999 cauchy_born_energy(hardyData_["displacement_gradient"].quantity(), 1000 hardyData_["cauchy_born_energy"].set_quantity(), temp); 1001 } 1002 // 1st PK transformed to cauchy (lag) or cauchy transformed to 1st PK (eul) 1003 if (fieldFlags_(TRANSFORMED_STRESS)) { 1004 compute_transformed_stress(hardyData_["transformed_stress"].set_quantity(), 1005 hardyData_["stress"].quantity(), 1006 hardyData_["displacement_gradient"].quantity()); 1007 } 1008 if (fieldFlags_(VACANCY_CONCENTRATION)) { 1009 compute_vacancy_concentration(hardyData_["vacancy_concentration"].set_quantity(), 1010 hardyData_["displacement_gradient"].quantity(), 1011 hardyData_["number_density"].quantity()); 1012 } 1013 if (fieldFlags_(ELECTRIC_POTENTIAL)) { 1014 compute_electric_potential( 1015 hardyData_[field_to_string(ELECTRIC_POTENTIAL)].set_quantity()); 1016 } 1017 // compute: rotation and/or stretch from deformation gradient 1018 if (fieldFlags_(ROTATION) || fieldFlags_(STRETCH)) { 1019 compute_polar_decomposition(hardyData_["rotation"].set_quantity(), 1020 hardyData_["stretch"].set_quantity(), 1021 hardyData_["displacement_gradient"].quantity()); 1022 } 1023 // compute: rotation and/or stretch from deformation gradient 1024 if (fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT)) { 1025 compute_elastic_deformation_gradient2(hardyData_["elastic_deformation_gradient"].set_quantity(), 1026 hardyData_["stress"].quantity(), 1027 hardyData_["displacement_gradient"].quantity()); 1028 } 1029 1030 // (3) computes 1031 lammpsInterface_->computes_clearstep(); 1032 map <string,int>::const_iterator iter; 1033 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 1034 string tag = iter->first; 1035 COMPUTE_POINTER cmpt = lammpsInterface_->compute_pointer(tag); 1036 int projection = iter->second; 1037 int ncols = lammpsInterface_->compute_ncols_peratom(cmpt);; 1038 DENS_MAT atomicData(nLocal_,ncols); 1039 if (ncols == 1) { 1040 double * atomData = lammpsInterface_->compute_vector_peratom(cmpt); 1041 for (int i = 0; i < nLocal_; i++) { 1042 int atomIdx = internalToAtom_(i); 1043 atomicData(i,0) = atomData[atomIdx]; 1044 } 1045 } 1046 else { 1047 double ** atomData = lammpsInterface_->compute_array_peratom(cmpt); 1048 for (int i = 0; i < nLocal_; i++) { 1049 int atomIdx = internalToAtom_(i); 1050 for (int k = 0; k < ncols; k++) { 1051 atomicData(i,k) = atomData[atomIdx][k]; 1052 } 1053 } 1054 } 1055 // REFACTOR -- make dep manage 1056 if (projection == NO_NORMALIZATION) { 1057 project(atomicData,hardyData_[tag].set_quantity()); 1058 } 1059 else if (projection == VOLUME_NORMALIZATION) { 1060 project_volume_normalized(atomicData,hardyData_[tag].set_quantity()); 1061 } 1062 else if (projection == NUMBER_NORMALIZATION) { 1063 project_count_normalized(atomicData,hardyData_[tag].set_quantity()); 1064 } 1065 else if (projection == MASS_NORMALIZATION) { 1066 throw ATC_Error("unimplemented normalization"); 1067 } 1068 else { 1069 throw ATC_Error("unimplemented normalization"); 1070 } 1071 #ifdef ESHELBY_VIRIAL 1072 if (tag == "virial" && fieldFlags_(ESHELBY_STRESS)) { 1073 if (atomToElementMapType_ == LAGRANGIAN) { 1074 DENS_MAT tmp = hardyData_[tag].quantity(); 1075 DENS_MAT & myData(hardyData_[tag].set_quantity()); 1076 myData.reset(nNodes_,FieldSizes[STRESS]); 1077 DENS_MAT F(3,3),FT(3,3),FTINV(3,3),CAUCHY(3,3),PK1(3,3); 1078 const DENS_MAT& H(hardyData_["displacement_gradient"].quantity()); 1079 for (int k = 0; k < nNodes_; k++ ) { 1080 vector_to_symm_matrix(k,tmp,CAUCHY); 1081 vector_to_matrix(k,H,F); 1082 F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0; 1083 FT = F.transpose(); 1084 FTINV = inv(FT); 1085 1086 // volumes are already reference volumes. 1087 PK1 = CAUCHY*FTINV; 1088 matrix_to_vector(k,PK1,myData); 1089 } 1090 } 1091 compute_eshelby_stress(hardyData_["eshelby_virial"].set_quantity(), 1092 hardyData_["internal_energy"].quantity(),hardyData_[tag].quantity(), 1093 hardyData_["displacement_gradient"].quantity()); 1094 } 1095 #endif 1096 } 1097 1098 }// end of compute_fields routine 1099 1100 //------------------------------------------------------------------- time_filter_pre(double dt)1101 void ATC_Transfer::time_filter_pre(double dt) 1102 { 1103 sampleCounter_++; 1104 string name; 1105 double delta_t = dt*sampleFrequency_; 1106 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 1107 if (fieldFlags_(index)) { 1108 name = field_to_string((FieldName) index); 1109 timeFilters_(index)->apply_pre_step1(filteredData_[name].set_quantity(), 1110 hardyData_[name].quantity(), delta_t); 1111 } 1112 } 1113 map <string,int>::const_iterator iter; 1114 int index = NUM_TOTAL_FIELDS; 1115 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 1116 string tag = iter->first; 1117 timeFilters_(index)->apply_pre_step1(filteredData_[tag].set_quantity(), 1118 hardyData_[tag].quantity(), delta_t); 1119 index++; 1120 } 1121 } 1122 1123 //------------------------------------------------------------------- time_filter_post(double dt)1124 void ATC_Transfer::time_filter_post(double dt) 1125 { 1126 sampleCounter_++; 1127 string name; 1128 double delta_t = dt*sampleFrequency_; 1129 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 1130 if (fieldFlags_(index)) { 1131 name = field_to_string((FieldName) index); 1132 timeFilters_(index)->apply_post_step2(filteredData_[name].set_quantity(), 1133 hardyData_[name].quantity(), delta_t); 1134 } 1135 } 1136 if (rateFlags_.has_member(true)) { 1137 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 1138 if (rateFlags_(index)) { 1139 string field= field_to_string((FieldName) index); 1140 string rate_field = field + "_rate"; 1141 timeFilters_(index)->rate(hardyData_[rate_field].set_quantity(), 1142 filteredData_[field].quantity(), 1143 hardyData_[field].quantity(), delta_t); 1144 } 1145 } 1146 } 1147 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 1148 if (rateFlags_(index)) { 1149 name = field_to_string((FieldName) index); 1150 string rate_field = name + "_rate"; 1151 filteredData_[rate_field] = hardyData_[rate_field]; 1152 } 1153 if (gradFlags_(index)) { 1154 name = field_to_string((FieldName) index); 1155 string grad_field = name + "_gradient"; 1156 filteredData_[grad_field] = hardyData_[grad_field]; 1157 } 1158 } 1159 1160 // lists/accessing of fields ( & computes) 1161 map <string,int>::const_iterator iter; 1162 int index = NUM_TOTAL_FIELDS; 1163 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 1164 string tag = iter->first; 1165 timeFilters_(index)->apply_post_step2(filteredData_[tag].set_quantity(), 1166 hardyData_[tag].quantity(), delta_t); 1167 #ifdef ESHELBY_VIRIAL 1168 if (tag == "virial" && fieldFlags_(ESHELBY_STRESS)) { 1169 filteredData_["eshelby_virial"] = hardyData_["eshelby_virial"]; 1170 } 1171 #endif 1172 index++; 1173 } 1174 } 1175 1176 //------------------------------------------------------------------- output()1177 void ATC_Transfer::output() 1178 { 1179 feEngine_->departition_mesh(); 1180 1181 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 1182 if (outputFlags_(index)) { 1183 FieldName fName = (FieldName) index; 1184 string name= field_to_string(fName); 1185 fields_[fName] = filteredData_[name]; 1186 } 1187 } 1188 1189 ATC_Method::output(); 1190 if (lammpsInterface_->comm_rank() == 0) { 1191 // data 1192 OUTPUT_LIST output_data; 1193 #ifdef REFERENCE_PE_OUTPUT 1194 output_data["reference_potential_energy"] = nodalRefPotentialEnergy_->quantity(); 1195 #endif 1196 for(int index=0; index < NUM_TOTAL_FIELDS; ++index) { 1197 if (outputFlags_(index)) { 1198 string name= field_to_string((FieldName) index); 1199 output_data[name] = & ( filteredData_[name].set_quantity()); 1200 } 1201 if (rateFlags_(index)) { 1202 string name= field_to_string((FieldName) index); 1203 string rate_name = name + "_rate"; 1204 output_data[rate_name] = & ( filteredData_[rate_name].set_quantity()); 1205 } 1206 if (gradFlags_(index)) { 1207 string name= field_to_string((FieldName) index); 1208 string grad_name = name + "_gradient"; 1209 output_data[grad_name] = & ( filteredData_[grad_name].set_quantity()); 1210 } 1211 } 1212 1213 // lists/accessing of fields ( & computes) 1214 map <string,int>::const_iterator iter; 1215 for (iter = computes_.begin(); iter != computes_.end(); iter++) { 1216 string tag = iter->first; 1217 output_data[tag] = & ( filteredData_[tag].set_quantity()); 1218 #ifdef ESHELBY_VIRIAL 1219 if (tag == "virial" && fieldFlags_(ESHELBY_STRESS)) { 1220 output_data["eshelby_virial"] = & ( filteredData_["eshelby_virial"].set_quantity() ); 1221 } 1222 #endif 1223 } 1224 1225 DENS_MAT nodalInverseVolumes = CLON_VEC(accumulantInverseVolumes_->quantity()); 1226 output_data["NodalInverseVolumes"] = &nodalInverseVolumes; 1227 1228 // output 1229 feEngine_->write_data(output_index(), & output_data); 1230 } 1231 feEngine_->partition_mesh(); 1232 } 1233 1234 /////// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1235 1236 /////// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1237 //------------------------------------------------------------------- 1238 // computes nodeData = N*atomData project(const DENS_MAT & atomData,DENS_MAT & nodeData)1239 void ATC_Transfer::project(const DENS_MAT & atomData, 1240 DENS_MAT & nodeData) 1241 { 1242 if (! kernelOnTheFly_ ) { 1243 nodeData.reset(nNodes_,atomData.nCols(),true); 1244 DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols()); 1245 if (nLocal_>0) workNodeArray = (accumulant_->quantity()).transMat(atomData); 1246 int count = nodeData.nRows()*nodeData.nCols(); 1247 lammpsInterface_->allsum(workNodeArray.ptr(),nodeData.ptr(),count); 1248 } 1249 else { 1250 compute_projection(atomData,nodeData); 1251 } 1252 } 1253 1254 //------------------------------------------------------------------- 1255 // computes nodeData = N*molData specially for molecules project_molecule(const DENS_MAT & molData,DENS_MAT & nodeData)1256 void ATC_Transfer::project_molecule(const DENS_MAT & molData, 1257 DENS_MAT & nodeData) 1258 { 1259 if (! kernelOnTheFly_ ) { 1260 nodeData.reset(nNodes_,molData.nCols(),true); 1261 DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols()); 1262 if (nLocal_>0) workNodeArray = (accumulantMol_->quantity()).transMat(molData); 1263 int count = nodeData.nRows()*nodeData.nCols(); 1264 lammpsInterface_->allsum(workNodeArray.ptr(),nodeData.ptr(),count); 1265 } 1266 else { 1267 compute_projection(molData,nodeData); 1268 } 1269 } 1270 1271 //------------------------------------------------------------------- 1272 // computes nodeData = gradient of N*molData specially for molecules project_molecule_gradient(const DENS_MAT & molData,DENS_MAT & nodeData)1273 void ATC_Transfer::project_molecule_gradient(const DENS_MAT & molData, 1274 DENS_MAT & nodeData) 1275 { 1276 if (! kernelOnTheFly_ ) { 1277 nodeData.reset(nNodes_,molData.nCols(),true); 1278 DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols()); 1279 if (nLocal_>0) workNodeArray = (accumulantMolGrad_->quantity()).transMat(molData); 1280 int count = nodeData.nRows()*nodeData.nCols(); 1281 lammpsInterface_->allsum(workNodeArray.ptr(),nodeData.ptr(),count); 1282 } 1283 else { 1284 compute_projection(molData,nodeData); 1285 } 1286 } 1287 1288 //------------------------------------------------------------------- 1289 // count normalized project_count_normalized(const DENS_MAT & atomData,DENS_MAT & nodeData)1290 void ATC_Transfer::project_count_normalized(const DENS_MAT & atomData, 1291 DENS_MAT & nodeData) 1292 { 1293 DENS_MAT tmp; 1294 project(atomData,tmp); 1295 nodeData = (accumulantWeights_->quantity())*tmp; 1296 } 1297 1298 //------------------------------------------------------------------- 1299 // volume normalized project_volume_normalized(const DENS_MAT & atomData,DENS_MAT & nodeData)1300 void ATC_Transfer::project_volume_normalized(const DENS_MAT & atomData, 1301 DENS_MAT & nodeData) 1302 { 1303 DENS_MAT tmp; 1304 project(atomData,tmp); 1305 nodeData = (accumulantInverseVolumes_->quantity())*tmp; 1306 } 1307 1308 //------------------------------------------------------------------- 1309 // volume normalized molecule project_volume_normalized_molecule(const DENS_MAT & molData,DENS_MAT & nodeData)1310 void ATC_Transfer::project_volume_normalized_molecule(const DENS_MAT & molData, 1311 DENS_MAT & nodeData) 1312 { 1313 DENS_MAT tmp; 1314 project_molecule(molData,tmp); 1315 nodeData = (accumulantInverseVolumes_->quantity())*tmp; 1316 } 1317 1318 //------------------------------------------------------------------- 1319 // volume normalized molecule_gradient project_volume_normalized_molecule_gradient(const DENS_MAT & molData,DENS_MAT & nodeData)1320 void ATC_Transfer::project_volume_normalized_molecule_gradient(const DENS_MAT & molData, 1321 DENS_MAT & nodeData) 1322 { 1323 DENS_MAT tmp; 1324 project_molecule_gradient(molData,tmp); 1325 nodeData = (accumulantInverseVolumes_->quantity())*tmp; 1326 } 1327 1328 1329 //------------------------------------------------------------------- gradient_compute(const DENS_MAT & inNodeData,DENS_MAT & outNodeData)1330 void ATC_Transfer::gradient_compute(const DENS_MAT & inNodeData, 1331 DENS_MAT & outNodeData) 1332 { 1333 int nrows = inNodeData.nRows(); 1334 int ncols = inNodeData.nCols(); 1335 outNodeData.reset(nrows,ncols*nsd_); 1336 int index = 0; 1337 for (int n = 0; n < ncols; n++) { //output v1,1 v1,2 v1,3 ... 1338 for (int m = 0; m < nsd_; m++) { 1339 CLON_VEC inData(inNodeData,CLONE_COL,n); 1340 CLON_VEC outData(outNodeData,CLONE_COL,index); 1341 outData = (*((gradientMatrix_->quantity())[m]))*inData; 1342 ++index; 1343 } 1344 } 1345 } 1346 1347 1348 1349 //------------------------------------------------------------------- compute_force_matrix()1350 void ATC_Transfer::compute_force_matrix() 1351 { 1352 atomicForceMatrix_ = pairVirial_->quantity(); 1353 } 1354 1355 //------------------------------------------------------------------- 1356 // computes "virial" part of heat flux 1357 // This is correct ONLY for pair potentials. compute_heat_matrix()1358 void ATC_Transfer::compute_heat_matrix() 1359 { 1360 atomicHeatMatrix_ = pairHeatFlux_->quantity(); 1361 } 1362 1363 //------------------------------------------------------------------- 1364 // set xPointer_ to xref or xatom depending on Lagrangian/Eulerian analysis set_xPointer()1365 void ATC_Transfer::set_xPointer() 1366 { 1367 xPointer_ = xref_; 1368 if (atomToElementMapType_ == EULERIAN) { 1369 xPointer_ = lammpsInterface_->xatom(); 1370 } 1371 } 1372 1373 //------------------------------------------------------------------- 1374 // SOON TO BE OBSOLETE 1375 // check consistency of fieldFlags_ check_field_dependencies()1376 void ATC_Transfer::check_field_dependencies() 1377 { 1378 fieldFlags_ = outputFlags_; 1379 if (fieldFlags_(TRANSFORMED_STRESS)) { 1380 fieldFlags_(STRESS) = true; 1381 fieldFlags_(DISPLACEMENT) = true; 1382 } 1383 if (fieldFlags_(ESHELBY_STRESS)) { 1384 fieldFlags_(STRESS) = true; 1385 fieldFlags_(INTERNAL_ENERGY) = true; 1386 fieldFlags_(DISPLACEMENT) = true; 1387 } 1388 if (fieldFlags_(CAUCHY_BORN_STRESS) 1389 || fieldFlags_(CAUCHY_BORN_ENERGY) 1390 || fieldFlags_(CAUCHY_BORN_ESHELBY_STRESS) 1391 || fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT)) { 1392 if (! (cauchyBornStress_) ) { 1393 throw ATC_Error("can't compute cauchy-born stress w/o cauchy born model"); 1394 } 1395 } 1396 if (fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT)) { 1397 fieldFlags_(STRESS) = true; 1398 } 1399 if (fieldFlags_(CAUCHY_BORN_STRESS) 1400 || fieldFlags_(CAUCHY_BORN_ENERGY)) { 1401 fieldFlags_(TEMPERATURE) = true; 1402 fieldFlags_(DISPLACEMENT) = true; 1403 } 1404 if (fieldFlags_(CAUCHY_BORN_ESHELBY_STRESS)) { 1405 fieldFlags_(TEMPERATURE) = true; 1406 fieldFlags_(DISPLACEMENT) = true; 1407 fieldFlags_(STRESS) = true; 1408 } 1409 if (fieldFlags_(VACANCY_CONCENTRATION)) { 1410 fieldFlags_(DISPLACEMENT) = true; 1411 fieldFlags_(NUMBER_DENSITY) = true; 1412 } 1413 if (fieldFlags_(INTERNAL_ENERGY)) { 1414 fieldFlags_(POTENTIAL_ENERGY) = true; 1415 fieldFlags_(THERMAL_ENERGY) = true; 1416 } 1417 if (fieldFlags_(ENERGY)) { 1418 fieldFlags_(POTENTIAL_ENERGY) = true; 1419 fieldFlags_(KINETIC_ENERGY) = true; 1420 } 1421 if (fieldFlags_(TEMPERATURE) || fieldFlags_(HEAT_FLUX) || 1422 fieldFlags_(KINETIC_ENERGY) || fieldFlags_(THERMAL_ENERGY) || 1423 fieldFlags_(ENERGY) || fieldFlags_(INTERNAL_ENERGY) || 1424 fieldFlags_(KINETIC_ENERGY) || (fieldFlags_(STRESS) && 1425 atomToElementMapType_ == EULERIAN) ) { 1426 fieldFlags_(VELOCITY) = true; 1427 fieldFlags_(MASS_DENSITY) = true; 1428 } 1429 1430 if (fieldFlags_(VELOCITY)) { 1431 fieldFlags_(MASS_DENSITY) = true; 1432 fieldFlags_(MOMENTUM) = true; 1433 } 1434 if (fieldFlags_(DISPLACEMENT)) { 1435 fieldFlags_(MASS_DENSITY) = true; 1436 } 1437 if (fieldFlags_(TEMPERATURE) ) { 1438 fieldFlags_(NUMBER_DENSITY) = true; 1439 } 1440 1441 if (fieldFlags_(ROTATION) || 1442 fieldFlags_(STRETCH)) { 1443 fieldFlags_(DISPLACEMENT) = true; 1444 } 1445 if (fieldFlags_(ESHELBY_STRESS) 1446 || fieldFlags_(CAUCHY_BORN_STRESS) 1447 || fieldFlags_(CAUCHY_BORN_ENERGY) 1448 || fieldFlags_(CAUCHY_BORN_ESHELBY_STRESS) 1449 || fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT) 1450 || fieldFlags_(VACANCY_CONCENTRATION) 1451 || fieldFlags_(ROTATION) 1452 || fieldFlags_(STRETCH) ) { 1453 gradFlags_(DISPLACEMENT) = true; 1454 } 1455 1456 // check whether single_enable==0 for stress/heat flux calculation 1457 if (fieldFlags_(STRESS) || fieldFlags_(HEAT_FLUX)) { 1458 if (lammpsInterface_->single_enable()==0) { 1459 throw ATC_Error("Calculation of stress field not possible with selected pair type."); 1460 } 1461 } 1462 1463 } 1464 1465 //============== THIN WRAPPERS ==================================== 1466 // OBSOLETE 1467 // HARDY COMPUTES 1468 // ***************UNCONVERTED************************** 1469 1470 //------------------------------------------------------------------- 1471 // MOLECULE 1472 //------------------------------------------------------------------- compute_dipole_moment(DENS_MAT & dipole_moment)1473 void ATC_Transfer::compute_dipole_moment(DENS_MAT & dipole_moment) 1474 { 1475 const DENS_MAT & molecularVector(dipoleMoment_->quantity()); 1476 project_volume_normalized_molecule(molecularVector,dipole_moment); // KKM add 1477 // 1478 } 1479 //------------------------------------------------------------------- compute_quadrupole_moment(DENS_MAT & quadrupole_moment)1480 void ATC_Transfer::compute_quadrupole_moment(DENS_MAT & quadrupole_moment) 1481 { 1482 const DENS_MAT & molecularVector(quadrupoleMoment_->quantity()); 1483 project_volume_normalized_molecule_gradient(molecularVector,quadrupole_moment); // KKM add 1484 // 1485 } 1486 //------------------------------------------------------------------- compute_stress(DENS_MAT & stress)1487 void ATC_Transfer::compute_stress(DENS_MAT & stress) 1488 { 1489 // table of bond functions already calculated in initialize function 1490 // get conversion factor for nktV to p units 1491 double nktv2p = lammpsInterface_->nktv2p(); 1492 1493 // calculate kinetic energy tensor part of stress for Eulerian analysis 1494 if (atomToElementMapType_ == EULERIAN && nLocal_>0) { 1495 compute_kinetic_stress(stress); 1496 } 1497 else { 1498 // zero stress table for Lagrangian analysis or if nLocal_ = 0 1499 stress.zero(); 1500 } 1501 // add-in potential part of stress tensor 1502 int nrows = stress.nRows(); 1503 int ncols = stress.nCols(); 1504 DENS_MAT local_potential_hardy_stress(nrows,ncols); 1505 if (nLocal_>0) { 1506 if (bondOnTheFly_) { 1507 compute_potential_stress(local_potential_hardy_stress); 1508 } 1509 else { 1510 // compute table of force & position dyad 1511 compute_force_matrix(); 1512 // calculate force part of stress tensor 1513 local_potential_hardy_stress = atomicBondMatrix_*atomicForceMatrix_; 1514 local_potential_hardy_stress *= 0.5; 1515 } 1516 } 1517 // global summation of potential part of stress tensor 1518 DENS_MAT potential_hardy_stress(nrows,ncols); 1519 int count = nrows*ncols; 1520 lammpsInterface_->allsum(local_potential_hardy_stress.ptr(), 1521 potential_hardy_stress.ptr(), count); 1522 stress += potential_hardy_stress; 1523 stress = nktv2p*stress; 1524 } 1525 1526 //------------------------------------------------------------------- 1527 // kinetic energy portion of stress compute_kinetic_stress(DENS_MAT & stress)1528 void ATC_Transfer::compute_kinetic_stress(DENS_MAT& stress) 1529 { 1530 const DENS_MAT& density = hardyData_["mass_density"].quantity(); 1531 const DENS_MAT& velocity = hardyData_["velocity"].quantity(); 1532 1533 int * type = lammpsInterface_->atom_type(); 1534 double * mass = lammpsInterface_->atom_mass(); 1535 double * rmass = lammpsInterface_->atom_rmass(); 1536 double ** vatom = lammpsInterface_->vatom(); 1537 double mvv2e = lammpsInterface_->mvv2e(); // [MV^2]-->[Energy] 1538 1539 atomicTensor_.reset(nLocal_,6); 1540 for (int i = 0; i < nLocal_; i++) { 1541 int atomIdx = internalToAtom_(i); 1542 double ma = mass ? mass[type[atomIdx]]: rmass[atomIdx]; 1543 ma *= mvv2e; // convert mass to appropriate units 1544 double* v = vatom[atomIdx]; 1545 atomicTensor_(i,0) -= ma*v[0]*v[0]; 1546 atomicTensor_(i,1) -= ma*v[1]*v[1]; 1547 atomicTensor_(i,2) -= ma*v[2]*v[2]; 1548 atomicTensor_(i,3) -= ma*v[0]*v[1]; 1549 atomicTensor_(i,4) -= ma*v[0]*v[2]; 1550 atomicTensor_(i,5) -= ma*v[1]*v[2]; 1551 } 1552 project_volume_normalized(atomicTensor_, stress); 1553 1554 for (int i = 0; i < nNodes_; i++) { 1555 double rho_i = mvv2e*density(i,0); 1556 stress(i,0) += rho_i*velocity(i,0)*velocity(i,0); 1557 stress(i,1) += rho_i*velocity(i,1)*velocity(i,1); 1558 stress(i,2) += rho_i*velocity(i,2)*velocity(i,2); 1559 stress(i,3) += rho_i*velocity(i,0)*velocity(i,1); 1560 stress(i,4) += rho_i*velocity(i,0)*velocity(i,2); 1561 stress(i,5) += rho_i*velocity(i,1)*velocity(i,2); 1562 } 1563 } 1564 1565 //------------------------------------------------------------------- compute_heatflux(DENS_MAT & flux)1566 void ATC_Transfer::compute_heatflux(DENS_MAT & flux) 1567 { 1568 // calculate kinetic part of heat flux 1569 if (atomToElementMapType_ == EULERIAN && nLocal_>0) { 1570 compute_kinetic_heatflux(flux); 1571 } 1572 else { 1573 flux.zero(); // zero stress table for Lagrangian analysis 1574 } 1575 // add potential part of heat flux vector 1576 int nrows = flux.nRows(); 1577 int ncols = flux.nCols(); 1578 DENS_MAT local_hardy_heat(nrows,ncols); 1579 if (nLocal_>0) { 1580 if (bondOnTheFly_) { 1581 compute_potential_heatflux(local_hardy_heat); 1582 } 1583 else { 1584 // calculate force/potential-derivative part of heat flux 1585 compute_heat_matrix(); 1586 local_hardy_heat = atomicBondMatrix_*atomicHeatMatrix_; 1587 } 1588 } 1589 // global summation of potential part of heat flux vector 1590 DENS_MAT hardy_heat(nrows,ncols); 1591 int count = nrows*ncols; 1592 lammpsInterface_->allsum(local_hardy_heat.ptr(), 1593 hardy_heat.ptr(), count); 1594 flux += hardy_heat; 1595 } 1596 1597 //------------------------------------------------------------------- 1598 // compute kinetic part of heat flux compute_kinetic_heatflux(DENS_MAT & flux)1599 void ATC_Transfer::compute_kinetic_heatflux(DENS_MAT& flux) 1600 { 1601 const DENS_MAT& velocity = hardyData_["velocity"].quantity(); 1602 const DENS_MAT& energy = hardyData_["mass_density"].quantity(); 1603 const DENS_MAT& stress = hardyData_["stress"].quantity(); 1604 1605 int * type = lammpsInterface_->atom_type(); 1606 double * mass = lammpsInterface_->atom_mass(); 1607 double * rmass = lammpsInterface_->atom_rmass(); 1608 double ** vatom = lammpsInterface_->vatom(); 1609 double mvv2e = lammpsInterface_->mvv2e(); 1610 double * atomPE = lammpsInterface_->compute_pe_peratom(); 1611 double atomKE, atomEnergy; 1612 atomicVector_.reset(nLocal_,3); 1613 for (int i = 0; i < nLocal_; i++) { 1614 int atomIdx = internalToAtom_(i); 1615 double ma = mass ? mass[type[atomIdx]]: rmass[atomIdx]; 1616 ma *= mvv2e; // convert mass to appropriate units 1617 double* v = vatom[atomIdx]; 1618 atomKE = 0.0; 1619 for (int k = 0; k < nsd_; k++) { atomKE += v[k]*v[k]; } 1620 atomKE *= 0.5*ma; 1621 atomEnergy = atomKE + atomPE[atomIdx]; 1622 for (int j = 0; j < nsd_; j++) { 1623 atomicVector_(i,j) += atomEnergy*v[j]; 1624 } 1625 } 1626 project_volume_normalized(atomicVector_,flux); 1627 1628 // - e^0_I v_I + \sigma^T_I v_I 1629 for (int i = 0; i < nNodes_; i++) { 1630 double e_i = energy(i,0); 1631 flux(i,0) += (e_i + stress(i,0))*velocity(i,0) 1632 + stress(i,3)*velocity(i,1)+ stress(i,4)*velocity(i,2); 1633 flux(i,1) += (e_i + stress(i,1))*velocity(i,1) 1634 + stress(i,3)*velocity(i,0)+ stress(i,5)*velocity(i,2); 1635 flux(i,2) += (e_i + stress(i,2))*velocity(i,2) 1636 + stress(i,4)*velocity(i,0)+ stress(i,5)*velocity(i,1); 1637 } 1638 } 1639 //-------------------------------------------------------------------- compute_electric_potential(DENS_MAT & phi)1640 void ATC_Transfer::compute_electric_potential(DENS_MAT & phi) 1641 { 1642 // Poisson solve with insulating bcs 1643 const DENS_MAT & rho = (restrictedCharge_->quantity()); 1644 SPAR_MAT K; 1645 feEngine_->stiffness_matrix(K); 1646 double permittivity = lammpsInterface_->dielectric(); 1647 permittivity *= LammpsInterface::instance()->epsilon0(); 1648 K *= permittivity; 1649 BC_SET bcs; 1650 bcs.insert(pair<int,int>(0,0)); 1651 LinearSolver solver(K,bcs); 1652 CLON_VEC x = column(phi,0); 1653 CLON_VEC b = column(rho,0); 1654 solver.solve(x,b); 1655 //x.print("x:phi"); 1656 //b.print("b:rho"); 1657 //LinearSolver solver(K,AUTO_SOLVE,true); 1658 } 1659 //-------------------------------------------------------------------- compute_vacancy_concentration(DENS_MAT & Cv,const DENS_MAT & H,const DENS_MAT & rhoN)1660 void ATC_Transfer::compute_vacancy_concentration(DENS_MAT & Cv, 1661 const DENS_MAT & H, const DENS_MAT & rhoN) 1662 { 1663 int * type = lammpsInterface_->atom_type(); 1664 DENS_MAT new_rho(nNodes_,1); 1665 DENS_MAT atomCnt(nLocal_,1); 1666 double atomic_weight_sum = 0.0; 1667 double site_weight_sum = 0.0; 1668 int number_atoms = 0; 1669 const DIAG_MAT & myAtomicWeights(atomVolume_->quantity()); 1670 1671 for (int i = 0; i < nLocal_; i++) { 1672 int atomIdx = internalToAtom_(i); 1673 if (type[atomIdx] != 13) { 1674 atomCnt(i,0) = myAtomicWeights(i,i); 1675 atomic_weight_sum += myAtomicWeights(i,i); 1676 number_atoms++; 1677 } 1678 site_weight_sum += myAtomicWeights(i,i); 1679 } 1680 project_volume_normalized(atomCnt, new_rho); 1681 DENS_MAT F(3,3); 1682 for (int i = 0; i < nNodes_; i++) { 1683 if (atomToElementMapType_ == EULERIAN) { 1684 // for Eulerian analysis: F = (1-H)^{-1} 1685 DENS_MAT G(3,3); 1686 vector_to_matrix(i,H,G); 1687 G *= -1.; 1688 G(0,0) += 1.0; G(1,1) += 1.0; G(2,2) += 1.0; 1689 F = inv(G); 1690 } 1691 else if (atomToElementMapType_ == LAGRANGIAN) { 1692 // for Lagrangian analysis: F = (1+H) 1693 vector_to_matrix(i,H,F); 1694 F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0; 1695 } 1696 double J = det(F); 1697 double volume_per_atom = lammpsInterface_->volume_per_atom(); 1698 J *= volume_per_atom; 1699 Cv(i,0) = 1.0 - J*new_rho(i,0); 1700 } 1701 } 1702 1703 //-------------------------------------------------------------------- compute_eshelby_stress(DENS_MAT & M,const DENS_MAT & E,const DENS_MAT & S,const DENS_MAT & H)1704 void ATC_Transfer::compute_eshelby_stress(DENS_MAT & M, 1705 const DENS_MAT & E, const DENS_MAT & S, const DENS_MAT & H) 1706 { 1707 // eshelby stress:M, energy:E, stress:S, displacement gradient: H 1708 // eshelby stress = W I - F^T.P = W I - C.S [energy] 1709 // symmetric if isotropic S = a_0 I + a_1 C + a_2 C^2 1710 M.reset(nNodes_,FieldSizes[ESHELBY_STRESS]); 1711 double nktv2p = lammpsInterface_->nktv2p(); 1712 DENS_MAT P(3,3),F(3,3),FT(3,3),FTP(3,3),ESH(3,3); 1713 for (int i = 0; i < nNodes_; i++) { 1714 double W = E(i,0); 1715 ESH.identity(); 1716 ESH *= W; 1717 1718 // copy to local 1719 if (atomToElementMapType_ == LAGRANGIAN) { 1720 // Stress notation convention:: 0:11 1:12 2:13 3:21 4:22 5:23 6:31 7:32 8:33 1721 vector_to_matrix(i,S,P); 1722 1723 1724 vector_to_matrix(i,H,F); 1725 #ifndef H_BASED 1726 F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0; 1727 #endif 1728 FT = F.transpose(); 1729 } 1730 else if (atomToElementMapType_ == EULERIAN) { 1731 vector_to_symm_matrix(i,S,P); 1732 vector_to_matrix(i,H,F); 1733 FT = F.transpose(); 1734 } 1735 FTP = (1.0/nktv2p)*FT*P; 1736 ESH -= FTP; 1737 if (atomToElementMapType_ == EULERIAN) { 1738 // For Eulerian analysis, M = F^T*(w-H^T.CauchyStress) 1739 DENS_MAT Q(3,3); 1740 Q.identity(); 1741 // Q stores (1-H) 1742 Q -= FT.transpose(); 1743 DENS_MAT F(3,3); 1744 F = inv(Q); 1745 FT = F.transpose(); 1746 ESH = FT*ESH; 1747 } 1748 // copy to global 1749 matrix_to_vector(i,ESH,M); 1750 } 1751 } 1752 //--------------------------------------------------------------------------- 1753 // Computes the Cauchy Born stress tensor, T given displacement gradient, H 1754 // and optional temperature argument (passed by pointer), TEMP 1755 //--------------------------------------------------------------------------- cauchy_born_stress(const DENS_MAT & H,DENS_MAT & T,const DENS_MAT * temp)1756 void ATC_Transfer::cauchy_born_stress(const DENS_MAT &H, DENS_MAT &T, const DENS_MAT *temp) 1757 { 1758 FIELD_MATS uField; // uField should contain temperature. 1759 DENS_MAT_VEC tField; 1760 GRAD_FIELD_MATS hField; 1761 DENS_MAT_VEC &h = hField[DISPLACEMENT]; 1762 h.assign(nsd_, DENS_MAT(nNodes_,nsd_)); 1763 tField.assign(nsd_, DENS_MAT(nNodes_,nsd_)); 1764 // each row is the CB stress at a node stored in voigt form 1765 T.reset(nNodes_,FieldSizes[CAUCHY_BORN_STRESS]); 1766 const double nktv2p = lammpsInterface_->nktv2p(); 1767 const double fact = -lammpsInterface_->mvv2e()*nktv2p; 1768 1769 // reshape H (#nodes,9) into h [3](#nodes,3) displacement gradient 1770 vector_to_dens_mat_vec(H,h); 1771 1772 // if temperature is provided, then set it 1773 if (temp) uField[TEMPERATURE] = *temp; 1774 1775 // Computes the stress at each node. 1776 cauchyBornStress_->stress(uField, hField, tField); 1777 1778 // reshapes the stress, T to a (#node,6) DenseMatrix. 1779 DENS_MAT S(nNodes_,6); 1780 symm_dens_mat_vec_to_vector(tField,S); 1781 S *= fact; 1782 1783 // tField/S holds the 2nd P-K stress tensor. Transform to 1784 // Cauchy for EULERIAN analysis, transform to 1st P-K 1785 // for LAGRANGIAN analysis. 1786 DENS_MAT PK2(3,3),G(3,3),F(3,3),FT(3,3),STRESS(3,3); 1787 for (int i = 0; i < nNodes_; i++) { 1788 1789 vector_to_symm_matrix(i,S,PK2); 1790 1791 if (atomToElementMapType_ == EULERIAN) { 1792 1793 // for Eulerian analysis: F = (1-H)^{-1} 1794 vector_to_matrix(i,H,G); 1795 G *= -1.; 1796 1797 G(0,0) += 1.0; G(1,1) += 1.0; G(2,2) += 1.0; 1798 F = inv(G); 1799 FT = transpose(F); 1800 double J = det(F); 1801 STRESS = F*PK2*FT; 1802 STRESS *= 1/J; 1803 symm_matrix_to_vector(i,STRESS,T); 1804 } 1805 else{ 1806 // for Lagrangian analysis: F = 1 + H 1807 vector_to_matrix(i,H,F); 1808 1809 F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0; 1810 STRESS = F*PK2; 1811 matrix_to_vector(i,STRESS,T); 1812 } 1813 1814 } 1815 } 1816 //--------------------------------------------------------------------------- 1817 // Computes the Cauchy Born energy density, E given displacement gradient, H 1818 // and optional temperature argument (passed by pointer), TEMP 1819 //--------------------------------------------------------------------------- cauchy_born_energy(const DENS_MAT & H,DENS_MAT & E,const DENS_MAT * temp)1820 void ATC_Transfer::cauchy_born_energy(const DENS_MAT &H, DENS_MAT &E, const DENS_MAT *temp) 1821 { 1822 FIELD_MATS uField; // uField should contain temperature. 1823 GRAD_FIELD_MATS hField; 1824 DENS_MAT_VEC &h = hField[DISPLACEMENT]; 1825 h.assign(nsd_, DENS_MAT(nNodes_,nsd_)); 1826 1827 // reshape H (#nodes,9) into h [3](#nodes,3) displacement gradient 1828 vector_to_dens_mat_vec(H,h); 1829 1830 // if temperature is provided, then set it 1831 if (temp) uField[TEMPERATURE] = *temp; 1832 1833 // Computes the free/potential energy at each node. 1834 cauchyBornStress_->elastic_energy(uField, hField, E); 1835 1836 // convert back to energy units for ( ATC coupling uses MLT units) 1837 double mvv2e = lammpsInterface_->mvv2e(); // [MV^2]-->[Energy] 1838 E *= mvv2e; 1839 1840 // for Eulerian analysis, convert energy density to per-unit deformed volume 1841 if (atomToElementMapType_ == EULERIAN) { 1842 DENS_MAT G(3,3),F(3,3); 1843 for (int i = 0; i < nNodes_; i++) { 1844 // for Eulerian analysis: F = (1-H)^{-1} 1845 vector_to_matrix(i,H,G); 1846 G *= -1.; 1847 1848 G(0,0) += 1.0; G(1,1) += 1.0; G(2,2) += 1.0; 1849 F = inv(G); 1850 double J = det(F); 1851 E(i,0) *= 1/J; 1852 } 1853 } 1854 // subtract zero point energy 1855 if (nodalRefPotentialEnergy_) 1856 E -= nodalRefPotentialEnergy_->quantity(); 1857 } 1858 //--------------------------------------------------------------------------- 1859 // Computes the M/LH entropic energy density 1860 //--------------------------------------------------------------------------- cauchy_born_entropic_energy(const DENS_MAT & H,DENS_MAT & E,const DENS_MAT & T)1861 void ATC_Transfer::cauchy_born_entropic_energy(const DENS_MAT &H, DENS_MAT &E, const DENS_MAT &T) 1862 { 1863 FIELD_MATS uField; // uField should contain temperature. 1864 uField[TEMPERATURE] = T; 1865 GRAD_FIELD_MATS hField; 1866 DENS_MAT_VEC &h = hField[DISPLACEMENT]; 1867 h.assign(nsd_, DENS_MAT(nNodes_,nsd_)); 1868 1869 // reshape H (#nodes,9) into h [3](#nodes,3) displacement gradient 1870 vector_to_dens_mat_vec(H,h); 1871 1872 // Computes the free/potential energy at each node. 1873 cauchyBornStress_->entropic_energy(uField, hField, E); 1874 1875 // convert back to energy units for ( ATC coupling uses MLT units) 1876 double mvv2e = lammpsInterface_->mvv2e(); // [MV^2]-->[Energy] 1877 E *= mvv2e; 1878 1879 // for Eulerian analysis, convert energy density to per-unit deformed volume 1880 if (atomToElementMapType_ == EULERIAN) { 1881 DENS_MAT G(3,3),F(3,3); 1882 for (int i = 0; i < nNodes_; i++) { 1883 // for Eulerian analysis: F = (1-H)^{-1} 1884 vector_to_matrix(i,H,G); 1885 G *= -1.; 1886 1887 G(0,0) += 1.0; G(1,1) += 1.0; G(2,2) += 1.0; 1888 F = inv(G); 1889 double J = det(F); 1890 E(i,0) *= 1/J; 1891 } 1892 } 1893 1894 } 1895 //-------------------------------------------------------------------- compute_transformed_stress(DENS_MAT & stress,const DENS_MAT & T,const DENS_MAT & H)1896 void ATC_Transfer::compute_transformed_stress(DENS_MAT & stress, 1897 const DENS_MAT & T, const DENS_MAT & H) 1898 { 1899 stress.reset(nNodes_,FieldSizes[TRANSFORMED_STRESS]); 1900 DENS_MAT S(3,3),FT(3,3),P(3,3); 1901 for (int i = 0; i < nNodes_; i++) { 1902 if (atomToElementMapType_ == EULERIAN) { 1903 vector_to_symm_matrix(i,T,P); 1904 // for Eulerian analysis: F^T = (1-H)^{-T} 1905 DENS_MAT G(3,3); 1906 vector_to_matrix(i,H,G); 1907 G *= -1.; 1908 1909 G(0,0) += 1.0; G(1,1) += 1.0; G(2,2) += 1.0; 1910 FT = inv(G.transpose()); 1911 } 1912 else{ 1913 vector_to_matrix(i,T,P); 1914 // for Lagrangian analysis: F^T = (1+H)^T 1915 DENS_MAT F(3,3); 1916 vector_to_matrix(i,H,F); 1917 1918 F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0; 1919 FT = F.transpose(); 1920 } 1921 // 1922 double J = det(FT); 1923 FT *= 1/J; 1924 if (atomToElementMapType_ == EULERIAN) { 1925 FT = inv(FT); 1926 } 1927 S = P*FT; 1928 matrix_to_vector(i,S,stress); 1929 } 1930 } 1931 //-------------------------------------------------------------------- compute_polar_decomposition(DENS_MAT & rotation,DENS_MAT & stretch,const DENS_MAT & H)1932 void ATC_Transfer::compute_polar_decomposition(DENS_MAT & rotation, 1933 DENS_MAT & stretch, const DENS_MAT & H) 1934 { 1935 DENS_MAT F(3,3),R(3,3),U(3,3); 1936 for (int i = 0; i < nNodes_; i++) { 1937 vector_to_matrix(i,H,F); 1938 F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0; 1939 if (atomToElementMapType_ == EULERIAN) { 1940 polar_decomposition(F,R,U,false); } // F = V R 1941 else { 1942 polar_decomposition(F,R,U); } // F = R U 1943 // copy to local 1944 if ( fieldFlags_(ROTATION) ) { 1945 1946 matrix_to_vector(i,R,rotation); 1947 } 1948 if ( fieldFlags_(STRETCH) ) { 1949 matrix_to_vector(i,U,stretch); 1950 } 1951 } 1952 } 1953 //-------------------------------------------------------------------- compute_elastic_deformation_gradient(DENS_MAT & Fe,const DENS_MAT & P,const DENS_MAT & H)1954 void ATC_Transfer::compute_elastic_deformation_gradient(DENS_MAT & Fe, 1955 const DENS_MAT & P, const DENS_MAT & H) 1956 1957 { 1958 // calculate Fe for every node 1959 const double nktv2p = lammpsInterface_->nktv2p(); 1960 const double fact = 1.0/ ( lammpsInterface_->mvv2e()*nktv2p ); 1961 for (int i = 0; i < nNodes_; i++) { 1962 DENS_VEC Pv = global_vector_to_vector(i,P); 1963 Pv *= fact; 1964 CBElasticTangentOperator tangent(cauchyBornStress_, Pv); 1965 NonLinearSolver solver(&tangent); 1966 DENS_VEC Fv = global_vector_to_vector(i,H); // pass in initial guess 1967 add_identity_voigt_unsymmetric(Fv); 1968 solver.solve(Fv); 1969 vector_to_global_vector(i,Fv,Fe); 1970 } 1971 } 1972 //-------------------------------------------------------------------- compute_elastic_deformation_gradient2(DENS_MAT & Fe,const DENS_MAT & P,const DENS_MAT & H)1973 void ATC_Transfer::compute_elastic_deformation_gradient2(DENS_MAT & Fe, 1974 const DENS_MAT & P, const DENS_MAT & H) 1975 { 1976 // calculate Fe for every node 1977 const double nktv2p = lammpsInterface_->nktv2p(); 1978 const double fact = 1.0/ ( lammpsInterface_->mvv2e()*nktv2p ); 1979 DENS_MAT F(3,3),R(3,3),U(3,3),PP(3,3),S(3,3); 1980 for (int i = 0; i < nNodes_; i++) { 1981 // get F = RU 1982 vector_to_matrix(i,H,F); 1983 F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0; 1984 if (atomToElementMapType_ == EULERIAN) { 1985 polar_decomposition(F,R,U,false); } // F = V R 1986 else { 1987 polar_decomposition(F,R,U); } // F = R U 1988 // get S 1989 vector_to_matrix(i,P,PP); 1990 //S = PP*transpose(F); 1991 S = inv(F)*PP; 1992 1993 S += S.transpose(); S *= 0.5; // symmetrize 1994 DENS_VEC Sv = to_voigt(S); 1995 Sv *= fact; 1996 // solve min_U || S - S_CB(U) || 1997 CB2ndElasticTangentOperator tangent(cauchyBornStress_, Sv); 1998 NonLinearSolver solver(&tangent); 1999 //DENS_VEC Uv = to_voigt_unsymmetric(U); // pass in initial guess 2000 DENS_VEC Uv = to_voigt(U); // pass in initial guess 2001 //DENS_VEC Uv(6); Uv(0)=1;Uv(1)=1;Uv(2)=1;Uv(3)=0;Uv(4)=0;Uv(5)=0; 2002 solver.solve(Uv); 2003 DENS_MAT Ue = from_voigt(Uv); 2004 DENS_MAT FFe = R*Ue; 2005 matrix_to_vector(i,FFe,Fe); 2006 } 2007 } 2008 // =========== Analytical solutions ========================== 2009 } // end namespace ATC 2010 2011 2012