1 // ATC headers 2 #include "ATC_Method.h" 3 #include "LammpsInterface.h" 4 #include "FE_Engine.h" 5 #include "Array.h" 6 #include "Array2D.h" 7 #include "ATC_Error.h" 8 #include "Function.h" 9 #include "PrescribedDataManager.h" 10 #include "TimeIntegrator.h" 11 #include "PhysicsModel.h" 12 #include "PerAtomQuantityLibrary.h" 13 #include "TransferLibrary.h" 14 #include "KernelFunction.h" 15 #include "Utility.h" 16 #include "FieldManager.h" 17 18 #include <fstream> 19 #include <sstream> 20 #include <iostream> 21 22 using ATC_Utility::sgn; 23 using ATC_Utility::to_string; 24 using ATC_Utility::is_dbl; 25 26 using std::stringstream; 27 using std::ifstream; 28 using std::ofstream; 29 using std::string; 30 using std::map; 31 using std::set; 32 using std::vector; 33 using std::pair; 34 35 namespace ATC { 36 ATC_Method(string groupName,double ** & perAtomArray,LAMMPS_NS::Fix * thisFix)37 ATC_Method::ATC_Method(string groupName, double ** & perAtomArray, LAMMPS_NS::Fix * thisFix) : 38 nodalAtomicVolume_(NULL), 39 needReset_(true), 40 lammpsInterface_(LammpsInterface::instance()), 41 interscaleManager_(this), 42 timeFilterManager_(this), 43 integrateInternalAtoms_(false), 44 atomTimeIntegrator_(NULL), 45 ghostManager_(this), 46 feEngine_(NULL), 47 initialized_(false), 48 meshDataInitialized_(false), 49 localStep_(0), 50 sizeComm_(8), // 3 positions + 1 material id * 2 for output 51 atomCoarseGrainingPositions_(NULL), 52 atomGhostCoarseGrainingPositions_(NULL), 53 atomProcGhostCoarseGrainingPositions_(NULL), 54 atomReferencePositions_(NULL), 55 nsd_(lammpsInterface_->dimension()), 56 xref_(NULL), 57 readXref_(false), 58 needXrefProcessorGhosts_(false), 59 trackDisplacement_(false), 60 needsAtomToElementMap_(true), 61 atomElement_(NULL), 62 atomGhostElement_(NULL), 63 internalElementSet_(""), 64 atomMasses_(NULL), 65 atomPositions_(NULL), 66 atomVelocities_(NULL), 67 atomForces_(NULL), 68 parallelConsistency_(true), 69 outputNow_(false), 70 outputTime_(true), 71 outputFrequency_(0), 72 sampleFrequency_(0), 73 sampleCounter_(0), 74 peScale_(1./(lammpsInterface_->mvv2e())), 75 keScale_(1.), 76 scalarFlag_(0), 77 vectorFlag_(0), 78 sizeVector_(0), 79 scalarVectorFreq_(0), 80 sizePerAtomCols_(4), 81 perAtomOutput_(NULL), 82 perAtomArray_(perAtomArray), 83 extScalar_(0), 84 extVector_(0), 85 extList_(NULL), 86 thermoEnergyFlag_(0), 87 atomVolume_(NULL), 88 atomicWeightsWriteFlag_(false), 89 atomicWeightsWriteFrequency_(0), 90 atomWeightType_(LATTICE), 91 domainDecomposition_(REPLICATED_MEMORY), 92 groupbit_(0), 93 groupbitGhost_(0), 94 needProcGhost_(false), 95 groupTag_(groupName), 96 nLocal_(0), 97 nLocalTotal_(0), 98 nLocalGhost_(0), 99 atomToElementMapType_(LAGRANGIAN), 100 atomToElementMapFrequency_(0), 101 regionID_(-1), 102 mdMassNormalization_(false), 103 kernelBased_(false), 104 kernelOnTheFly_(false), 105 kernelFunction_(NULL), 106 bondOnTheFly_(false), 107 accumulant_(NULL), 108 accumulantMol_(NULL), 109 accumulantMolGrad_(NULL), 110 accumulantWeights_(NULL), 111 accumulantInverseVolumes_(&invNodeVolumes_), 112 accumulantBandwidth_(0), 113 useRestart_(false), 114 hasRefPE_(false), 115 setRefPE_(false), 116 setRefPEvalue_(false), 117 refPEvalue_(0.), 118 readRefPE_(false), 119 nodalRefPotentialEnergy_(NULL), 120 simTime_(0.0), 121 stepCounter_(0) 122 { 123 lammpsInterface_->print_msg_once("version "+version()); 124 lammpsInterface_->set_fix_pointer(thisFix); 125 interscaleManager_.set_lammps_data_prefix(); 126 grow_arrays(lammpsInterface_->nmax()); 127 feEngine_ = new FE_Engine(lammpsInterface_->world()); 128 129 130 lammpsInterface_->create_compute_pe_peratom(); 131 } 132 ~ATC_Method()133 ATC_Method::~ATC_Method() 134 { 135 lammpsInterface_->destroy_2d_double_array(xref_); 136 lammpsInterface_->destroy_2d_double_array(perAtomOutput_); 137 if (atomTimeIntegrator_) delete atomTimeIntegrator_; 138 if (feEngine_) delete feEngine_; 139 } 140 141 //-------------------------------------------------- 142 // pack_fields 143 // bundle all allocated field matrices into a list 144 // for output needs 145 //-------------------------------------------------- pack_fields(RESTART_LIST & data)146 void ATC_Method::pack_fields(RESTART_LIST & data) 147 { 148 map<FieldName,int>::const_iterator field; 149 for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) { 150 FieldName thisField = field->first; 151 string fieldName = field_to_string(thisField); 152 string matrixName; 153 // copy all fields from ATC_Method.h 154 matrixName = "fields_" + fieldName; 155 data[matrixName] = & fields_[thisField].set_quantity(); 156 matrixName = "dot_fields_" + fieldName; 157 data[matrixName] = & dot_fields_[thisField].set_quantity(); 158 matrixName = "ddot_fields_" + fieldName; 159 data[matrixName] = & ddot_fields_[thisField].set_quantity(); 160 matrixName = "dddot_fields_" + fieldName; 161 data[matrixName] = & dddot_fields_[thisField].set_quantity(); 162 matrixName = "NodalAtomicFieldsRoc_" + fieldName; 163 data[matrixName] = & nodalAtomicFieldsRoc_[thisField].set_quantity(); 164 } 165 } 166 167 //-------------------------------------------------- 168 // write_restart_file 169 // bundle matrices that need to be saved and call 170 // fe_engine to write the file 171 //-------------------------------------------------- write_restart_data(string fileName,RESTART_LIST & data)172 void ATC_Method::write_restart_data(string fileName, RESTART_LIST & data) 173 { 174 pack_fields(data); 175 feEngine_->write_restart_file(fileName,&data); 176 } 177 178 //-------------------------------------------------- 179 // read_restart_file 180 // bundle matrices that need to be saved and call 181 // fe_engine to write the file 182 //-------------------------------------------------- read_restart_data(string fileName,RESTART_LIST & data)183 void ATC_Method::read_restart_data(string fileName, RESTART_LIST & data) 184 { 185 pack_fields(data); 186 feEngine_->read_restart_file(fileName,&data); 187 } 188 189 //-------------------------------------------------- 190 // Interactions with LAMMPS fix commands 191 // parse input command and pass on to finite element engine 192 // or physics specific transfers if necessary 193 // revert to physics-specific transfer if no command matches input 194 // first keyword is unique to particular class 195 // base class keyword matching must apply to ALL physics 196 // order: derived, base, owned objects 197 //-------------------------------------------------- modify(int narg,char ** arg)198 bool ATC_Method::modify(int narg, char **arg) 199 { 200 int argIdx=0; 201 bool match = false; 202 203 // gateways to other modules e.g. extrinsic, control, mesh 204 // pass off to fe engine 205 if (strcmp(arg[argIdx],"mesh")==0) { 206 match = feEngine_->modify(narg, arg); 207 if (feEngine_->has_mesh() && !meshDataInitialized_) 208 this->initialize_mesh_data(); 209 } 210 // pass off to time filter 211 else if (strcmp(arg[argIdx],"filter")==0) { 212 argIdx++; 213 match = timeFilterManager_.modify(narg-argIdx,&arg[argIdx]); 214 215 // consistentInitialization_ = false; 216 } 217 // pass off to kernel function manager 218 else if (strcmp(arg[argIdx],"kernel")==0) { 219 argIdx++; 220 221 if (kernelFunction_) { 222 //delete kernelFunction_; 223 //resetKernelFunction_ = true; 224 } 225 kernelFunction_ = KernelFunctionMgr::instance()->function(&arg[argIdx],narg-argIdx); 226 if (kernelFunction_) match = true; 227 else ATC_Error("no matching kernel found"); 228 feEngine_->set_kernel(kernelFunction_); 229 230 accumulantMol_=&kernelAccumulantMol_; // KKM add 231 accumulantMolGrad_=&kernelAccumulantMolGrad_; // KKM add 232 } 233 // pass off to ghost manager 234 else if (strcmp(arg[argIdx],"boundary_dynamics")==0) { 235 argIdx++; 236 match = ghostManager_.modify(narg-argIdx,&arg[argIdx]); 237 } 238 239 // parsing handled here 240 else { 241 if (strcmp(arg[argIdx],"parallel_consistency")==0) { 242 argIdx++; 243 //if (!kernelFunction_) { throw ATC_Error("on_the_fly requires a kernel function"); } 244 if (strcmp(arg[argIdx],"off")==0) parallelConsistency_ = false; 245 else parallelConsistency_ = true; 246 match = true; 247 } 248 /*! \page man_hardy_on_the_fly fix_modify AtC on_the_fly 249 \section syntax 250 fix_modify AtC on_the_fly <bond | kernel> <optional on | off> \n - bond | kernel (keyword) = specifies on-the-fly calculation of bond or 251 kernel matrix elements \n 252 - on | off (keyword) = activate or discontinue on-the-fly mode \n 253 \section examples 254 <TT> fix_modify AtC on_the_fly bond on </TT> \n <TT> fix_modify AtC on_the_fly kernel </TT> \n 255 <TT> fix_modify AtC on_the_fly kernel off </TT> \n 256 \section description 257 Overrides normal mode of pre-calculating and storing bond pair-to-node a 258 nd 259 kernel atom-to-node matrices. If activated, will calculate elements of t 260 hese 261 matrices during repeated calls of field computations (i.e. "on-the-fly") and not store them for 262 future use. \n on flag is optional - if omitted, on_the_fly will be activated for the s 263 pecified 264 matrix. Can be deactivated using off flag. \n 265 \section restrictions 266 Must be used with the hardy/field type of AtC fix 267 ( see \ref man_fix_atc ) 268 \section related 269 \section default 270 By default, on-the-fly calculation is not active (i.e. off). However, code does a memory allocation check to determine if it can store all needed bond and kernel matrix ele ments. If this allocation fails, on-the-fly is activated. \n 271 */ 272 273 else if (strcmp(arg[argIdx],"on_the_fly")==0) { 274 argIdx++; 275 //if (!kernelFunction_) { throw ATC_Error("on_the_fly requires a kernel function"); } 276 if (strcmp(arg[argIdx],"bond")==0) { 277 argIdx++; 278 bondOnTheFly_ = true; 279 if (narg > argIdx && strcmp(arg[argIdx],"off")==0) bondOnTheFly_ = false; 280 } 281 else if (strcmp(arg[argIdx],"kernel")==0) { 282 argIdx++; 283 kernelOnTheFly_ = true; 284 if (narg > argIdx && strcmp(arg[argIdx],"off")==0) kernelOnTheFly_ = false; 285 } 286 else { throw ATC_Error("unsupported on_the_fly type"); } 287 match = true; 288 } 289 290 /*! \page man_output fix_modify AtC output 291 \section syntax 292 fix_modify AtC output <filename_prefix> <frequency> 293 [text | full_text | binary | vector_components | tensor_components ] 294 fix_modify AtC output index [step | time ] 295 - filename_prefix (string) = prefix for data files 296 - frequency (integer) = frequency of output in time-steps 297 - options (keyword/s): \n 298 text = creates text output of index, step and nodal variable values for unique nodes \n 299 full_text = creates text output index, nodal id, step, nodal coordinates and nodal variable values for unique and image nodes \n 300 binary = creates binary Ensight output \n 301 vector_components = outputs vectors as scalar components \n 302 tensor_components = outputs tensor as scalar components 303 (use this for Paraview)\n 304 305 \section examples 306 <TT> fix_modify AtC output heatFE 100 </TT> \n 307 <TT> fix_modify AtC output hardyFE 1 text tensor_components </TT> \n 308 <TT> fix_modify AtC output hardyFE 10 text binary tensor_components </TT> \n 309 <TT> fix_modify AtC output index step </TT> \n 310 \section description 311 Creates text and/or binary (Ensight, "gold" format) output of nodal/mesh data 312 which is transfer/physics specific. Output indexed by step or time is possible. 313 \section restrictions 314 \section related 315 see \ref man_fix_atc 316 \section default 317 no default format 318 output indexed by time 319 */ 320 else if (strcmp(arg[argIdx],"output")==0) { 321 argIdx++; 322 /*! \page man_output_nodeset fix_modify AtC output nodeset 323 \section syntax 324 fix_modify AtC output nodeset <nodeset_name> <operation> 325 - nodeset_name (string) = name of nodeset to be operated on 326 - operation (keyword/s): \n 327 sum = creates nodal sum over nodes in specified nodeset \n 328 \section examples 329 <TT> fix_modify AtC output nodeset nset1 sum </TT> \n 330 \section description 331 Performs operation over the nodes belonging to specified nodeset 332 and outputs resulting variable values to GLOBALS file. 333 \section restrictions 334 \section related 335 see \ref man_fix_atc 336 \section default 337 none 338 */ 339 if (strcmp(arg[argIdx],"nodeset")==0) { 340 argIdx++; 341 string nset = arg[argIdx++]; 342 if (strcmp(arg[argIdx],"sum")==0) { 343 argIdx++; 344 string field = arg[argIdx]; 345 pair < string, FieldName > id(nset,string_to_field(field)); 346 nsetData_[id] = NODESET_SUM; 347 match = true; 348 } 349 else if (strcmp(arg[argIdx],"average")==0) { 350 argIdx++; 351 string field = arg[argIdx]; 352 pair < string, FieldName > id(nset,string_to_field(field)); 353 nsetData_[id] = NODESET_AVERAGE; 354 match = true; 355 } 356 } 357 358 359 /*! \page man_boundary_integral fix_modify AtC output boundary_integral 360 \section syntax 361 fix_modify AtC output boundary_integral [field] faceset [name] 362 - field (string) : name of hardy field 363 - name (string) : name of faceset 364 \section examples 365 <TT> fix_modify AtC output boundary_integral stress faceset loop1 </TT> \n 366 \section description 367 Calculates a surface integral of the given field dotted with the 368 outward normal of the faces and puts output in the "GLOBALS" file 369 \section restrictions 370 Must be used with the hardy/field type of AtC fix 371 ( see \ref man_fix_atc ) 372 \section related 373 \section default 374 none 375 */ 376 377 /*! \page man_contour_integral fix_modify AtC output contour_integral 378 \section syntax 379 fix_modify AtC output contour_integral [field] faceset [name] <axis [x | y | z 380 ]> 381 - field (string) : name of hardy field 382 - name (string) : name of faceset 383 - axis (string) : x or y or z 384 \section examples 385 <TT> fix_modify AtC output contour_integral stress faceset loop1 </TT> \n 386 \section description 387 Calculates a surface integral of the given field dotted with the 388 outward normal of the faces and puts output in the "GLOBALS" file 389 \section restrictions 390 Must be used with the hardy/field type of AtC fix 391 ( see \ref man_fix_atc ) 392 \section related 393 \section default 394 none 395 */ 396 397 else if ( (strcmp(arg[argIdx],"boundary_integral")==0) 398 || (strcmp(arg[argIdx],"contour_integral")==0) ) { 399 FacesetIntegralType type = BOUNDARY_INTEGRAL; 400 if (strcmp(arg[argIdx],"contour_integral")==0) 401 type = CONTOUR_INTEGRAL; 402 argIdx++; 403 string field(arg[argIdx++]); 404 if(strcmp(arg[argIdx],"faceset")==0) { 405 argIdx++; 406 string name(arg[argIdx++]); 407 pair <string,string> pair_name(name,field); 408 fsetData_[pair_name] = type; 409 match = true; 410 } 411 } // end "boundary_integral" || "contour_integral" 412 413 /*! \page man_output_elementset fix_modify AtC output elementset 414 \section syntax 415 fix_modify AtC output volume_integral <eset_name> <field> {` 416 - set_name (string) = name of elementset to be integrated over 417 - fieldname (string) = name of field to integrate 418 csum = creates nodal sum over nodes in specified nodeset \n 419 \section examples 420 <TT> fix_modify AtC output eset1 mass_density </TT> \n 421 \section description 422 Performs volume integration of specified field over elementset 423 and outputs resulting variable values to GLOBALS file. 424 \section restrictions 425 \section related 426 see \ref man_fix_atc 427 \section default 428 none 429 */ 430 431 else if ( (strcmp(arg[argIdx],"volume_integral")==0) ) { 432 argIdx++; 433 string name(arg[argIdx++]); 434 string field(arg[argIdx++]); 435 pair <string,FieldName> pair_name(name,string_to_field(field)); 436 if (++argIdx < narg) { // keyword average 437 esetData_[pair_name] = ELEMENTSET_AVERAGE; 438 } 439 else { 440 esetData_[pair_name] = ELEMENTSET_TOTAL; 441 } 442 match = true; 443 } 444 445 else if (strcmp(arg[argIdx],"now")==0) { 446 argIdx++; 447 double dt = 1.0; 448 if (argIdx < narg) { 449 dt = atof(arg[argIdx++]); 450 } 451 update_time(dt); 452 update_step(); 453 outputNow_ = true; 454 this->output(); 455 outputNow_ = false; 456 match = true; 457 } 458 else 459 if (strcmp(arg[argIdx],"index")==0) { 460 argIdx++; 461 if (strcmp(arg[argIdx],"step")==0) { outputTime_ = false; } 462 else { outputTime_ = true; } 463 match = true; 464 } 465 else { 466 outputPrefix_ = arg[argIdx++]; 467 outputFrequency_ = atoi(arg[argIdx++]); 468 bool ensight_output = false, full_text_output = false; 469 bool text_output = false, vect_comp = false, tensor_comp = false; 470 int rank = lammpsInterface_->comm_rank(); 471 for (int i = argIdx; i<narg; ++i) { 472 if (strcmp(arg[i],"full_text")==0) full_text_output = true; 473 else if (strcmp(arg[i],"text")==0) text_output = true; 474 else if (strcmp(arg[i],"binary")==0) ensight_output = true; 475 else if (strcmp(arg[i],"vector_components")==0) vect_comp = true; 476 else if (strcmp(arg[i],"tensor_components")==0) tensor_comp = true; 477 else { throw ATC_Error(" output: unknown keyword "); } 478 } 479 if (outputFrequency_>0) { 480 set<int> otypes; 481 if (full_text_output || text_output) { 482 lammpsInterface_->print_msg_once("Warning : text output can create _LARGE_ files"); 483 } 484 if (full_text_output) otypes.insert(FULL_GNUPLOT); 485 if (text_output) otypes.insert(GNUPLOT); 486 if (ensight_output) otypes.insert(ENSIGHT); 487 if (ntracked() > 0) { 488 string fstem = field_to_string(SPECIES_CONCENTRATION); 489 string istem = field_to_intrinsic_name(SPECIES_CONCENTRATION); 490 vector<string> tnames = tracked_names(); 491 vector<string> fnames; 492 vector<string> inames; 493 for (unsigned int i = 0; i < tnames.size(); i++) { 494 fnames.push_back(fstem+tnames[i]); 495 inames.push_back(istem+tnames[i]); 496 } 497 feEngine_->add_field_names(fstem,fnames); 498 feEngine_->add_field_names(istem,inames); 499 } 500 feEngine_->initialize_output(rank,outputPrefix_,otypes); 501 if (vect_comp) 502 feEngine_->output_manager() 503 ->set_option(OUTPUT_VECTOR_COMPONENTS,true); 504 if (tensor_comp) 505 feEngine_->output_manager() 506 ->set_option(OUTPUT_TENSOR_COMPONENTS,true); 507 } 508 match = true; 509 } 510 } 511 // add a species for tracking 512 /*! \page man_add_species fix_modify AtC add_species 513 \section syntax 514 fix_modify_AtC add_species <TAG> <group|type> <ID> \n 515 - <TAG> = tag for tracking a species \n 516 - group|type = LAMMPS defined group or type of atoms \n 517 - <ID> = name of group or type number \n 518 \section examples 519 <TT> fix_modify AtC add_species gold type 1 </TT> \n 520 <TT> group GOLDGROUP type 1 </TT> \n 521 <TT> fix_modify AtC add_species gold group GOLDGROUP </TT> 522 \section description 523 Associates a tag with all atoms of a specified type or within a specified group. \n 524 \section restrictions 525 \section related 526 \section default 527 No defaults for this command. 528 */ 529 else if (strcmp(arg[argIdx],"add_species")==0) { 530 argIdx++; 531 string speciesTag = arg[argIdx]; 532 string tag = arg[argIdx]; 533 argIdx++; 534 if (strcmp(arg[argIdx],"group")==0) { 535 if (narg-argIdx == 2) { 536 string name = arg[++argIdx]; 537 int id = lammpsInterface_->group_bit(name); 538 groupList_.push_back(id); 539 groupNames_.push_back(tag); 540 } 541 else { 542 while (++argIdx < narg) { 543 string name = arg[argIdx]; 544 int id = lammpsInterface_->group_bit(name); 545 string tag = speciesTag+"-"+name; 546 groupList_.push_back(id); 547 groupNames_.push_back(tag); 548 } 549 } 550 } 551 else if (strcmp(arg[argIdx],"type")==0) { 552 if (narg-argIdx == 2) { 553 int id = atoi(arg[++argIdx]); 554 typeList_.push_back(id); 555 typeNames_.push_back(tag); 556 } 557 else { 558 while (++argIdx < narg) { 559 int id = atoi(arg[argIdx]); 560 string tag = speciesTag+"_"+to_string(id); 561 typeList_.push_back(id); 562 typeNames_.push_back(tag); 563 } 564 } 565 } 566 else { 567 throw ATC_Error("ATC_Method: add_species only handles groups or types"); } 568 match = true; 569 } 570 571 // remove species from tracking 572 573 /*! \page man_remove_species fix_modify AtC remove_species 574 \section syntax 575 fix_modify_AtC delete_species <TAG> \n 576 577 - <TAG> = tag for tracking a species \n 578 \section examples 579 <TT> fix_modify AtC remove_species gold </TT> \n 580 \section description 581 Removes tag designated for tracking a specified species. \n 582 \section restrictions 583 \section related 584 \section default 585 No defaults for this command. 586 */ 587 else if (strcmp(arg[argIdx],"delete_species")==0) { 588 argIdx++; 589 string tag = arg[argIdx++]; 590 if (strcmp(arg[argIdx],"group")==0) { 591 for (unsigned int j = 0; j < groupList_.size(); j++) { 592 if (tag == groupNames_[j]) { 593 groupList_.erase(groupList_.begin()+j); 594 groupNames_.erase(groupNames_.begin()+j); 595 break; 596 } 597 } 598 } 599 else if (strcmp(arg[argIdx],"type")==0) { 600 for (unsigned int j = 0; j < typeList_.size(); j++) { 601 if (tag == typeNames_[j]) { 602 typeList_.erase(typeList_.begin()+j); 603 typeNames_.erase(typeNames_.begin()+j); 604 break; 605 } 606 } 607 } 608 else { 609 throw ATC_Error("ATC_Method: delete_species only handles groups or types"); } 610 match = true; 611 612 } 613 614 // add a molecule for tracking 615 /*! \page man_add_molecule fix_modify AtC add_molecule 616 \section syntax 617 fix_modify_AtC add_molecule <small|large> <TAG> <GROUP_NAME> \n 618 619 - small|large = can be small if molecule size < cutoff radius, must be large otherwise \n 620 - <TAG> = tag for tracking a species \n 621 - <GROUP_NAME> = name of group that tracking will be applied to \n 622 \section examples 623 <TT> group WATERGROUP type 1 2 </TT> \n 624 <TT> fix_modify AtC add_molecule small water WATERGROUP </TT> \n 625 \section description 626 Associates a tag with all molecules corresponding to a specified group. \n 627 \section restrictions 628 \section related 629 \section default 630 No defaults for this command. 631 */ 632 else if (strcmp(arg[argIdx],"add_molecule")==0) { 633 argIdx++; 634 MolSize size; 635 if (strcmp(arg[argIdx],"small")==0) { 636 size = MOL_SMALL; 637 //needXrefProcessorGhosts_ = true; 638 needProcGhost_ = true; 639 } 640 else 641 throw ATC_Error("ATC_CouplingMass: Bad molecule size in add_molecule"); 642 argIdx++; 643 string moleculeTag = arg[argIdx]; 644 645 argIdx++; 646 int groupBit = lammpsInterface_->group_bit(arg[argIdx]); 647 moleculeIds_[moleculeTag] = pair<MolSize,int>(size,groupBit); 648 match = true; 649 } 650 651 // remove molecule from tracking 652 /*! \page man_remove_molecule fix_modify AtC remove_molecule 653 \section syntax 654 fix_modify_AtC remove_molecule <TAG> \n 655 656 - <TAG> = tag for tracking a molecule type \n 657 \section examples 658 <TT> fix_modify AtC remove_molecule water </TT> \n 659 \section description 660 Removes tag designated for tracking a specified set of molecules. \n 661 \section restrictions 662 \section related 663 \section default 664 No defaults for this command. 665 */ 666 else if (strcmp(arg[argIdx],"remove_molecule")==0) { 667 argIdx++; 668 string moleculeTag = arg[argIdx]; 669 moleculeIds_.erase(moleculeTag); 670 671 taggedDensMan_.erase(moleculeTag); 672 } 673 674 /*! \page man_boundary fix_modify AtC boundary 675 \section syntax 676 fix_modify AtC boundary type <atom-type-id> 677 - <atom-type-id> = type id for atoms that represent a ficticious 678 boundary internal to the FE mesh 679 \section examples 680 <TT> fix_modify AtC boundary type ghost_atoms </TT> 681 \section description 682 Command to define the atoms that represent the ficticious 683 boundary internal to the FE mesh. For fully overlapped MD/FE 684 domains with periodic boundary conditions no boundary atoms should 685 be defined. 686 \section restrictions 687 \section default 688 none 689 */ 690 else if (strcmp(arg[argIdx],"boundary")==0) { 691 argIdx++; 692 groupTagGhost_ = arg[argIdx++]; 693 match = true; 694 } 695 696 /*! \page man_internal_atom_integrate fix_modify AtC internal_atom_integrate 697 \section syntax 698 fix_modify AtC internal_atom_integrate <on | off> 699 <TT> fix_modify AtC internal_atom_integrate on </TT> 700 \section description 701 Has AtC perform time integration for the atoms in the group on which it operates. This does not include boundary atoms. 702 \section restrictions 703 AtC must be created before any fixes doing time integration. It must be on for coupling methods which impose constraints on velocities during the first verlet step, e.g. control momentum glc_velocity. 704 \section default 705 on for coupling methods, off for post-processors 706 off 707 */ 708 else if (strcmp(arg[argIdx],"internal_atom_integrate")==0) { 709 argIdx++; 710 if (strcmp(arg[argIdx],"off")==0) { 711 integrateInternalAtoms_ = false; 712 match = true; 713 } 714 else { 715 integrateInternalAtoms_ = true; 716 match = true; 717 } 718 } 719 720 /*! \page man_internal_element_set fix_modify AtC internal_element_set 721 \section syntax 722 fix_modify AtC internal_element_set <element-set-name> 723 - <element-set-name> = name of element set defining internal region, or off 724 \section examples 725 <TT> fix_modify AtC internal_element_set myElementSet </TT> 726 <TT> fix_modify AtC internal_element_set off </TT> 727 \section description 728 Enables AtC to base the region for internal atoms to be an element set. 729 If no ghost atoms are used, all the AtC atoms must be constrained to remain 730 in this element set by the user, e.g., with walls. If boundary atoms are 731 used in conjunction with Eulerian atom maps 732 AtC will partition all atoms of a boundary or internal type to be of type internal 733 if they are in the internal region or to be of type boundary otherwise. 734 \section restrictions 735 If boundary atoms are used in conjunction with Eulerian atom maps, the Eulerian 736 reset frequency must be an integer multiple of the Lammps reneighbor frequency 737 \section related 738 see \ref atom_element_map_type and \ref boundary 739 \section default 740 off 741 */ 742 else if (strcmp(arg[argIdx],"internal_element_set")==0) { 743 argIdx++; 744 if (strcmp(arg[argIdx],"off")==0) { 745 internalElementSet_ = ""; 746 match = true; 747 } 748 else { 749 internalElementSet_ = string(arg[argIdx]); 750 const set<int> & elementSet((feEngine_->fe_mesh())->elementset(internalElementSet_)); // check it exists and is not trivial 751 if (elementSet.size()==0) throw ATC_Error("internal_element_set - element set " + internalElementSet_ + " has no elements"); 752 match = true; 753 } 754 } 755 756 /*! \page man_atom_weight fix_modify AtC atom_weight 757 \section syntax 758 fix_modify AtC atom_weight <method> <arguments> 759 - <method> = \n 760 value: atoms in specified group assigned constant value given \n 761 lattice: volume per atom for specified lattice type (e.g. fcc) and parameter \n 762 element: element volume divided among atoms within element \n 763 region: volume per atom determined based on the atom count in the MD regions and their volumes. Note: meaningful only if atoms completely fill all the regions. \n 764 group: volume per atom determined based on the atom count in a group and its volume\n 765 read_in: list of values for atoms are read-in from specified file \n 766 \section examples 767 <TT> fix_modify atc atom_weight constant myatoms 11.8 </TT> \n 768 <TT> fix_modify atc atom_weight lattice </TT> \n 769 <TT> fix_modify atc atom_weight read-in atm_wt_file.txt </TT> \n 770 \section description 771 Command for assigning the value of atomic weights used for atomic integration in 772 atom-continuum coupled simulations. 773 \section restrictions 774 Use of lattice option requires a lattice type and parameter is already specified. 775 \section related 776 \section default 777 lattice 778 */ 779 else if (strcmp(arg[argIdx],"atom_weight")==0) { 780 argIdx++; 781 if (strcmp(arg[argIdx],"constant")==0) { 782 argIdx++; 783 atomWeightType_ = USER; 784 int groupbit = -1; 785 if (strcmp(arg[argIdx],"all")==0) { 786 } 787 else { 788 groupbit = lammpsInterface_->group_bit(arg[argIdx]); 789 } 790 argIdx++; 791 double value = atof(arg[argIdx]); 792 Valpha_[groupbit] = value; 793 match = true; 794 } 795 else if (strcmp(arg[argIdx],"lattice")==0) { 796 atomWeightType_ = LATTICE; 797 match = true; 798 } 799 else if (strcmp(arg[argIdx],"element")==0) { 800 atomWeightType_ = ELEMENT; 801 match = true; 802 } 803 else if (strcmp(arg[argIdx],"region")==0) { 804 atomWeightType_ = REGION; 805 match = true; 806 } 807 else if (strcmp(arg[argIdx],"group")==0) { 808 atomWeightType_ = GROUP; 809 match = true; 810 } 811 else if (strcmp(arg[argIdx],"multiscale")==0) { 812 atomWeightType_ = MULTISCALE; 813 match = true; 814 } 815 else if (strcmp(arg[argIdx],"node")==0) { 816 atomWeightType_ = NODE; 817 match = true; 818 } 819 else if (strcmp(arg[argIdx],"node_element")==0) { 820 atomWeightType_ = NODE_ELEMENT; 821 match = true; 822 } 823 else if (strcmp(arg[argIdx],"read_in")==0) { 824 atomWeightType_ = READ_IN; 825 argIdx++; 826 atomicWeightsFile_ = arg[argIdx]; 827 match = true; 828 } 829 if (match) { 830 needReset_ = true; 831 } 832 } 833 834 /*! \page man_decomposition fix_modify AtC decomposition 835 \section syntax 836 fix_modify AtC decomposition <type> 837 - <type> = \n 838 replicated_memory: nodal information replicated on each processor \n 839 distributed_memory: only owned nodal information on processor \n 840 \section examples 841 <TT> fix_modify atc decomposition distributed_memory </TT> \n 842 \section description 843 Command for assigning the distribution of work and memory for parallel runs. 844 \section restrictions 845 replicated_memory is appropriate for simulations were the number of nodes << number of atoms 846 \section related 847 \section default 848 replicated_memory 849 */ 850 else if (strcmp(arg[argIdx],"decomposition")==0) { 851 argIdx++; 852 if (strcmp(arg[argIdx],"replicated_memory")==0) { 853 domainDecomposition_ = REPLICATED_MEMORY; 854 match = true; 855 } 856 else if (strcmp(arg[argIdx],"distributed_memory")==0) { 857 domainDecomposition_ = DISTRIBUTED_MEMORY; 858 match = true; 859 } 860 } 861 862 /*! \page man_write_atom_weights fix_modify AtC write_atom_weights 863 \section syntax 864 fix_modify AtC write_atom_weights <filename> <frequency> 865 - <filename> = name of file that atomic weights are written to \n 866 - <frequency> = how often writes will occur \n 867 \section examples 868 <TT> fix_modify atc write_atom_weights atm_wt_file.txt 10 </TT> \n 869 \section description 870 Command for writing the values of atomic weights to a specified file. 871 \section restrictions 872 \section related 873 \section default 874 */ 875 else if (strcmp(arg[argIdx],"write_atom_weights")==0) { 876 argIdx++; 877 atomicWeightsFile_ = arg[argIdx]; 878 argIdx++; 879 atomicWeightsWriteFrequency_ = atoi(arg[argIdx]); 880 atomicWeightsWriteFlag_ = true; 881 match = true; 882 } 883 884 885 /*! \page man_reset_time fix_modify AtC reset_time 886 \section syntax 887 fix_modify AtC reset_time <value> 888 \section examples 889 <TT> fix_modify atc reset_time 0.0 </TT> \n 890 \section description 891 Resets the simulation time counter. 892 \section restrictions 893 \section related 894 \section default 895 */ 896 else if (strcmp(arg[argIdx],"reset_time")==0) { 897 argIdx++; 898 set_time(); 899 if (narg > argIdx) { 900 double time = atof(arg[argIdx]); 901 set_time(time); 902 } 903 match = true; 904 } 905 906 /*! \page man_reset_time fix_modify AtC kernel_bandwidth 907 \section syntax 908 fix_modify AtC kernel_bandwidth <value> 909 \section examples 910 <TT> fix_modify atc reset_time 8 </TT> \n 911 \section description 912 Sets a maximum parallel bandwidth for the kernel functions during parallel communication. If the command is not issued, the default will be to assume the bandwidth of the kernel matrix corresponds to the number of sampling locations. 913 \section restrictions 914 Only is used if kernel functions are being used. 915 \section related 916 \section default 917 Number of sample locations. 918 */ 919 else if (strcmp(arg[argIdx],"kernel_bandwidth")==0) { 920 argIdx++; 921 accumulantBandwidth_ = atoi(arg[argIdx]); 922 match = true; 923 } 924 925 /*! \page man_reset_atomic_reference_positions fix_modify AtC reset_atomic_reference_positions 926 \section syntax 927 fix_modify AtC reset_atomic_reference_positions 928 \section examples 929 <TT> fix_modify atc reset_atomic_reference_positions 930 \section description 931 Resets the atomic positions ATC uses to perform point to field operations. 932 In can be used to use perfect lattice sites in ATC but a thermalized or 933 deformed lattice in LAMMPS. 934 \section restrictions 935 \section related 936 \section default 937 Default is off 938 */ 939 else if (strcmp(arg[argIdx],"reset_atomic_reference_positions")==0) { 940 argIdx++; 941 xRefFile_ = arg[argIdx]; 942 readXref_ = true; 943 match = true; 944 } 945 946 /*! \page man_set fix_modify AtC set 947 \section syntax 948 fix_modify AtC set reference_potential_energy <value_or_filename(optional)> 949 - value (double) : optional user specified zero point for PE in native LAMMPS energy units \n 950 - filename (string) : optional user specified string for file of nodal PE values to be read-in 951 \section examples 952 <TT> fix_modify AtC set reference_potential_energy </TT> \n 953 <TT> fix_modify AtC set reference_potential_energy -0.05 </TT> \n 954 <TT> fix_modify AtC set reference_potential_energy myPEvalues </TT> \n 955 \section description 956 Used to set various quantities for the post-processing algorithms. 957 It sets the zero point for the potential energy density using 958 the value provided for all nodes, or from the current 959 configuration of the lattice if no value is provided, or 960 values provided within the specified filename. 961 \section restrictions 962 Must be used with the hardy/field type of AtC fix 963 ( see \ref man_fix_atc ) 964 \section related 965 \section default 966 Defaults to lammps zero point i.e. isolated atoms 967 */ 968 else if (strcmp(arg[argIdx],"set")==0) { 969 argIdx++; 970 if (strcmp(arg[argIdx],"reference_potential_energy")==0) { 971 argIdx++; 972 setRefPE_ = true; 973 if (narg > argIdx) { 974 string a(arg[argIdx]); 975 if (is_dbl(a)) { 976 double value = atof(arg[argIdx]); 977 refPEvalue_ = value; 978 setRefPEvalue_ = true; 979 } 980 else { 981 nodalRefPEfile_ = arg[argIdx]; 982 readRefPE_ = true; 983 } 984 } 985 match = true; 986 } 987 } // end "set" 988 989 990 991 /*! \page man_atom_element_map fix_modify AtC atom_element_map 992 \section syntax 993 fix_modify AtC atom_element_map <eulerian|lagrangian> <frequency> \n 994 - frequency (int) : frequency of updating atom-to-continuum maps based on the 995 current configuration - only for eulerian 996 \section examples 997 <TT> fix_modify atc atom_element_map eulerian 100 </TT> 998 \section description 999 Changes frame of reference from eulerian to lagrangian and sets the 1000 frequency for which the map from atoms to elements is reformed and 1001 all the attendant data is recalculated. 1002 \section restrictions 1003 Cannot change map type after initialization. 1004 \section related 1005 \section default 1006 lagrangian 1007 */ 1008 else if (strcmp(arg[argIdx],"atom_element_map")==0) { 1009 argIdx++; 1010 if (strcmp(arg[argIdx],"eulerian")==0) { 1011 atomToElementMapType_ = EULERIAN; 1012 argIdx++; 1013 atomToElementMapFrequency_ = atoi(arg[argIdx]); 1014 } 1015 else { 1016 atomToElementMapType_ = LAGRANGIAN; 1017 atomToElementMapFrequency_ = 0; 1018 } 1019 match = true; 1020 needReset_ = true; 1021 } 1022 1023 /*! \page man_read_restart fix_modify AtC read_restart 1024 \section syntax 1025 fix_modify AtC read_restart [file_name] \n 1026 \section examples 1027 <TT> fix_modify AtC read_restart ATC_state </TT> \n 1028 \section description 1029 Reads the current state of the fields from a named text-based restart 1030 file. 1031 \section restrictions 1032 The restart file only contains fields and their time derivatives. 1033 The reference positions of the atoms and the commands that initialize 1034 the fix are not saved e.g. an identical mesh containing the same atoms 1035 will have to be recreated. 1036 \section related 1037 see write_restart \ref man_write_restart 1038 \section default 1039 none 1040 */ 1041 else if (strcmp(arg[argIdx],"read_restart")==0) { 1042 argIdx++; 1043 restartFileName_ = arg[argIdx]; 1044 useRestart_ = true; 1045 match = true; 1046 } 1047 1048 /*! \page man_write_restart fix_modify AtC write_restart 1049 \section syntax 1050 fix_modify AtC write_restart [file_name] \n 1051 \section examples 1052 <TT> fix_modify AtC write_restart restart.mydata </TT> \n 1053 \section description 1054 Dumps the current state of the fields to a named text-based restart file. 1055 This done when the command is invoked and not repeated, unlike the 1056 similar lammps command. 1057 \section restrictions 1058 The restart file only contains fields and their time derivatives. 1059 The reference positions of the atoms and the commands that initialize 1060 the fix are not saved e.g. an identical mesh containing the same atoms 1061 will have to be recreated. 1062 \section related 1063 see read_restart \ref man_read_restart 1064 \section default 1065 none 1066 */ 1067 else if (strcmp(arg[argIdx],"write_restart")==0) { 1068 argIdx++; 1069 string restartFileName(arg[argIdx]); 1070 RESTART_LIST data; 1071 write_restart_data(restartFileName,data); 1072 match = true; 1073 } 1074 1075 } // end else 1076 1077 return match; // return to FixATC 1078 1079 } 1080 1081 //-------------------------------------------------- 1082 // helper function for parser 1083 // handles : "displacement x" and "temperature" by indexing argIdx 1084 // for fluxes : only normal fluxes can be prescribed 1085 //-------------------------------------------------- parse_field(char ** args,int & argIdx,FieldName & thisField,int & thisIndex)1086 void ATC_Method::parse_field(/*const*/ char ** args, int & argIdx, 1087 FieldName & thisField, int & thisIndex) 1088 { 1089 string thisName = args[argIdx++]; 1090 if (args[argIdx] == NULL) { 1091 throw ATC_Error("Need to give field '"+thisName+"' more args"); 1092 } 1093 thisField = string_to_field(thisName); 1094 map<FieldName,int>::const_iterator iter = fieldSizes_.find(thisField); 1095 if (iter == fieldSizes_.end()) { 1096 throw ATC_Error("Bad field name: "+thisName); 1097 } 1098 string thisDim = args[argIdx]; 1099 thisIndex = 0; 1100 if (string_to_index(thisDim,thisIndex)) { 1101 if ( !( thisIndex < fieldSizes_[thisField]) ) { 1102 throw ATC_Error("Bad field dimension "+thisDim); 1103 } 1104 argIdx++; 1105 } 1106 } 1107 1108 //------------------------------------------------------------------- 1109 // this sets the peratom output 1110 update_peratom_output()1111 void ATC_Method::update_peratom_output() 1112 { 1113 perAtomArray_ = perAtomOutput_; 1114 // copy values 1115 for (int i = 0; i < lammpsInterface_->nlocal(); i++) { 1116 for (int j = 0; j < nsd_; j++) { 1117 perAtomOutput_[i][j] = xref_[i][j]; 1118 } 1119 for (int j = nsd_; j < sizePerAtomCols_; j++) { 1120 perAtomOutput_[i][j] = 0; 1121 } 1122 } 1123 int indx = nsd_; 1124 if (atomVolume_->nRows() > 0) { // kernel Hardy does not compute these 1125 const DIAG_MAT & myAtomicWeights(atomVolume_->quantity()); 1126 for (int i = 0; i < nLocal_; i++) { 1127 double wg = myAtomicWeights(i,i); 1128 if (wg > 0) { 1129 int ii = internalToAtom_(i); 1130 perAtomOutput_[ii][indx] = 1./wg; 1131 } 1132 } 1133 } 1134 } 1135 adjust_xref_pbc()1136 void ATC_Method::adjust_xref_pbc() 1137 { 1138 1139 int nlocal = lammpsInterface_->nlocal(); 1140 int xperiodic = lammpsInterface_->xperiodic(); 1141 int yperiodic = lammpsInterface_->yperiodic(); 1142 int zperiodic = lammpsInterface_->zperiodic(); 1143 double **x = lammpsInterface_->xatom(); 1144 double boxxlo,boxxhi; 1145 double boxylo,boxyhi; 1146 double boxzlo,boxzhi; 1147 1148 lammpsInterface_->box_bounds(boxxlo,boxxhi, 1149 boxylo,boxyhi, 1150 boxzlo,boxzhi); 1151 // bool changed = false; 1152 for (int i = 0; i < nlocal; i++) { 1153 if (xperiodic) { 1154 if (x[i][0] < boxxlo) { 1155 xref_[i][0] += Xprd_; 1156 // changed = true; 1157 } 1158 if (x[i][0] >= boxxhi) { 1159 xref_[i][0] -= Xprd_; 1160 // changed = true; 1161 } 1162 } 1163 1164 if (yperiodic) { 1165 if (x[i][1] < boxylo) { 1166 xref_[i][1] += Yprd_; 1167 // changed = true; 1168 } 1169 if (x[i][1] >= boxyhi) { 1170 xref_[i][1] -= Yprd_; 1171 // changed = true; 1172 } 1173 } 1174 1175 if (zperiodic) { 1176 if (x[i][2] < boxzlo) { 1177 xref_[i][2] += Zprd_; 1178 // changed = true; 1179 } 1180 if (x[i][2] >= boxzhi) { 1181 xref_[i][2] -= Zprd_; 1182 // changed = true; 1183 } 1184 } 1185 } 1186 1187 // propagate reset if needed 1188 if (atomToElementMapType_ == LAGRANGIAN) { 1189 if (atomCoarseGrainingPositions_) { 1190 atomCoarseGrainingPositions_->force_reset(); 1191 } 1192 } 1193 else if (atomReferencePositions_) { 1194 atomReferencePositions_->force_reset(); 1195 } 1196 1197 } 1198 //------------------------------------------------------------------- initialize()1199 void ATC_Method::initialize() 1200 { 1201 feEngine_->partition_mesh(); 1202 // initialized_ is set to true by derived class initialize() 1203 // localStep_ is a counter within a run or minimize 1204 localStep_ = 0; 1205 // STEP 1) get basic information data from Lammps/fix 1206 // 1a) group ids for all internal atoms 1207 groupbit_ = lammpsInterface_->group_bit(groupTag_); 1208 1209 // 1b) group ids for ghost atoms 1210 groupbitGhost_ = 0; 1211 if (!groupTagGhost_.empty()) { 1212 groupbitGhost_ = lammpsInterface_->group_bit(groupTagGhost_); 1213 } 1214 1215 // 1c) periodicity and box bounds/lengths 1216 if (!initialized_) { 1217 1218 lammpsInterface_->box_periodicity(periodicity[0], 1219 periodicity[1], 1220 periodicity[2]); 1221 lammpsInterface_->box_bounds(box_bounds[0][0],box_bounds[1][0], 1222 box_bounds[0][1],box_bounds[1][1], 1223 box_bounds[0][2],box_bounds[1][2]); 1224 for (int k = 0; k < nsd_; k++) { 1225 box_length[k] = box_bounds[1][k] - box_bounds[0][k]; 1226 } 1227 1228 lammpsInterface_->set_reference_box(); 1229 1230 // get periodicity data from lammps for parallel exchange to adjust for periodicity 1231 Xprd_ = lammpsInterface_->domain_xprd(); 1232 Yprd_ = lammpsInterface_->domain_yprd(); 1233 Zprd_ = lammpsInterface_->domain_zprd(); 1234 // box_length[0] = Xprd_; 1235 // box_length[1] = Yprd_; 1236 // box_length[2] = Zprd_; 1237 XY_ = lammpsInterface_->domain_xy(); 1238 XZ_ = lammpsInterface_->domain_xz(); 1239 YZ_ = lammpsInterface_->domain_yz(); 1240 } 1241 1242 // STEP 2 computational geometry 1243 // 2a) get basic information from continuum/FE 1244 this->set_continuum_data(); 1245 1246 // STEP 2b) set up data structures for computational geometry 1247 if (this->reset_methods()) { 1248 // clear memory manager 1249 interscaleManager_.clear_temporary_data(); 1250 atomVolume_ = NULL; 1251 1252 // reference positions and energy 1253 if (!initialized_) { 1254 double **x = lammpsInterface_->xatom(); 1255 for (int i = 0; i < lammpsInterface_->nmax(); i++) { 1256 for (int j = 0; j < nsd_; j++) { 1257 xref_[i][j] = x[i][j]; 1258 } 1259 } 1260 1261 // re-write non-ghosts' xref with values from a file 1262 if (readXref_) { 1263 bool success = read_atomic_ref_positions(xRefFile_.c_str()); 1264 if (!success) 1265 throw ATC_Error("Error reading atomic reference positions"); 1266 readXref_ = false; 1267 } 1268 1269 // set up maps from lammps to atc indexing 1270 reset_nlocal(); 1271 } 1272 1273 this->set_computational_geometry(); 1274 } 1275 1276 // 2c) basic data regarding atomic system, e.g. atom coordinates 1277 if (atomToElementMapType_ == EULERIAN) { 1278 reset_coordinates(); 1279 } 1280 1281 // STEP 3) set up variables which will be integrated in time 1282 this->construct_time_integration_data(); 1283 1284 // STEP 4) instantiate all the various specific algorithms and methods 1285 this->construct_methods(); 1286 1287 // STEP 5) construct dependency-managed data 1288 // 5b) all other transfer operators 1289 // needs to be done before every run in case options have changed or the atoms have been changed by the user 1290 1291 if (this->reset_methods()) { 1292 // construct all the needed data structures 1293 this->construct_transfers(); 1294 1295 // allocate all space needed for lammps arrays 1296 interscaleManager_.grow_arrays(lammpsInterface_->nmax()); 1297 } 1298 // reset all computes invoked flags and lammps data 1299 interscaleManager_.lammps_force_reset(); 1300 1301 // STEP 6) initialize data 1302 // 6b) size quantities which use pack_comm 1303 interscaleManager_.size_comm_quantities(); 1304 1305 // 6c) set coarse-graining functions and atomic weights 1306 if (!initialized_) { 1307 // FE_Engine allocates all required memory 1308 // assume initial atomic position is the reference position for now 1309 1310 // \int_\Omega N_I dV : static if the mesh is 1311 NodeVolumes_.reset(nNodes_,nNodes_); 1312 invNodeVolumes_.reset(nNodes_,nNodes_); 1313 feEngine_->compute_lumped_mass_matrix(NodeVolumes_); 1314 invNodeVolumes_.set_quantity() = NodeVolumes_.inv(); 1315 } 1316 atomVolume_->set_reset(); 1317 1318 // 6d) reference values 1319 this->set_reference_potential_energy(); 1320 1321 // 6e) atomic output for 0th step 1322 update_peratom_output(); 1323 1324 // clear need for resets 1325 needReset_ = false; 1326 1327 } 1328 //------------------------------------------------------------------- set_continuum_data()1329 void ATC_Method::set_continuum_data() 1330 { 1331 // initialize finite element engine and get basic properties 1332 if (!initialized_) { 1333 feEngine_->initialize(); 1334 if (nsd_!=feEngine_->nsd()) { 1335 throw ATC_Error("Spatial dimensions inconsistent between LAMMPS and ATC"); 1336 } 1337 nNodes_ = feEngine_->num_nodes(); 1338 } 1339 } 1340 1341 //------------------------------------------------------------------- set_computational_geometry()1342 void ATC_Method::set_computational_geometry() 1343 { 1344 // set positions used for coarse-graining operators 1345 1346 1347 1348 1349 if (!initialized_) { 1350 if (atomToElementMapType_ == EULERIAN) { 1351 FundamentalAtomQuantity * atomPositionsAll = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION,ALL); 1352 ClonedAtomQuantity<double> * myAtomPositions = 1353 new ClonedAtomQuantity<double>(this,atomPositionsAll,INTERNAL); 1354 atomCoarseGrainingPositions_ = myAtomPositions; 1355 interscaleManager_.add_per_atom_quantity(myAtomPositions, 1356 "AtomicCoarseGrainingPositions"); 1357 1358 if (trackDisplacement_) { 1359 XrefWrapper * myAtomReferencePositions = new XrefWrapper(this); 1360 atomReferencePositions_ = myAtomReferencePositions; 1361 interscaleManager_.add_per_atom_quantity(myAtomReferencePositions, 1362 "AtomicReferencePositions"); 1363 atomReferencePositions_->set_memory_type(PERSISTENT); 1364 } 1365 1366 if (groupbitGhost_) { 1367 myAtomPositions = new ClonedAtomQuantity<double>(this,atomPositionsAll,GHOST); 1368 atomGhostCoarseGrainingPositions_ = myAtomPositions; 1369 interscaleManager_.add_per_atom_quantity(myAtomPositions, 1370 "AtomicGhostCoarseGrainingPositions"); 1371 } 1372 if(needProcGhost_){ 1373 FundamentalAtomQuantity * atomPositionsAll = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION,PROC_GHOST); 1374 ClonedAtomQuantity<double> * myAtomPositions = 1375 new ClonedAtomQuantity<double>(this,atomPositionsAll,PROC_GHOST); 1376 atomProcGhostCoarseGrainingPositions_ = myAtomPositions; 1377 interscaleManager_.add_per_atom_quantity(myAtomPositions, 1378 "AtomicProcGhostCoarseGrainingPositions"); 1379 } 1380 } 1381 else { 1382 XrefWrapper * myAtomPositions = new XrefWrapper(this); 1383 atomCoarseGrainingPositions_ = myAtomPositions; 1384 interscaleManager_.add_per_atom_quantity(myAtomPositions, 1385 "AtomicCoarseGrainingPositions"); 1386 atomReferencePositions_ = atomCoarseGrainingPositions_; 1387 1388 if (groupbitGhost_) { 1389 myAtomPositions = new XrefWrapper(this,GHOST); 1390 atomGhostCoarseGrainingPositions_ = myAtomPositions; 1391 interscaleManager_.add_per_atom_quantity(myAtomPositions, 1392 "AtomicGhostCoarseGrainingPositions"); 1393 } 1394 if (needProcGhost_) { 1395 XrefWrapper * myAtomPositions = new XrefWrapper(this); 1396 atomProcGhostCoarseGrainingPositions_ = myAtomPositions; 1397 interscaleManager_.add_per_atom_quantity(myAtomPositions, 1398 "AtomicProcGhostCoarseGrainingPositions"); 1399 } 1400 } 1401 atomCoarseGrainingPositions_->set_memory_type(PERSISTENT); 1402 if (atomGhostCoarseGrainingPositions_) atomGhostCoarseGrainingPositions_->set_memory_type(PERSISTENT); 1403 if (atomProcGhostCoarseGrainingPositions_) atomProcGhostCoarseGrainingPositions_->set_memory_type(PERSISTENT); 1404 } 1405 1406 // Add in atom to element map if using shape functions 1407 if (needsAtomToElementMap_) { 1408 atomElement_ = new AtomToElementMap(this); 1409 interscaleManager_.add_per_atom_int_quantity(atomElement_,"AtomElement"); 1410 } 1411 } 1412 1413 //------------------------------------------------------------------- construct_methods()1414 void ATC_Method::construct_methods() 1415 { 1416 1417 if (this->reset_methods()) { 1418 if (atomTimeIntegrator_) delete atomTimeIntegrator_; 1419 if (integrateInternalAtoms_) { 1420 atomTimeIntegrator_ = new AtomTimeIntegratorType(this,INTERNAL); 1421 } 1422 else { 1423 atomTimeIntegrator_ = new AtomTimeIntegrator(); 1424 } 1425 1426 // set up integration schemes for ghosts 1427 ghostManager_.construct_methods(); 1428 } 1429 } 1430 1431 //------------------------------------------------------------------- construct_transfers()1432 void ATC_Method::construct_transfers() 1433 { 1434 this->construct_interpolant(); 1435 1436 this->construct_molecule_transfers(); 1437 1438 atomTimeIntegrator_->construct_transfers(); 1439 ghostManager_.construct_transfers(); 1440 1441 1442 1443 } 1444 //------------------------------------------------------------------- create_atom_volume()1445 PerAtomDiagonalMatrix<double> * ATC_Method::create_atom_volume() 1446 { 1447 if (atomVolume_) { 1448 return atomVolume_; 1449 } 1450 else { 1451 // set variables to compute atomic weights 1452 DENS_MAN * nodalVolume(NULL); 1453 switch (atomWeightType_) { 1454 case USER: 1455 atomVolume_ = new AtomVolumeUser(this,Valpha_); 1456 break; 1457 case LATTICE: 1458 atomVolume_ = new AtomVolumeLattice(this); 1459 break; 1460 case ELEMENT: 1461 atomVolume_ = new AtomVolumeElement(this); 1462 break; 1463 case REGION: 1464 atomVolume_ = new AtomVolumeRegion(this); 1465 break; 1466 case GROUP: 1467 atomVolume_ = new AtomVolumeGroup(this,Valpha_); 1468 break; 1469 case MULTISCALE: 1470 if (!shpFcn_) { 1471 throw ATC_Error("ATC_Method::create_atom_volume - Multiscale algorithm requires an interpolant"); 1472 } 1473 nodalVolume = new NodalAtomVolume(this,shpFcn_); 1474 interscaleManager_.add_dense_matrix(nodalVolume,"NodalAtomVolume"); 1475 atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_); 1476 break; 1477 case NODE: 1478 if (!shpFcn_) { 1479 throw ATC_Error("ATC_Method::create_atom_volume - Node algorithm requires an interpolant"); 1480 } 1481 nodalVolume = new NodalVolume(this,shpFcn_); 1482 interscaleManager_.add_dense_matrix(nodalVolume,"NodalVolume"); 1483 atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_); 1484 break; 1485 case NODE_ELEMENT: 1486 if (!shpFcn_) { 1487 throw ATC_Error("ATC_Method::create_atom_volume - Node-Element algorithm requires an interpolant"); 1488 } 1489 nodalVolume = new NodalAtomVolumeElement(this,shpFcn_); 1490 interscaleManager_.add_dense_matrix(nodalVolume,"NodalAtomVolumeElement"); 1491 atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_); 1492 break; 1493 case READ_IN: 1494 atomVolume_ = new AtomVolumeFile(this,atomicWeightsFile_); 1495 break; 1496 } 1497 if (atomVolume_) { 1498 interscaleManager_.add_per_atom_diagonal_matrix(atomVolume_,"AtomVolume"); 1499 } 1500 else { 1501 throw ATC_Error("ATC_Method::create_atom_volume - bad option for atom volume algorithm"); 1502 } 1503 1504 return atomVolume_; 1505 } 1506 } 1507 //-------------------------------------------------------- init_integrate()1508 void ATC_Method::init_integrate() 1509 { 1510 atomTimeIntegrator_->init_integrate_velocity(dt()); 1511 ghostManager_.init_integrate_velocity(dt()); 1512 // account for other fixes doing time integration 1513 interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_VELOCITY); 1514 1515 atomTimeIntegrator_->init_integrate_position(dt()); 1516 ghostManager_.init_integrate_position(dt()); 1517 // account for other fixes doing time integration 1518 interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_POSITION); 1519 } 1520 //------------------------------------------------------------------- post_init_integrate()1521 void ATC_Method::post_init_integrate() 1522 { 1523 ghostManager_.post_init_integrate(); 1524 } 1525 //------------------------------------------------------------------- pre_exchange()1526 void ATC_Method::pre_exchange() 1527 { 1528 adjust_xref_pbc(); 1529 // call interscale manager to sync atc per-atom data with lammps array ahead of parallel communication 1530 interscaleManager_.prepare_exchange(); 1531 1532 // change types based on moving from internal region to ghost region 1533 if ((atomToElementMapType_ == EULERIAN) && (step() % atomToElementMapFrequency_ == 0)) { 1534 ghostManager_.pre_exchange(); 1535 } 1536 } 1537 //------------------------------------------------------------------- setup_pre_exchange()1538 void ATC_Method::setup_pre_exchange() 1539 { 1540 adjust_xref_pbc(); 1541 // call interscale manager to sync atc per-atom data with lammps array ahead of parallel communication 1542 interscaleManager_.prepare_exchange(); 1543 } 1544 //------------------------------------------------------------------- pre_neighbor()1545 void ATC_Method::pre_neighbor() 1546 { 1547 // reset quantities arising from atom exchange 1548 reset_nlocal(); 1549 1550 interscaleManager_.post_exchange(); 1551 1552 // forward_comm should go here 1553 } 1554 //------------------------------------------------------------------- min_post_force()1555 void ATC_Method::min_post_force() 1556 { 1557 post_force(); 1558 } 1559 //------------------------------------------------------------------- post_force()1560 void ATC_Method::post_force() 1561 { 1562 // this resets allow for the possibility of other fixes modifying positions and velocities, e.g. walls, but reduces efficiency 1563 interscaleManager_.lammps_force_reset(); 1564 } 1565 //-------------------------------------------------------- final_integrate()1566 void ATC_Method::final_integrate() 1567 { 1568 atomTimeIntegrator_->final_integrate(dt()); 1569 ghostManager_.final_integrate(dt()); 1570 // account for other fixes doing time integration 1571 interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_VELOCITY); 1572 } 1573 //------------------------------------------------------------------- post_final_integrate()1574 void ATC_Method::post_final_integrate() 1575 { 1576 if (atomicWeightsWriteFlag_ && (step() % atomicWeightsWriteFrequency_ == 0)) { 1577 write_atomic_weights(atomicWeightsFile_,atomVolume_->quantity()); 1578 } 1579 } 1580 //------------------------------------------------------------------- end_of_step()1581 void ATC_Method::end_of_step() 1582 { 1583 localStep_ += 1; 1584 } 1585 //-------------------------------------------------------------- finish()1586 void ATC_Method::finish() 1587 { 1588 // FE Engine 1589 if (feEngine_) feEngine_->finish(); 1590 feEngine_->departition_mesh(); 1591 } 1592 1593 //-------------------------------------------------------------- 1594 /** method to add new fields to the included list */ 1595 //-------------------------------------------------------------- add_fields(map<FieldName,int> & newFieldSizes)1596 void ATC_Method::add_fields(map<FieldName,int> & newFieldSizes) 1597 { 1598 map<FieldName,int>::const_iterator field; 1599 for (field = newFieldSizes.begin(); field!=newFieldSizes.end(); field++) { 1600 FieldName thisField = field->first; 1601 int thisSize = field->second; 1602 if (fieldSizes_.find(thisField)==fieldSizes_.end()) { 1603 fieldSizes_[thisField] = thisSize; 1604 } 1605 } 1606 } 1607 1608 //------------------------------------------------------------------- set_reference_potential_energy(void)1609 void ATC_Method::set_reference_potential_energy(void) 1610 { 1611 if (setRefPE_) { 1612 if (setRefPEvalue_) { 1613 nodalRefPotentialEnergy_->set_quantity() = refPEvalue_; 1614 setRefPEvalue_ = false; 1615 } 1616 else if (readRefPE_) { 1617 if (LammpsInterface::instance()->rank_zero()) { 1618 stringstream ss; 1619 ss << "reading reference potential energy from " << nodalRefPEfile_; 1620 LammpsInterface::instance()->print_msg(ss.str()); 1621 } 1622 (nodalRefPotentialEnergy_->set_quantity()).from_file(nodalRefPEfile_); 1623 readRefPE_ = false; 1624 } 1625 else { 1626 hasRefPE_ = false; 1627 SPAR_MAN * referenceAccumulant = interscaleManager_.sparse_matrix("ReferenceAccumulant"); 1628 if (referenceAccumulant) { 1629 referenceAccumulant->set_quantity() = accumulant_->quantity(); 1630 } 1631 DIAG_MAN * referenceAccumulantInverseVolumes = interscaleManager_.diagonal_matrix("ReferenceAccumulantInverseVolumes"); 1632 if (referenceAccumulantInverseVolumes) { 1633 referenceAccumulantInverseVolumes->set_quantity() = accumulantInverseVolumes_->quantity(); 1634 } 1635 PAQ * atomicRefPe = interscaleManager_.per_atom_quantity("AtomicReferencePotential"); 1636 if (!atomicRefPe) { 1637 throw ATC_Error("ATC_Method::set_reference_potential_energy - atomic reference PE object was not created during construct_transfers"); 1638 } 1639 PAQ* pe = interscaleManager_.per_atom_quantity("AtomicPotentialEnergy"); 1640 if (!pe) { 1641 throw ATC_Error("ATC_Method::set_reference_potential_energy - atomic PE object was not created during construct_transfers"); 1642 } 1643 atomicRefPe->set_quantity() = pe->quantity(); 1644 atomicRefPe->fix_quantity(); 1645 } 1646 setRefPE_ = false; 1647 hasRefPE_ = true; 1648 } 1649 } 1650 //------------------------------------------------------------------- 1651 1652 1653 //================================================================= 1654 // memory management and processor information exchange 1655 //================================================================= 1656 1657 1658 //----------------------------------------------------------------- 1659 // number of doubles 1660 //----------------------------------------------------------------- doubles_per_atom() const1661 int ATC_Method::doubles_per_atom() const 1662 { 1663 1664 int doubles = 4; 1665 doubles += interscaleManager_.memory_usage(); 1666 return doubles; 1667 } 1668 1669 //----------------------------------------------------------------- 1670 // memory usage of local atom-based arrays 1671 //----------------------------------------------------------------- memory_usage()1672 int ATC_Method::memory_usage() 1673 { 1674 int bytes = doubles_per_atom(); 1675 bytes *= lammpsInterface_->nmax() * sizeof(double); 1676 return bytes; 1677 } 1678 1679 //----------------------------------------------------------------- 1680 // allocate local atom-based arrays 1681 //----------------------------------------------------------------- grow_arrays(int nmax)1682 void ATC_Method::grow_arrays(int nmax) 1683 { 1684 xref_ = 1685 lammpsInterface_->grow_2d_double_array(xref_,nmax,3,"fix_atc:xref"); 1686 1687 perAtomOutput_ = 1688 lammpsInterface_->grow_2d_double_array(perAtomOutput_,nmax,sizePerAtomCols_,"fix_atc:perAtomOutput"); 1689 interscaleManager_.grow_arrays(nmax); 1690 } 1691 1692 //----------------------------------------------------------------- 1693 // copy values within local atom-based arrays 1694 //----------------------------------------------------------------- copy_arrays(int i,int j)1695 void ATC_Method::copy_arrays(int i, int j) 1696 { 1697 xref_[j][0] = xref_[i][0]; 1698 xref_[j][1] = xref_[i][1]; 1699 xref_[j][2] = xref_[i][2]; 1700 1701 for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) { 1702 perAtomOutput_[j][ii] = perAtomOutput_[i][ii]; 1703 } 1704 interscaleManager_.copy_arrays(i,j); 1705 } 1706 1707 //----------------------------------------------------------------- 1708 // pack values in local atom-based arrays for exchange with another proc 1709 //----------------------------------------------------------------- pack_exchange(int i,double * buf)1710 int ATC_Method::pack_exchange(int i, double *buf) 1711 { 1712 buf[0] = xref_[i][0]; 1713 buf[1] = xref_[i][1]; 1714 buf[2] = xref_[i][2]; 1715 1716 int j = 4; 1717 for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) { 1718 buf[j++] = perAtomOutput_[i][ii]; 1719 } 1720 int interscaleSizeComm = interscaleManager_.pack_exchange(i,&buf[j]); 1721 return sizeComm_ + interscaleSizeComm; 1722 } 1723 1724 //----------------------------------------------------------------- 1725 // unpack values in local atom-based arrays from exchange with another proc 1726 //----------------------------------------------------------------- unpack_exchange(int nlocal,double * buf)1727 int ATC_Method::unpack_exchange(int nlocal, double *buf) 1728 { 1729 xref_[nlocal][0] = buf[0]; 1730 xref_[nlocal][1] = buf[1]; 1731 xref_[nlocal][2] = buf[2]; 1732 1733 int j = 4; 1734 for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) { 1735 perAtomOutput_[nlocal][ii] = buf[j++]; 1736 } 1737 int interscaleSizeComm = interscaleManager_.unpack_exchange(nlocal,&buf[j]); 1738 return sizeComm_ + interscaleSizeComm; 1739 } 1740 1741 //----------------------------------------------------------------- 1742 // pack values in local atom-based arrays from exchange with another proc 1743 //----------------------------------------------------------------- pack_comm(int n,int * list,double * buf,int pbc_flag,int * pbc)1744 int ATC_Method::pack_comm(int n, int *list, double *buf, 1745 int pbc_flag, int *pbc) 1746 { 1747 int i,j,m; 1748 double dx = 0,dy = 0,dz = 0; 1749 1750 int * num_bond = lammpsInterface_->num_bond(); 1751 int ** bond_atom = lammpsInterface_->bond_atom(); 1752 1753 m = 0; 1754 if (pbc_flag == 0) { 1755 for (i = 0; i < n; i++) { 1756 j = list[i]; 1757 buf[m++] = xref_[j][0]; 1758 buf[m++] = xref_[j][1]; 1759 buf[m++] = xref_[j][2]; 1760 1761 if (num_bond) { 1762 buf[m++] = num_bond[j]; 1763 for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) { 1764 buf[m++] = bond_atom[j][ii]; 1765 } 1766 } 1767 } 1768 } 1769 else { 1770 if (lammpsInterface_->domain_triclinic() == 0) { 1771 dx = pbc[0]*Xprd_; 1772 dy = pbc[1]*Yprd_; 1773 dz = pbc[2]*Zprd_; 1774 } 1775 else { 1776 dx = pbc[0]*Xprd_ + pbc[5]*XY_ + pbc[4]*XZ_; 1777 dy = pbc[1]*Yprd_ + pbc[3]*YZ_; 1778 dz = pbc[2]*Zprd_; 1779 } 1780 for (i = 0; i < n; i++) { 1781 j = list[i]; 1782 buf[m++] = xref_[j][0] + dx; 1783 buf[m++] = xref_[j][1] + dy; 1784 buf[m++] = xref_[j][2] + dz; 1785 1786 if (num_bond) { 1787 buf[m++] = num_bond[j]; 1788 for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) { 1789 buf[m++] = bond_atom[j][ii]; 1790 } 1791 } 1792 1793 } 1794 } 1795 1796 int mySize = 3; 1797 if (num_bond) 1798 mySize += 1 + lammpsInterface_->bond_per_atom(); 1799 return mySize; 1800 } 1801 1802 //----------------------------------------------------------------- 1803 // unpack values in local atom-based arrays from exchange with another proc 1804 //----------------------------------------------------------------- unpack_comm(int n,int first,double * buf)1805 void ATC_Method::unpack_comm(int n, int first, double *buf) 1806 { 1807 int i,m,last; 1808 1809 int * num_bond = lammpsInterface_->num_bond(); 1810 int ** bond_atom = lammpsInterface_->bond_atom(); 1811 1812 m = 0; 1813 last = first + n; 1814 for (i = first; i < last; i++) { 1815 xref_[i][0] = buf[m++]; 1816 xref_[i][1] = buf[m++]; 1817 xref_[i][2] = buf[m++]; 1818 1819 if (num_bond) { 1820 num_bond[i] = static_cast<int>(buf[m++]); 1821 for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) { 1822 bond_atom[i][ii] = static_cast<int>(buf[m++]); 1823 } 1824 } 1825 1826 } 1827 } 1828 1829 //----------------------------------------------------------------- 1830 // 1831 //----------------------------------------------------------------- reset_nlocal()1832 void ATC_Method::reset_nlocal() 1833 { 1834 nLocalTotal_ = lammpsInterface_->nlocal(); 1835 const int * mask = lammpsInterface_->atom_mask(); 1836 nLocal_ = 0; 1837 nLocalGhost_ = 0; 1838 1839 for (int i = 0; i < nLocalTotal_; ++i) { 1840 if (mask[i] & groupbit_) nLocal_++; 1841 if (mask[i] & groupbitGhost_) nLocalGhost_++; 1842 } 1843 1844 // set up internal & ghost maps 1845 1846 if (nLocal_>0) { 1847 // set map 1848 internalToAtom_.resize(nLocal_); 1849 int j = 0; 1850 // construct internalToAtom map 1851 // : internal index -> local lammps atom index 1852 for (int i = 0; i < nLocalTotal_; ++i) { 1853 if (mask[i] & groupbit_) internalToAtom_(j++) = i; 1854 } 1855 #ifdef EXTENDED_ERROR_CHECKING 1856 stringstream ss; 1857 ss << "Nlocal = " << nLocal_ << " but only found " << j << "atoms"; 1858 if (j!=nLocal_) throw ATC_Error(ss.str()); 1859 #endif 1860 // construct reverse map 1861 atomToInternal_.clear(); 1862 for (int i = 0; i < nLocal_; ++i) { 1863 atomToInternal_[internalToAtom_(i)] = i; 1864 } 1865 } 1866 if (nLocalGhost_>0) { 1867 // set map 1868 ghostToAtom_.resize(nLocalGhost_); 1869 int j = 0; 1870 for (int i = 0; i < nLocalTotal_; ++i) { 1871 if (mask[i] & groupbitGhost_) ghostToAtom_(j++) = i; 1872 } 1873 } 1874 1875 //WIP_JAT this should not be needed at all, but a memory problem with sparse matrices requires them to be reset (possibly related to note in SparseMatrix-inl.h::_delete()) 1876 interscaleManager_.reset_nlocal(); 1877 1878 } 1879 1880 //------------------------------------------------------------------- reset_coordinates()1881 void ATC_Method::reset_coordinates() 1882 { 1883 // update coarse graining positions for internal and ghost atoms 1884 atomCoarseGrainingPositions_->unfix_quantity(); 1885 atomCoarseGrainingPositions_->quantity(); 1886 atomCoarseGrainingPositions_->fix_quantity(); 1887 if (atomGhostCoarseGrainingPositions_) { 1888 atomGhostCoarseGrainingPositions_->unfix_quantity(); 1889 atomGhostCoarseGrainingPositions_->quantity(); 1890 atomGhostCoarseGrainingPositions_->fix_quantity(); 1891 } 1892 if (atomProcGhostCoarseGrainingPositions_) { 1893 atomProcGhostCoarseGrainingPositions_->unfix_quantity(); 1894 atomProcGhostCoarseGrainingPositions_->quantity(); 1895 atomProcGhostCoarseGrainingPositions_->fix_quantity(); 1896 } 1897 } 1898 1899 //----------------------------------------------------------------- 1900 // 1901 //----------------------------------------------------------------- write_atomic_weights(const string filename,const DIAG_MAT & atomicVolumeMatrix)1902 void ATC_Method::write_atomic_weights(const string filename, const DIAG_MAT & atomicVolumeMatrix) 1903 { 1904 int nlocal = lammpsInterface_->nlocal(); 1905 int nlocalmax; 1906 LammpsInterface::instance()->int_allmax(&nlocal,&nlocalmax); 1907 int natoms = int(lammpsInterface_->natoms()); 1908 ofstream out; 1909 const char* fname = &filename[0]; 1910 1911 // create tag to local id map for this processor 1912 map <int,int> id2tag; 1913 map <int,int>::const_iterator itr; 1914 int * atom_tag = lammpsInterface_->atom_tag(); 1915 for (int i = 0; i < nlocal; ++i) { 1916 id2tag[i] = atom_tag[i]; 1917 } 1918 1919 int comm_rank = LammpsInterface::instance()->comm_rank(); 1920 int nprocs; 1921 LammpsInterface::instance()->int_allmax(&comm_rank,&nprocs); 1922 nprocs += 1; 1923 1924 if (comm_rank == 0) { 1925 out.open(fname); 1926 // print header lines 1927 out << "Atomic Weights for LAMMPS/atc analysis\n"; 1928 out << " \n"; 1929 out << natoms << " Atoms in system\n"; 1930 out << " \n"; 1931 // print atomic weights from proc 0 1932 for(int i = 0; i < nlocal; i++) { 1933 out << id2tag[i] << " " << atomicVolumeMatrix(i,i) << "\n"; 1934 } 1935 } 1936 1937 if (nprocs > 1) { 1938 int max_size,send_size; 1939 send_size = nlocal; 1940 LammpsInterface::instance()->int_allmax(&send_size,&max_size); 1941 1942 if (comm_rank == 0) { 1943 int intbuf[max_size]; 1944 double buf[max_size]; 1945 for (int iproc = 1; iproc < nprocs; iproc++) { 1946 LammpsInterface::instance()->int_recv(intbuf,max_size,iproc); 1947 LammpsInterface::instance()->recv(buf,max_size,iproc); 1948 for (int i = 0; i < max_size; i++) { 1949 out << intbuf[i] << " " << buf[i] << "\n"; 1950 } 1951 } 1952 } else { 1953 int intbuf[send_size]; 1954 double buf[send_size]; 1955 for (int i = 0; i < send_size; i++) { 1956 intbuf[i] = id2tag[i]; 1957 buf[i] = atomicVolumeMatrix(i,i); 1958 } 1959 LammpsInterface::instance()->int_send(intbuf,send_size); 1960 LammpsInterface::instance()->send(buf,send_size); 1961 } 1962 } 1963 1964 if (comm_rank == 0) { 1965 out.close(); 1966 } 1967 } 1968 1969 //----------------------------------------------------------------- 1970 // 1971 //----------------------------------------------------------------- compute_consistent_md_mass_matrix(const SPAR_MAT & shapeFunctionMatrix,SPAR_MAT & mdMassMatrix) const1972 void ATC_Method::compute_consistent_md_mass_matrix(const SPAR_MAT & shapeFunctionMatrix, 1973 SPAR_MAT & mdMassMatrix) const 1974 { 1975 1976 int nCols = shapeFunctionMatrix.nCols(); 1977 DENS_MAT massMatrixLocal(nCols,nCols); 1978 DENS_MAT denseMdMassMatrix(nCols,nCols); 1979 if (nLocal_>0) 1980 massMatrixLocal = shapeFunctionMatrix.transMat(shapeFunctionMatrix); 1981 1982 lammpsInterface_->allsum(massMatrixLocal.ptr(), 1983 denseMdMassMatrix.ptr(), 1984 denseMdMassMatrix.size()); 1985 mdMassMatrix.reset(denseMdMassMatrix,1.e-10); 1986 } 1987 1988 //================================================================= 1989 // Interscale operators 1990 //================================================================= 1991 // in the spirit of the current design of ATC: atoms local, nodes global 1992 1993 1994 1995 nodal_influence(const int groupbit,set<int> & nset,set<int> & aset,double tol)1996 bool ATC_Method::nodal_influence(const int groupbit, 1997 set<int> & nset, set<int> & aset, double tol) 1998 { 1999 int nghost = nodal_influence(groupbit,nset,aset,true,tol); 2000 int local_nghost = nghost; 2001 lammpsInterface_->int_allsum(&local_nghost,&nghost); 2002 if (nghost == 0) { 2003 nodal_influence(groupbit,nset,aset,false,tol); 2004 } 2005 return (nghost > 0); 2006 } nodal_influence(const int groupbit,set<int> & nset,set<int> & aset,bool ghost,double tol)2007 int ATC_Method::nodal_influence(const int groupbit, 2008 set<int> & nset, set<int> & aset, bool ghost, double tol) 2009 { 2010 Array<int> & amap = (ghost) ? ghostToAtom_ : internalToAtom_; 2011 int natoms = (ghost) ? nLocalGhost_ : nLocal_; 2012 DENS_MAT influence(nNodes_,1); 2013 DENS_MAT atomInfluence(natoms,1); 2014 const int *mask = lammpsInterface_->atom_mask(); 2015 for (int i = 0; i < natoms; i++) { 2016 if (mask[amap(i)] & groupbit) { 2017 atomInfluence(i,0) = 1; 2018 aset.insert(i); 2019 } 2020 } 2021 // relies on shape functions 2022 if (ghost) { 2023 restrict_volumetric_quantity(atomInfluence,influence,(interscaleManager_.per_atom_sparse_matrix("InterpolantGhost"))->quantity()); 2024 } 2025 else { 2026 restrict_volumetric_quantity(atomInfluence,influence); 2027 } 2028 2029 DENS_MAT localInfluence = influence; 2030 lammpsInterface_->allsum(localInfluence.ptr(), 2031 influence.ptr(), 2032 influence.size()); 2033 2034 for (int i = 0; i < nNodes_; i++) { 2035 if (fabs(influence(i,0)) > tol) { nset.insert(i); } 2036 } 2037 return aset.size(); 2038 } 2039 2040 2041 //-------------------------------------------------------- restrict_volumetric_quantity(const MATRIX & atomData,MATRIX & nodeData,const SPAR_MAT & shpFcn)2042 void ATC_Method::restrict_volumetric_quantity(const MATRIX & atomData, 2043 MATRIX & nodeData, 2044 const SPAR_MAT & shpFcn) 2045 { 2046 // computes nodeData = N*DeltaVAtom*atomData where N are the shape functions 2047 DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols()); 2048 //DENS_MAT workNodeArray; 2049 2050 2051 if (atomData.nRows() > 0) { // or shpFcn_??? 2052 workNodeArray = shpFcn.transMat(atomData); 2053 } 2054 int count = nodeData.nRows()*nodeData.nCols(); 2055 lammpsInterface_->allsum(workNodeArray.ptr(),nodeData.ptr(),count); 2056 return; 2057 } 2058 2059 2060 //-------------------------------------------------------- restrict_volumetric_quantity(const MATRIX & atomData,MATRIX & nodeData)2061 void ATC_Method::restrict_volumetric_quantity(const MATRIX & atomData, 2062 MATRIX & nodeData) 2063 { 2064 restrict_volumetric_quantity(atomData,nodeData,shpFcn_->quantity()); 2065 return; 2066 } 2067 2068 2069 //-------------------------------------------------------- prolong(const MATRIX & nodeData,MATRIX & atomData)2070 void ATC_Method::prolong(const MATRIX & nodeData, 2071 MATRIX & atomData) 2072 { 2073 // computes the finite element interpolation at atoms atomData = N*nodeData 2074 if (nLocal_>0) { 2075 atomData = (shpFcn_->quantity())*nodeData; 2076 } 2077 return; 2078 } 2079 2080 //======================================================== 2081 // FE functions 2082 //======================================================== 2083 2084 //-------------------------------------------------------- output()2085 void ATC_Method::output() 2086 { 2087 // if (lammpsInterface_->comm_rank() == 0) { 2088 compute_nodeset_output(); 2089 compute_faceset_output(); 2090 compute_elementset_output(); 2091 // } 2092 } 2093 //-------------------------------------------------------- compute_nodeset_output(void)2094 void ATC_Method::compute_nodeset_output(void) 2095 { 2096 map< pair <string, FieldName>, NodesetOperationType >::const_iterator iter; 2097 for (iter = nsetData_.begin(); iter != nsetData_.end();iter++){ 2098 pair <string, FieldName> id = iter->first; 2099 string nsetName = id.first; 2100 FieldName field = id.second; 2101 double sum = 0.0; 2102 const set<int> nset = feEngine_->fe_mesh()->nodeset(nsetName); 2103 const DENS_MAT & thisField = (fields_.find(field)->second).quantity(); 2104 set< int >::const_iterator itr; 2105 for (itr = nset.begin(); itr != nset.end();itr++){ 2106 int node = *itr; 2107 sum += thisField(node,0); 2108 } 2109 string name = nsetName + "_" + field_to_string(field); 2110 if (iter->second == NODESET_AVERAGE) { 2111 sum /= nset.size(); 2112 name = "average_"+name; 2113 } 2114 feEngine_->add_global(name, sum); 2115 } 2116 } 2117 //-------------------------------------------------------- compute_faceset_output(void)2118 void ATC_Method::compute_faceset_output(void) 2119 { 2120 map < pair<string,string>, FacesetIntegralType >::const_iterator iter; 2121 DENS_MAT values; 2122 for (iter = fsetData_.begin(); iter != fsetData_.end(); iter++) { 2123 string bndyName = (iter->first).first; 2124 string fieldName = (iter->first).second; 2125 const set< PAIR > & faceSet = (feEngine_->fe_mesh())->faceset(bndyName); 2126 ATOMIC_DATA::iterator data_iter = filteredData_.find(fieldName); 2127 if (data_iter == filteredData_.end()) { 2128 string msg = "Specified fieldName "+fieldName+ 2129 " not found in filteredData_ while attempting surface integration"; 2130 throw ATC_Error(msg); 2131 } 2132 const DENS_MAT & data = ((data_iter->second).quantity()); 2133 string stem = bndyName + "_" + fieldName + "_"; 2134 bool tf = false; 2135 if (iter->second == CONTOUR_INTEGRAL) { 2136 stem = "contour_"+stem; 2137 tf = true; 2138 } 2139 feEngine_->field_surface_flux(data,faceSet,values,tf); 2140 for (int i = 0; i < values.nRows() ; ++i ) { 2141 string name = stem + to_string(i+1); 2142 feEngine_->add_global(name, values(i,0)); 2143 } 2144 } 2145 } 2146 //-------------------------------------------------------- compute_elementset_output(void)2147 void ATC_Method::compute_elementset_output(void) 2148 { 2149 map< pair <string, FieldName>, ElementsetOperationType >::const_iterator iter; 2150 for (iter = esetData_.begin(); iter != esetData_.end();iter++){ 2151 pair <string, FieldName> id = iter->first; 2152 string esetName = id.first; 2153 FieldName field = id.second; 2154 const ESET eset = feEngine_->fe_mesh()->elementset(esetName); 2155 const DENS_MAT & thisField = (fields_.find(field)->second).quantity(); 2156 DENS_VEC total = feEngine_->integrate(thisField,eset); 2157 string name = esetName + "_" + field_to_string(field); 2158 if (iter->second == ELEMENTSET_AVERAGE) { 2159 DENS_MAT ones(nNodes_,0); ones = 1; 2160 DENS_VEC V = feEngine_->integrate(ones,eset); 2161 total /= V[0]; 2162 name = "average_"+name; 2163 } 2164 if (total.size() == 1) { 2165 feEngine_->add_global(name, total[0]); 2166 } 2167 else { 2168 for (int i = 0; i < total.size(); i++) { 2169 feEngine_->add_global(name+to_string(i), total[i]); 2170 } 2171 } 2172 } 2173 } 2174 2175 2176 //================================================================= 2177 // 2178 //================================================================= 2179 //-------------------------------------------------------- read_atomic_ref_positions(const char * filename)2180 bool ATC_Method::read_atomic_ref_positions(const char * filename) 2181 { 2182 int nlocal = lammpsInterface_->nlocal(); 2183 ifstream in; 2184 const int lineSize = 256; 2185 char line[lineSize]; 2186 2187 // create tag to local id map for this processor 2188 map <int,int> tag2id; 2189 map <int,int>::const_iterator itr; 2190 int * atom_tag = lammpsInterface_->atom_tag(); 2191 for (int i = 0; i < nlocal; ++i) { 2192 tag2id[atom_tag[i]] = i; 2193 } 2194 2195 // get number of atoms 2196 int natoms = 0; 2197 if (LammpsInterface::instance()->rank_zero()) { 2198 in.open(filename); 2199 string msg; 2200 string name = filename; 2201 msg = "no "+name+" reference file found"; 2202 if (! in.good()) throw ATC_Error(msg); 2203 in.getline(line,lineSize); // header 2204 in.getline(line,lineSize); // blank line 2205 in >> natoms; 2206 in.close(); 2207 stringstream ss; 2208 ss << "found " << natoms << " atoms in reference file"; 2209 LammpsInterface::instance()->print_msg(ss.str()); 2210 } 2211 LammpsInterface::instance()->int_broadcast(&natoms); 2212 2213 // read atoms and assign 2214 if (LammpsInterface::instance()->rank_zero()) { 2215 in.open(filename); 2216 while(in.good()) { 2217 in.getline(line,lineSize); 2218 string str(line); 2219 int pos = str.find("Atoms"); 2220 if (pos > -1) { 2221 in.getline(line,lineSize); // blank line 2222 break; 2223 } 2224 } 2225 } 2226 int nread = 0,type = -1, tag = -1, count = 0; 2227 double x[3]={0,0,0}; 2228 while (nread < natoms) { 2229 if (LammpsInterface::instance()->rank_zero()) { 2230 in.getline(line,lineSize); 2231 stringstream ss (line,stringstream::in | stringstream::out); 2232 ss >> tag >> type >> x[0] >> x[1] >> x[2]; 2233 nread++; 2234 } 2235 LammpsInterface::instance()->int_broadcast(&nread); 2236 LammpsInterface::instance()->int_broadcast(&tag); 2237 LammpsInterface::instance()->broadcast(x,3); 2238 itr = tag2id.find(tag); 2239 if (itr != tag2id.end()) { 2240 int iatom = itr->second; 2241 xref_[iatom][0] = x[0]; 2242 xref_[iatom][1] = x[1]; 2243 xref_[iatom][2] = x[2]; 2244 count++; 2245 } 2246 } 2247 if (LammpsInterface::instance()->rank_zero()) { 2248 in.close(); 2249 stringstream ss; 2250 ss << "read " << natoms << " reference positions"; 2251 LammpsInterface::instance()->print_msg(ss.str()); 2252 } 2253 if (count != nlocal) 2254 throw ATC_Error("reset "+to_string(count)+" atoms vs "+to_string(nlocal)); 2255 2256 2257 return true; 2258 } 2259 2260 //-------------------------------------------------------- remap_ghost_ref_positions(void)2261 void ATC_Method::remap_ghost_ref_positions(void) 2262 { 2263 2264 int nlocal = lammpsInterface_->nlocal(); 2265 int nghost = lammpsInterface_->nghost(); 2266 2267 double box_bounds[2][3]; 2268 lammpsInterface_->box_bounds(box_bounds[0][0],box_bounds[1][0], 2269 box_bounds[0][1],box_bounds[1][1], 2270 box_bounds[0][2],box_bounds[1][2]); 2271 double xlo = box_bounds[0][0], xhi = box_bounds[1][0]; 2272 double ylo = box_bounds[0][1], yhi = box_bounds[1][1]; 2273 double zlo = box_bounds[0][2], zhi = box_bounds[1][2]; 2274 2275 double box_length[3]; 2276 for (int k = 0; k < 3; k++) { 2277 box_length[k] = box_bounds[1][k] - box_bounds[0][k]; 2278 } 2279 double Lx = box_length[0], Ly = box_length[1], Lz = box_length[2]; 2280 2281 // create tag to local id map for this processor 2282 map <int,int> tag2id; 2283 map <int,int>::const_iterator itr; 2284 int * atom_tag = lammpsInterface_->atom_tag(); 2285 for (int i = 0; i < nlocal; ++i) { 2286 tag2id[atom_tag[i]] = i; 2287 } 2288 2289 // loop over ghosts 2290 double ** x = lammpsInterface_->xatom(); 2291 for (int j = nlocal; j < nlocal + nghost; j++) { 2292 int tag = atom_tag[j]; 2293 int i = tag2id[tag]; 2294 //itr = tag2id.find(tag); 2295 //if (itr != tag2id.end()) 2296 double* xj = x[j]; 2297 double* Xj = xref_[j]; 2298 //double Xj[3]; 2299 double* Xi = xref_[i]; 2300 // the assumption is that xref_[j] has been shuffled 2301 // so make an image of xref_[i] that is close to x[j] 2302 if (xj[0] <= xlo) Xj[0] = Xi[0] -Lx; 2303 if (xj[0] >= xhi) Xj[0] = Xi[0] +Lx; 2304 if (xj[1] <= ylo) Xj[1] = Xi[1] -Ly; 2305 if (xj[1] >= yhi) Xj[1] = Xi[1] +Ly; 2306 if (xj[2] <= zlo) Xj[2] = Xi[2] -Lz; 2307 if (xj[2] >= zhi) Xj[2] = Xi[2] +Lz; 2308 } 2309 } 2310 }; 2311