// ATC_Transfer headers #include "ATC_Transfer.h" #include "FE_Engine.h" #include "Array.h" #include "Array2D.h" #include "ATC_Error.h" #include "CG.h" #include "XT_Function.h" #include "PrescribedDataManager.h" #include "TimeIntegrator.h" #include "PhysicsModel.h" #include "PhysicsModelThermal.h" #include "PhysicsModelTwoTemperature.h" namespace ATC { ATC_Transfer::ATC_Transfer(void) : lammpsInterface_(LammpsInterface::instance()), physicsModel_(NULL), feEngine_(NULL), prescribedDataMgr_(NULL), simTime_(0.0), stepCounter_(0), sampleCounter_(0), outputFrequency_(0), sampleFrequency_(1), outputFrequencyAtom_(0), nLocal_(0), nLocalTotal_(0), nLocalGhost_(0), nLocalMask_(0), nInternal_(0), nGhost_(0), nLocalLambda_(0), equilibriumStart_(false), mdOutputManager_(), bndyIntType_(NO_QUADRATURE), bndyFaceSet_(NULL), globalSwitch_(0), resetSwitch_(false), atomToElementMapType_(LAGRANGIAN), atomToElementMapFrequency_(0), neighborResetFrequency_(0), scalarFlag_(0), vectorFlag_(0), sizeVector_(0), globalFreq_(0), extScalar_(0), extVector_(0), extList_(NULL), timeFilterManager_(this), groupbit_(0), groupbitGhost_(0), atomSwitch_(false), initialized_(false), atomQuadForInternal_(true), useLocalizedLambda_(false), useLumpedLambda_(false), regionID_(-1), atomWeightType_(LATTICE), xref_(NULL), readXref_(false), timeIntegrator_(NULL), extrinsicModelManager_(this), useRestart_(false), trackCharge_(false) { // Defaults grow_arrays(lammpsInterface_->nmax()); fieldMask_.reset(NUM_FIELDS,NUM_FLUX); fieldMask_ = false; // fe_engine feEngine_ = new FE_Engine(this); // check to see if lammps has any charges double * atomicCharge = lammpsInterface_->atom_charge(); if (atomicCharge) trackCharge_ = true; } ATC_Transfer::~ATC_Transfer() { lammpsInterface_->destroy_2d_double_array(xref_); if (feEngine_) delete feEngine_; if (physicsModel_) delete physicsModel_; if (prescribedDataMgr_) delete prescribedDataMgr_; } //-------------------------------------------------- // pack_fields // bundle all allocated field matrices into a list // for output needs //-------------------------------------------------- void ATC_Transfer::pack_fields(OUTPUT_LIST & data) { map::const_iterator field; for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) { FieldName thisField = field->first; string fieldName = field_to_string(thisField); string matrixName; // copy all fields from ATC_Transfer.h matrixName = "fields_" + fieldName; data[matrixName] = & fields_[thisField]; matrixName = "dot_fields_" + fieldName; data[matrixName] = & dot_fields_[thisField]; matrixName = "ddot_fields_" + fieldName; data[matrixName] = &ddot_fields_[thisField]; matrixName = "dddot_fields_" + fieldName; data[matrixName] = & dddot_fields_[thisField]; matrixName = "dot_fieldsMD_" + fieldName; data[matrixName] = & dot_fieldsMD_[thisField]; matrixName = "ddot_fieldsMD_" + fieldName; data[matrixName] = & ddot_fieldsMD_[thisField]; matrixName = "dot_dot_fieldsMD_"; data[matrixName] = & dot_dot_fieldsMD_[thisField]; matrixName = "fieldNdOld_" + fieldName; data[matrixName] = & fieldNdOld_[thisField]; matrixName = "fieldNdFiltered_" + fieldName; data[matrixName] = & fieldNdFiltered_[thisField]; matrixName = "fieldRateNdOld_"; data[matrixName] = & fieldRateNdOld_[thisField]; matrixName = "fieldRateNdFiltered_" + fieldName; data[matrixName] = & fieldRateNdFiltered_[thisField]; matrixName = "dot_fieldRateNdOld_" + fieldName; data[matrixName] = & dot_fieldRateNdOld_[thisField]; matrixName = "dot_fieldRateNdFiltered_" + fieldName; data[matrixName] = & dot_fieldRateNdFiltered_[thisField]; matrixName = "auxStorage_" + fieldName; data[matrixName] = & auxStorage_[thisField]; } } //-------------------------------------------------- // write_restart_file // bundle matrices that need to be saved and call // fe_engine to write the file //-------------------------------------------------- void ATC_Transfer::write_restart_data(string fileName, OUTPUT_LIST & data) { pack_fields(data); feEngine_->write_restart_file(fileName,&data); } //-------------------------------------------------- // read_restart_file // bundle matrices that need to be saved and call // fe_engine to write the file //-------------------------------------------------- void ATC_Transfer::read_restart_data(string fileName, OUTPUT_LIST & data) { pack_fields(data); feEngine_->read_restart_file(fileName,&data); } //-------------------------------------------------- // Interactions with LAMMPS fix commands // parse input command and pass on to finite element engine // or physics specific transfers if necessary // revert to physics-specific transfer if no command matches input // first keyword is unique to particular class // base class keyword matching must apply to ALL physics // order: derived, base, owned objects //-------------------------------------------------- bool ATC_Transfer::modify(int narg, char **arg) { FieldName thisField; int thisIndex; int argIdx; bool match = false; if (strcmp(arg[0],"transfer")==0) { /*! \page man_transfer_output fix_modify AtC transfer output \section syntax fix_modify AtC transfer output [text | vector_components | tensor_components ] - filename_prefix (string) = prefix for data files - frequency (integer) = frequency of output in time-steps - options (keyword/s): \n text = creates text output as well as binary Ensight output \n vector_components = outputs vectors as scalar components \n tensor_components = outputs tensor as scalar components (use this for Paraview)\n \section examples fix_modify AtC transfer output heatFE 100 \n fix_modify AtC transfer output hardyFE 1 text tensor_components \n \section description Creates (binary, "gold" format) Ensight output of nodal/mesh data which is transfer/physics specific. \section restrictions \section related see \ref man_fix_atc \section default none */ if (strcmp(arg[1],"output")==0) { if (strcmp(arg[2],"nodeset")==0) { string nset = arg[3]; if (strcmp(arg[4],"sum")==0) { string field = arg[5]; pair < string, FieldName > id(nset,string_to_field(field)); nsetCompute_[id] = 0.0; match = true; } } else { outputPrefix_ = arg[2]; outputFrequency_ = atoi(arg[3]); bool text_output = false, vect_comp = false, tensor_comp = false; int rank = lammpsInterface_->comm_rank(); for (int i = 4; i0) { if (text_output) { feEngine_->initialize_output(rank,outputPrefix_,GNUPLOT); if (rank == 0) cout << " ATC:: Warning : text output can create _LARGE_ files\n"; } else feEngine_->initialize_output(rank,outputPrefix_); if (vect_comp) feEngine_->output_manager() ->set_option(OUTPUT_VECTOR_COMPONENTS,true); if (tensor_comp) feEngine_->output_manager() ->set_option(OUTPUT_TENSOR_COMPONENTS,true); } match = true; } } /*! \page man_transfer_atomic_output fix_modify AtC transfer atomic_output \section syntax fix_modify AtC transfer atomic_output [text] - filename_prefix (string) = prefix for data files - frequency (integer) = frequency of output in time-steps - option "text" = creates a text version of the output as well as the binary Ensight output \section examples fix_modify AtC transfer atomic_output heatMD 100 \n fix_modify AtC transfer atomic_output nanoMD 1 text \n \section description Creates (binary, "gold" format) Ensight output of atom (point) data which is transfer/physics specific. \section restrictions \section related \section default none */ else if (strcmp(arg[1],"atomic_output")==0) { outputPrefixAtom_ = arg[2]; outputFrequencyAtom_ = atoi(arg[3]); if (outputFrequencyAtom_>0) { if (narg == 5 && strcmp(arg[4],"text")==0) { mdOutputManager_.initialize(outputPrefixAtom_,GNUPLOT); cout << " ATC:: Warning : text output can create _LARGE_ files\n"; } else mdOutputManager_.initialize(outputPrefixAtom_); } match = true; } /*! \page man_transfer_internal fix_modify AtC transfer internal \section syntax fix_modify AtC transfer internal type - = type id for non-fixed atoms internal to the FE mesh \section examples fix_modify AtC transfer internal type internal_atoms \section description Command to define the atoms to couple to finite elements. This definition is required for a AtC simulation. \section restrictions \section related see \ref man_transfer_boundary \section default none */ else if ((strcmp(arg[1],"internal")==0) && (strcmp(arg[2],"type")==0)) { int igroup = lammpsInterface_->find_group(arg[3]); groupbit_ |= lammpsInterface_->group_bit(igroup); igroups_.insert(igroup); match = true; } /*! \page man_transfer_boundary fix_modify AtC transfer boundary \section syntax fix_modify AtC transfer boundary type - = type id for atoms that represent a ficticious boundary internal to the FE mesh \section examples fix_modify AtC transfer boundary type ghost_atoms \section description Command to define the atoms that represent the ficticious boundary internal to the FE mesh. For fully overlapped MD/FE domains with periodic boundary conditions no boundary atoms should be defined. \section restrictions \section related see \ref man_transfer_internal \section default none */ else if ((strcmp(arg[1],"boundary")==0) && (strcmp(arg[2],"type")==0)) { int igroup = lammpsInterface_->find_group(arg[3]); groupbitGhost_ |= lammpsInterface_->group_bit(igroup); igroupsGhost_.insert(igroup); match = true; } /*! \page man_internal_quadrature fix_modify AtC transfer internal_quadrature \section syntax fix_modify atc transfer internal_quadrature < on | off > \section examples fix_modify atc transfer internal_quadrature off \section description Command use or not use atomic quadrature on internal elements fully filled with atoms. By turning the internal quadrature off these elements do not contribute to the governing PDE and the fields at the internal nodes follow the weighted averages of the atomic data. \section restrictions \section related \section default on */ else if (strcmp(arg[1],"internal_quadrature")==0) { if (strcmp(arg[2],"on")==0) { atomQuadForInternal_ = true; match = true; } else if (strcmp(arg[2],"off")==0) { atomQuadForInternal_ = false; regionID_ = -1; match = true; } else if (strcmp(arg[3],"off")==0) { for (regionID_ = 0; regionID_ < lammpsInterface_->nregion(); regionID_++) if (strcmp(arg[2],lammpsInterface_->region_name(regionID_)) == 0) break; if (regionID_ < lammpsInterface_->nregion()) { atomQuadForInternal_ = false; match = true; } else regionID_ = -1; } } /*! \page man_initial fix_modify AtC transfer initial \section syntax fix_modify AtC transfer initial - = field name valid for type of physics, temperature | electron_temperature - = name of set of nodes to apply boundary condition - = value or name of function followed by its parameters \section examples fix_modify atc transfer initial temperature groupNAME 10. \section description Sets the initial values for the specified field at the specified nodes. \section restrictions keyword 'all' reserved in nodeset name \section related see \ref man_transfer_internal \section default none */ // set initial conditions else if (strcmp(arg[1],"initial")==0) { argIdx = 2; get_field(arg,argIdx,thisField,thisIndex); string nsetName(arg[argIdx++]); XT_Function * f = NULL; // parse constant if (narg == argIdx+1) { f = XT_Function_Mgr::instance()->get_constant_function(atof(arg[argIdx])); } // parse function else { f = XT_Function_Mgr::instance()->get_function(&(arg[argIdx]),narg-argIdx); } prescribedDataMgr_->fix_initial_field(nsetName,thisField,thisIndex,f); match = true; } /*! \page man_fix_nodes fix_modify AtC transfer fix \section syntax fix_modify AtC transfer fix - = field name valid for type of physics - = name of set of nodes to apply boundary condition - = value or name of function followed by its parameters \section examples fix_modify AtC transfer fix temperature groupNAME 10. \n fix_modify AtC transfer fix temperature groupNAME 0 0 0 10.0 0 0 1.0 \n \section description Creates a constraint on the values of the specified field at specified nodes. \section restrictions keyword 'all' reserved in nodeset name \section related see \ref man_unfix_nodes \section default none */ // fix and unfix nodes else if (strcmp(arg[1],"fix")==0) { argIdx = 2; get_field(arg,argIdx,thisField,thisIndex); string nsetName(arg[argIdx++]); XT_Function * f = NULL; // parse constant if (narg == argIdx+1) { f = XT_Function_Mgr::instance()->get_constant_function(atof(arg[argIdx])); } // parse function else { f = XT_Function_Mgr::instance()->get_function(&(arg[argIdx]),narg-argIdx); } prescribedDataMgr_->fix_field(nsetName,thisField,thisIndex,f); match = true; } /*! \page man_unfix_nodes fix_modify AtC transfer unfix \section syntax fix_modify AtC transfer unfix - = field name valid for type of physics - = name of set of nodes \section examples fix_modify AtC transfer unfix temperature groupNAME \section description Removes constraint on field values for specified nodes. \section restrictions keyword 'all' reserved in nodeset name \section related see \ref man_fix_nodes \section default none */ else if (strcmp(arg[1],"unfix")==0) { argIdx = 2; get_field(arg,argIdx,thisField,thisIndex); string nsetName(arg[argIdx++]); prescribedDataMgr_->unfix_field(nsetName,thisField,thisIndex); match = true; } /*! \page man_source fix_modify AtC transfer source \section syntax fix_modify AtC transfer source - = field name valid for type of physics - = name of set of elements \section examples fix_modify atc transfer source temperature middle temporal_ramp 10. 0. \section description Add domain sources to the mesh. The units are consistent with LAMMPS's units for mass, length and time and are defined by the PDE being solved, e.g. for thermal transfer the balance equation is for energy and source is energy per time. \section restrictions keyword 'all' reserved in element_set name \section related see \ref man_remove_source \section default none */ else if (strcmp(arg[1],"source")==0) { argIdx = 2; get_field(arg,argIdx,thisField,thisIndex); string esetName(arg[argIdx++]); XT_Function * f = NULL; // parse constant if (narg == argIdx+1) { f = XT_Function_Mgr::instance()->get_constant_function(atof(arg[argIdx])); } // parse function else { f = XT_Function_Mgr::instance()->get_function(&(arg[argIdx]),narg-argIdx); } prescribedDataMgr_->fix_source(esetName,thisField,thisIndex,f); fieldMask_(thisField,PRESCRIBED_SOURCE) = true; match = true; } /*! \page man_remove_source fix_modify AtC transfer remove_source \section syntax fix_modify AtC transfer remove_source - = field name valid for type of physics - = name of set of elements \section examples fix_modify atc transfer remove_source temperature groupNAME \section description Remove a domain source. \section restrictions keyword 'all' reserved in element_set name \section related see \ref man_source \section default */ else if (strcmp(arg[1],"remove_source")==0) { argIdx = 2; get_field(arg,argIdx,thisField,thisIndex); string esetName(arg[argIdx++]); prescribedDataMgr_->unfix_source(esetName,thisField,thisIndex); fieldMask_(thisField,PRESCRIBED_SOURCE) = false; match = true; } /*! \page man_fix_flux fix_modify AtC transfer fix_flux \section syntax fix_modify AtC transfer fix_flux - = field name valid for type of physics, temperature | electron_temperature - = name of set of element faces \section examples fix_modify atc transfer fix_flux temperature faceSet 10.0 \n \section description Command for fixing normal fluxes e.g. heat_flux. This command only prescribes the normal component of the physical flux, e.g. heat (energy) flux. The units are in AtC units, i.e. derived from the LAMMPS length, time, and mass scales. \section restrictions Only normal fluxes (Neumann data) can be prescribed. \section related see \ref man_unfix_flux \section default */ else if (strcmp(arg[1],"fix_flux")==0) { argIdx = 2; get_field(arg,argIdx,thisField,thisIndex); string fsetName(arg[argIdx++]); XT_Function * f = NULL; // parse constant if (narg == argIdx+1) { f = XT_Function_Mgr::instance()->get_constant_function(atof(arg[argIdx])); } // parse function else { f = XT_Function_Mgr::instance()->get_function(&(arg[argIdx]),narg-argIdx); } prescribedDataMgr_->fix_flux(fsetName,thisField,thisIndex,f); fieldMask_(thisField,PRESCRIBED_SOURCE) = true; match = true; } /*! \page man_unfix_flux fix_modify AtC transfer unfix_flux \section syntax fix_modify AtC transfer fix_flux - = field name valid for type of physics, temperature | electron_temperature - = name of set of element faces \section examples fix_modify atc transfer unfix_flux temperature faceSet \n \section description Command for removing prescribed normal fluxes e.g. heat_flux, stress. \section restrictions \section related see \ref man_unfix_flux \section default */ else if (strcmp(arg[1],"unfix_flux")==0) { argIdx = 2; get_field(arg,argIdx,thisField,thisIndex); string fsetName(arg[argIdx++]); prescribedDataMgr_->unfix_flux(fsetName,thisField,thisIndex); fieldMask_(thisField,PRESCRIBED_SOURCE) = false; match = true; } // pass off to time filter else if (strcmp(arg[1],"filter")==0) { match = timeFilterManager_.modify(narg-2,&arg[2]); if (match && timeIntegrator_) timeIntegrator_->force_reset(); } } // pass off to extrinsic model else if (strcmp(arg[0],"extrinsic")==0) { match = extrinsicModelManager_.modify(narg-1,&arg[1]); } /*! \page man_atom_element_map fix_modify AtC transfer atom_element_map \section syntax fix_modify AtC transfer atom_element_map \n - frequency (int) : frequency of updating atom-to-continuum maps based on the current configuration - only for eulerian \section examples fix_modify atc transfer atom_element_map eulerian 100 \section description Changes frame of reference from eulerian to lagrangian and sets the frequency for which the map from atoms to elements is reformed and all the attendant data is recalculated. \section restrictions Cannot change map type after initialization. \section related \section default lagrangian */ else if (strcmp(arg[0],"atom_element_map")==0) { if (strcmp(arg[1],"eulerian")==0) { atomToElementMapType_ = EULERIAN; atomToElementMapFrequency_ = atoi(arg[2]); } else { atomToElementMapType_ = LAGRANGIAN; atomToElementMapFrequency_ = 0; } match = true; } /*! \page man_neighbor_reset_frequency fix_modify AtC transfer neighbor_reset_frequency \section syntax fix_modify AtC transfer neighbor_reset_frequency - frequency (int) : reset neighbor map for AtC analysis every "frequency" steps \section examples fix_modify AtC transfer neighbor_reset_frequency 1000 \section description This command allows the control of how often the lists of which atoms pairs are in the support of a node are reformed \section restrictions Only currently relevant to Hardy post-processing see \ref man_fix_atc \section related This command is different from \ref man_atom_element_map in that it only concerns neighbor lists internal to ATC not the map of atom locations to elements which define the shape function values at atoms. \section default The default is for the internal neighbor lists to be recalculated on the next output step (or sample step if time filtering is used). */ else if (strcmp(arg[0],"neighbor_reset_frequency")==0) { neighborResetFrequency_ = atoi(arg[1]); match = true; } /*! \page man_read_restart fix_modify AtC transfer read_restart \section syntax fix_modify AtC transfer read_restart [file_name] \n \section examples fix_modify AtC transfer read_restart ATC_state \n \section description Reads the current state of the fields from a named text-based restart file. \section restrictions The restart file only contains fields and their time derivatives. The reference positions of the atoms and the commands that initialize the fix are not saved e.g. an identical mesh containing the same atoms will have to be recreated. \section related see write_restart \ref man_write_restart \section default none */ else if (strcmp(arg[0],"read_restart")==0) { restartFileName_ = arg[1]; useRestart_ = true; match = true; } /*! \page man_write_restart fix_modify AtC transfer write_restart \section syntax fix_modify AtC transfer write_restart [file_name] \n \section examples fix_modify AtC transfer write_restart restart.mydata \n \section description Dumps the current state of the fields to a named text-based restart file. This done when the command is invoked and not repeated, unlike the similar lammps command. \section restrictions The restart file only contains fields and their time derivatives. The reference positions of the atoms and the commands that initialize the fix are not saved e.g. an identical mesh containing the same atoms will have to be recreated. \section related see read_restart \ref man_read_restart \section default none */ else if (strcmp(arg[0],"write_restart")==0) { string restartFileName(arg[1]); OUTPUT_LIST data; write_restart_data(restartFileName,data); match = true; } // pass off to fe engine if (!match) { match = feEngine_->modify(narg, arg); } return match; } //-------------------------------------------------- // parse_boundary_integration // parses the boundary integration to determine // the type of boundary integration being used //-------------------------------------------------- ATC::ATC_Transfer::BoundaryIntegrationType ATC_Transfer::parse_boundary_integration(int narg, char **arg, const set< pair > * boundaryFaceSet) { int argIndex = 0; BoundaryIntegrationType myBoundaryIntegrationType; if (narg > 0) { if(strcmp(arg[argIndex],"faceset")==0) { argIndex++; myBoundaryIntegrationType = FE_QUADRATURE; string name(arg[argIndex]); boundaryFaceSet = & ( (feEngine_->get_feMesh())->get_faceset(name)); set_boundary_face_set(boundaryFaceSet); } else if (strcmp(arg[argIndex],"interpolate")==0) { myBoundaryIntegrationType = FE_INTERPOLATION; } else if (strcmp(arg[argIndex],"no_boundary")==0) { myBoundaryIntegrationType = NO_QUADRATURE; } else { throw ATC_Error(0,"Bad boundary integration type"); } } else { // default is interpolation myBoundaryIntegrationType = FE_INTERPOLATION; } set_boundary_integration_type(myBoundaryIntegrationType); return myBoundaryIntegrationType; } //-------------------------------------------------- // helper function for parser // handles : "displacement x" and "temperature" by indexing argIdx // for fluxes : only normal fluxes can be prescribed //-------------------------------------------------- void ATC_Transfer::get_field(/*const*/ char ** args, int & argIdx, FieldName & thisField, int & thisIndex) { string thisName = args[argIdx++]; thisField = string_to_field(thisName); map::const_iterator iter = fieldSizes_.find(thisField); if (iter == fieldSizes_.end()) { throw ATC_Error(0,"Bad field name: "+thisName); } string thisDim = args[argIdx]; thisIndex = 0; if (string_to_index(thisDim,thisIndex)) { if ( !( thisIndex < fieldSizes_[thisField]) ) { throw ATC_Error(0,"Bad field dimension "+thisDim); } argIdx++; } } void ATC_Transfer::init_filter() { if (timeFilterManager_.need_reset()) { map::const_iterator field; for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) { FieldName thisField = field->first; int thisSize = field->second; // Allocate restricted fields fieldRateNdFiltered_[thisField].reset(nNodes_,thisSize); fieldRateNdOld_[thisField].reset(nNodes_,thisSize); dot_fieldRateNdFiltered_[thisField].reset(nNodes_,thisSize); dot_fieldRateNdOld_[thisField].reset(nNodes_,thisSize); } } } void ATC_Transfer::set_xref() { double **x = lammpsInterface_->xatom(); for (int i = 0; i < lammpsInterface_->nmax(); i++) { for (int j = 0; j < nsd_; j++) { xref_[i][j] = x[i][j]; } } Xprd_ = lammpsInterface_->domain_xprd(); Yprd_ = lammpsInterface_->domain_yprd(); Zprd_ = lammpsInterface_->domain_zprd(); XY_ = lammpsInterface_->domain_xy(); XZ_ = lammpsInterface_->domain_xz(); YZ_ = lammpsInterface_->domain_yz(); } void ATC_Transfer::initialize() { // initialized_ is set to true by derived class initialize() int i, j; if (!initialized_) { if (!trackCharge_) { fieldSizes_.erase(CHARGE_DENSITY); prescribedDataMgr_->remove_field(CHARGE_DENSITY); for (i = 0; i < NUM_FLUX; i++) fieldMask_(CHARGE_DENSITY,i) = false; } // Initialize finite elements shape function, mapping variables feEngine_->initialize(); nNodes_ = feEngine_->get_nNodes(); nsd_ = feEngine_->get_nsd(); if (nsd_!=lammpsInterface_->dimension()) throw ATC_Error(0,"Spatial dimensions inconsistent between LAMMPS and ATC"); // set reference position of atoms set_xref(); } // determine which elements are strictly in the MD region and mask out for intrinsic fields // start by finding local atom counts reset_nlocal(); int nelts = feEngine_->get_feMesh()->get_nElements(); elementMask_.reset(nelts); if (atomQuadForInternal_) { elementMask_ = true; } else { // determine which elements contain internal atoms Array hasInternal(nelts); Array hasInternal_local(nelts); hasInternal_local = 0; for (int i = 0; i < nLocal_; ++i) { DENS_VEC coords(3); coords(0) = atomicCoords_(0,i); coords(1) = atomicCoords_(1,i); coords(2) = atomicCoords_(2,i); int eltID = feEngine_->get_feMesh()->map_to_element(coords); hasInternal_local(eltID) = 1; } // swap contributions lammpsInterface_->logical_or(hasInternal_local.get_ptr(), hasInternal.get_ptr(),hasInternal.size()); // determine which elements contain ghost atoms Array hasGhost(nelts); Array hasGhost_local(nelts); hasGhost_local = 0; for (int i = 0; i < nLocalGhost_; ++i) { DENS_VEC coords(3); coords(0) = ghostAtomCoords_(0,i); coords(1) = ghostAtomCoords_(1,i); coords(2) = ghostAtomCoords_(2,i); int eltID = feEngine_->get_feMesh()->map_to_element(coords); hasGhost_local(eltID) = 1; } //swap contributions lammpsInterface_->logical_or(hasGhost_local.get_ptr(), hasGhost.get_ptr(),hasGhost.size()); for (int i = 0; i < nelts; ++i) elementMask_(i) = !hasInternal(i) || hasGhost(i); } set & nullElements = feEngine_->null_elements(); set::const_iterator iset; for (iset = nullElements.begin(); iset != nullElements.end(); iset++) { int ielem = *iset; elementMask_(ielem) = false; } // set atomic weights reset_shape_functions(); // creates atom->element map set_atomic_weights(); // group atoms by material int nMatls = physicsModel_->get_nMaterials(); atomMaterialGroups_.reset(nMatls); for (int i = 0; i < nLocal_; i++) { int matId = elementToMaterialMap_(atomToElementMap_(i)); atomMaterialGroups_(matId).insert(i);// not internalToAtom_ } if (atomToElementMapType_ == EULERIAN && nMatls >1 ) { throw ATC_Error(0," multiple materials & eulerian map not supported"); } // compute necessary restructuring matrices maskMat_.reset(nNodes_,nNodes_); if (!useLocalizedLambda_) { for (i = 0; i < nNodes_; ++i) if (shpWeight_(i,i) > 0.) maskMat_(i,i) = 1.; } else { // select nodes for thermostat // 1) their shape functions have both mask and ghost atoms in their support DENS_VEC hasInternalLocal(nNodes_); DENS_VEC hasInternal(nNodes_); if (nLocal_>0) hasInternalLocal = shpFcn_.col_sum(); lammpsInterface_->allsum(hasInternalLocal.get_ptr(),hasInternal.get_ptr(),nNodes_); DENS_VEC hasGhostLocal(nNodes_); DENS_VEC hasGhost(nNodes_); if (nLocalGhost_>0) hasGhostLocal = shpFcnGhost_.col_sum(); lammpsInterface_->allsum(hasGhostLocal.get_ptr(),hasGhost.get_ptr(),nNodes_); for (i = 0; i < nNodes_; ++i) { if (hasInternal(i) > 0. && hasGhost(i) > 0.) { maskMat_(i,i) = 1.; } } // 2) they are in a specified boundary faceset set::const_iterator iter; for (iter = boundaryFaceNames_.begin(); iter != boundaryFaceNames_.end(); iter++) { set nodeSet; feEngine_->get_feMesh()->faceset_to_nodeset(*iter,nodeSet); set::const_iterator nodeIter; for (nodeIter = nodeSet.begin(); nodeIter != nodeSet.end(); nodeIter++) { maskMat_(*nodeIter,*nodeIter) = 1.; } } } // determine non-zero values and row-column indices // also used in creating overlap map and mapping from atoms to elements // note we assume the overlap map is time-independent nodeToOverlapMap_.reset(nNodes_); nNodeOverlap_ = 0; for (i = 0; i < nNodes_; i++) { if (maskMat_(i,i) > 0.) { nodeToOverlapMap_(i) = nNodeOverlap_; nNodeOverlap_++; } else { nodeToOverlapMap_(i) = -1; } } overlapToNodeMap_.reset(nNodeOverlap_); int counter = 0; for (i = 0; i < nNodes_; i++) { if (nodeToOverlapMap_(i)!=-1) { overlapToNodeMap_(counter) = i; counter++; } } // reset local overlap shape function storage reset_NhatOverlap(); // Shape sparse matrix used in Hoover-type thermostats // first get local pattern sized nNodeOverlap X nNodeOverlap //compute_consistent_md_mass_matrix(NhatOverlap_,tmp); DENS_MAT tmpLocal(nNodeOverlap_,nNodeOverlap_); DENS_MAT tmp(nNodeOverlap_,nNodeOverlap_); if (nLocalLambda_>0) { DENS_VEC T(nLocalLambda_); T = 1.0; for (i = 0; i < nNodeOverlap_; i++) { for (j = 0; j < nNodeOverlap_; j++) { for (int iatom = 0; iatom < nLocalLambda_; iatom++) { tmpLocal(i,j) += NhatOverlap_(iatom,i)*T(iatom)*NhatOverlap_(iatom,j); } } } } // second accumulate total pattern lammpsInterface_->allsum(tmpLocal.get_ptr(), tmp.get_ptr(), tmp.size()); // third extract non-zero entries & construct M_T_Template.reset(nNodeOverlap_,nNodeOverlap_); for (i = 0; i < nNodeOverlap_; i++) { for (j = 0; j < nNodeOverlap_; j++) { if (abs(tmp(i,j))>0) { M_T_Template.add(i,j,0.); } } } M_T_Template.compress(); // \int_\Omega N_I dV NodeVolumes_.reset(nNodes_,nNodes_); invNodeVolumes_.reset(nNodes_,nNodes_); feEngine_->compute_lumped_mass_matrix(NodeVolumes_); invNodeVolumes_ = NodeVolumes_.inv(); // this is exact only for uniform meshes and certain types of atomic weights // \int_{\Omega_MD} N_I dV = \sum_\alpha N_I\alpha V_\alpha // DENS_VEC atom_weights_vec(atomicWeights_.size(), // atomicWeights_.get_ptr()); DENS_VEC subsetNodeVolume(NodeVolumes_.size()); DENS_VEC unitVec(atomicWeights_.size()); unitVec = 1.0; restrict_unscaled(unitVec, subsetNodeVolume); fluxMask_.reset(invNodeVolumes_* subsetNodeVolume); DIAG_MAT id(fluxMask_.nRows(),fluxMask_.nCols()); id = 1.0; fluxMaskComplement_ = id + -1.0*fluxMask_; // set flux masks for nodes we can tell by geometry DENS_VEC hasInternalLocal(nNodes_); DENS_VEC hasInternal(nNodes_); if (nLocal_>0) hasInternalLocal = shpFcn_.col_sum(); lammpsInterface_->allsum(hasInternalLocal.get_ptr(),hasInternal.get_ptr(),nNodes_); DENS_VEC hasGhostLocal(nNodes_); DENS_VEC hasGhost(nNodes_); if (nLocalGhost_>0) hasGhostLocal = shpFcnGhost_.col_sum(); lammpsInterface_->allsum(hasGhostLocal.get_ptr(),hasGhost.get_ptr(),nNodes_); for (i = 0; i < nNodes_; ++i) { if (hasInternal(i) > 0.) { if (hasGhost(i) == 0.) { fluxMask_(i,i) = 1.; fluxMaskComplement_(i,i) = 0.; } } else { fluxMask_(i,i) = 0.; fluxMaskComplement_(i,i) = 1.; } } // Initialize fields and set initial conditions if (!initialized_) { nodeType_.reset(nNodes_); for (int iNode = 0; iNode < nNodes_; ++iNode) nodeType_(iNode) = check_shape_function_type(iNode); map::const_iterator field; for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) { FieldName thisField = field->first; int thisSize = field->second; // Allocate fields, initialize to default values, set up initial schedule fields_[thisField].reset(nNodes_,thisSize); dot_fields_[thisField].reset(nNodes_,thisSize); ddot_fields_[thisField].reset(nNodes_,thisSize); dddot_fields_[thisField].reset(nNodes_,thisSize); dot_fieldsMD_[thisField].reset(nNodes_,thisSize); ddot_fieldsMD_[thisField].reset(nNodes_,thisSize); dot_dot_fieldsMD_[thisField].reset(nNodes_,thisSize); // Allocate restricted fields fieldRateNdFiltered_[thisField].reset(nNodes_,thisSize); fieldRateNdOld_[thisField].reset(nNodes_,thisSize); dot_fieldRateNdFiltered_[thisField].reset(nNodes_,thisSize); dot_fieldRateNdOld_[thisField].reset(nNodes_,thisSize); fieldNdFiltered_[thisField].reset(nNodes_,thisSize); fieldNdOld_[thisField].reset(nNodes_,thisSize); // Allocate auxilliary storage auxStorage_[thisField].reset(nNodes_,thisSize); // Dimension finite element rhs matrix rhs_[thisField].reset(nNodes_,thisSize); rhsAtomDomain_[thisField].reset(nNodes_,thisSize); sources_[thisField].reset(nNodes_,thisSize); extrinsicSources_[thisField].reset(nNodes_,thisSize); boundaryFlux_[thisField].reset(nNodes_,thisSize); if (is_intrinsic(thisField)) { // compute FE mass matrix in full domain massMatInv_[thisField].reset(nNodes_,nNodes_); Array massMask(1); massMask(0) = thisField; feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel_, elementToMaterialMap_,massMats_, &elementMask_); // atomic quadrature for FE mass matrix in atomic domain map massMatAq; feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel_, atomMaterialGroups_,atomicWeightsMask_, shpFcnMasked_,massMatAq); // remove contributions from overlap by approximate quadrature // and fully remove contributions from internal massMats_[thisField] -= massMatAq[thisField]; for (int iNode = 0; iNode < nNodes_; iNode++) if (nodeType_(iNode)==MD_ONLY) massMats_[thisField](iNode,iNode) = 0.; // set up mass MD matrices massMatMDInv_[thisField].reset(nNodes_,nNodes_); compute_md_mass_matrix(thisField,massMatsMD_); massMats_[thisField] += massMatsMD_[thisField]; // compute inverse mass matrices since we're using lumped masses for (int iNode = 0; iNode < nNodes_; iNode++) { if (fabs(massMatsMD_[thisField](iNode,iNode))>0) { massMatMDInv_[thisField](iNode,iNode) = 1./massMatsMD_[thisField](iNode,iNode); } if (fabs(massMats_[thisField](iNode,iNode))>0) { massMatInv_[thisField](iNode,iNode) = 1./massMats_[thisField](iNode,iNode); } } } else { // no MD mass matrices needed, regular matrices computed in extrinsic model massMats_[thisField].reset(nNodes_,nNodes_); massMatInv_[thisField].reset(nNodes_,nNodes_); } } // compute inverse consistent mass for the continuity equation SPAR_MAT consistentMass; feEngine_->compute_mass_matrix(consistentMass); consistentMass.compress(); consistentMassInverse_ = inv(consistentMass.dense_copy()); // Apply integration masking and new ICs // initialize schedule derivatives try { set_initial_conditions(); } catch (ATC::ATC_Error& atcError) { if (!useRestart_) throw; } } // set lists of fixed nodes per field map::const_iterator field; for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) { FieldName thisField = field->first; int thisSize = field->second; isFixedNode_[thisField].reset(nNodes_,thisSize); for (int i = 0; i < nNodes_; i++) for (int j = 0; j < thisSize; j++) isFixedNode_[thisField](i,j) = prescribedDataMgr_->is_fixed(i,thisField,j); } return; } void ATC_Transfer::pre_final_integrate() { // atomSwtich is triggered by a lammps call to pre_exchange // (see ATC_Transfer.h) int localSwitch = 0; if (atomSwitch_ || resetSwitch_) localSwitch = 1; lammpsInterface_->int_allmax(&localSwitch,&globalSwitch_); if (globalSwitch_>0) { reset_nlocal(); reset_shape_functions(); set_atomic_weights(); reset_NhatOverlap(); atomSwitch_ = false; globalSwitch_ = 0; } else if (atomToElementMapType_ == EULERIAN && stepCounter_ % atomToElementMapFrequency_ == 0 ) { reset_coordinates(); reset_shape_functions(); set_atomic_weights(); reset_NhatOverlap(); } } void ATC_Transfer::init_integrate_velocity() { int nlocal = lammpsInterface_->nlocal(); int *type = lammpsInterface_->atom_type(); int *mask = lammpsInterface_->atom_mask(); double *mass = lammpsInterface_->atom_mass(); double **v = lammpsInterface_->vatom(); double **f = lammpsInterface_->fatom(); double dtf = 0.5 * lammpsInterface_->dt() * lammpsInterface_->ftm2v(); if (mass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit_) { double dtfm = dtf / mass[type[i]]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; } } } else { double *rmass = lammpsInterface_->atom_rmass(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit_) { double dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; } } } } void ATC_Transfer::init_integrate_position() { int nlocal = lammpsInterface_->nlocal(); double **x = lammpsInterface_->xatom(); double **v = lammpsInterface_->vatom(); double dtv = lammpsInterface_->dt(); int *mask = lammpsInterface_->atom_mask(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit_) { x[i][0] += dtv * v[i][0]; x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; } } } void ATC_Transfer::final_integrate() { int nlocal = lammpsInterface_->nlocal(); int *type = lammpsInterface_->atom_type(); int *mask = lammpsInterface_->atom_mask(); double *mass = lammpsInterface_->atom_mass(); double **v = lammpsInterface_->vatom(); double **f = lammpsInterface_->fatom(); double dtf = 0.5 * lammpsInterface_->dt() * lammpsInterface_->ftm2v(); if (mass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit_) { double dtfm = dtf / mass[type[i]]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; } } } else { double *rmass = lammpsInterface_->atom_rmass(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit_) { double dtfm = dtf / rmass[i]; v[i][0] += dtfm * f[i][0]; v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; } } } } void ATC_Transfer::finish() { // possibly resetting flags and deallocating corresponding memory // FE Engine feEngine_->finish(); return; } //-------------------------------------------------------------- // create_physics_model // - method to create physics model //-------------------------------------------------------------- void ATC_Transfer::create_physics_model(const PhysicsType & physicsType, string matFileName) { if (physicsModel_) { throw ATC_Error(0,"Attempted to create PhysicsModel multiple times in ATC_Transfer"); } // Create PhysicsModel based on physicsType switch (physicsType) { case NO_PHYSICS : break; case THERMAL : physicsModel_ = new PhysicsModelThermal(matFileName, this); break; default: throw ATC_Error(0,"Unknown physics type in ATC_Transfer::create_physics_model"); } } //-------------------------------------------------------------- /** method to trigger construction of mesh data after mesh construction */ //-------------------------------------------------------------- void ATC_Transfer::initialize_mesh_data(void) { int nelts = feEngine_->get_feMesh()->get_nElements(); elementToMaterialMap_.reset(nelts); elementToMaterialMap_ = 0; construct_prescribed_data_manager(); } //-------------------------------------------------------------- /** method to add new fields to the included list */ //-------------------------------------------------------------- void ATC_Transfer::add_fields(map & newFieldSizes) { map::const_iterator field; for (field = newFieldSizes.begin(); field!=newFieldSizes.end(); field++) { FieldName thisField = field->first; int thisSize = field->second; if (fieldSizes_.find(thisField)==fieldSizes_.end()) { fieldSizes_[thisField] = thisSize; } } } }; // namespace ATC