Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71423401
ATC_TransferUtility.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jul 11, 17:09
Size
77 KB
Mime Type
text/x-c
Expires
Sat, Jul 13, 17:09 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
18950279
Attached To
rLAMMPS lammps
ATC_TransferUtility.cpp
View Options
// 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 "LammpsInterface.h"
#include <fstream>
namespace ATC {
//=================================================================
// quadrature weight functions
//=================================================================
//-----------------------------------------------------------------
void ATC_Transfer::set_atomic_weights(void)
{
switch (atomWeightType_) {
case USER:
reset_atomicWeightsGroup();
break;
case LATTICE:
reset_atomicWeightsLattice();
break;
case ELEMENT:
reset_atomicWeightsElement();
break;
case REGION:
if (!initialized_) {
// Uses region bounding box as total volume
// note will *only* work if atoms exactly fill the region
int nregion = lammpsInterface_->nregion();
Array<double> volume(nregion);
for (int i = 0; i < nregion; ++i) {
volume(i) =
(lammpsInterface_->region_xhi(i)-lammpsInterface_->region_xlo(i))
*(lammpsInterface_->region_yhi(i)-lammpsInterface_->region_ylo(i))
*(lammpsInterface_->region_zhi(i)-lammpsInterface_->region_zlo(i));
}
DenseVector<int> localCount(nregion);
DenseVector<int> globalCount(nregion);
// loop over atoms
localCount = 0;
for (int i = 0; i < nLocal_; ++i) {
for (int j = 0; j < nregion; ++j) {
if (lammpsInterface_->region_match(j,
atomicCoords_(0,i),
atomicCoords_(1,i),
atomicCoords_(2,i))) {
localCount(j)++;
}
}
}
// communication to get total atom counts per group
lammpsInterface_->int_allsum(localCount.get_ptr(),
globalCount.get_ptr(),nregion);
for (int i = 0; i < nregion; ++i) {
if (globalCount(i) > 0)
Valpha_[i] = volume(i)/globalCount(i);
else
Valpha_[i] = 0;
}
}
reset_atomicWeightsRegion();
break;
case GROUP:
if (!initialized_) {
// Uses group bounding box as total volume
// note will *only* work if atoms exactly fill the box
int ngroup = lammpsInterface_->ngroup();
int *mask = lammpsInterface_->atom_mask();
double * bounds;
map<int, double> volume;
bounds = new double[6];
for (int i = 0; i < ngroup; ++i) {
lammpsInterface_->group_bounds(i, bounds);
volume[lammpsInterface_->group_bit(i)] =
(bounds[1]-bounds[0])*(bounds[3]-bounds[2])*(bounds[5]-bounds[4]);
}
delete [] bounds;
DenseVector<int> localCount(ngroup);
DenseVector<int> globalCount(ngroup);
// loop over atoms
localCount = 0;
for (int i = 0; i < nLocal_; ++i) {
for (int j = 0; j < ngroup; ++j) {
if (mask[internalToAtom_(i)] & lammpsInterface_->group_bit(j))
localCount(j)++;
}
}
// communication to get total atom counts per group
lammpsInterface_->int_allsum(localCount.get_ptr(),
globalCount.get_ptr(),ngroup);
for (int i = 0; i < ngroup; ++i) {
int iGroupBit = lammpsInterface_->group_bit(i);
if (globalCount(i) > 0)
Valpha_[iGroupBit] =
volume[iGroupBit]/globalCount(i);
else
Valpha_[iGroupBit] = 0;
}
}
reset_atomicWeightsGroup();
break;
case MULTISCALE:
reset_atomicWeightsMultiscale(shpFcn_,atomicWeights_);
break;
}
if( nLocalMask_ > 0 ) {
atomicWeightsMask_.reset(nLocalMask_,nLocalMask_);
int nmask = 0;
for (int i = 0; i < nLocal_; ++i) {
if (atomMask_(i)) {
atomicWeightsMask_(nmask,nmask) = atomicWeights_(i,i);
nmask++;
}
}
}
}
//-----------------------------------------------------------------
// check to see if shape function contains internal or ghost atoms
//-----------------------------------------------------------------
int ATC_Transfer::check_shape_function_type(int nodeIdx)
{
ShapeFunctionType thisShpType;
// NOTE change this check to not use shpWeight
if (shpWeight_(nodeIdx,nodeIdx) > 0.) { // NOTE need tolerance check
DENS_VEC tmpWeights(nNodes_);
DENS_VEC weights(nNodes_);
if (nLocalGhost_>0)
tmpWeights = shpFcnGhost_.col_sum();
lammpsInterface_->allsum(tmpWeights.get_ptr(),weights.get_ptr(),nNodes_);
if (weights(nodeIdx) > 0.) // NOTE need tolerance check
thisShpType = BOUNDARY;
else
thisShpType = MD_ONLY;
}
else
thisShpType = FE_ONLY;
return thisShpType;
}
//-----------------------------------------------------------------
// checks to see if an element is completely inside the MD domain
//-----------------------------------------------------------------
bool ATC_Transfer::check_internal(int eltIdx)
{
// NOTE we use a check only on the bounding boxes, not the actual shapes
int nNode;
DENS_MAT coords;
bool notInMD = false;
feEngine_->element_coordinates(eltIdx,coords);
nNode = coords.nCols();
if (regionID_ < 0) { // loop over groups to find bounds
double xhi,xlo,yhi,ylo,zhi,zlo;
double * tmpBounds = new double[6];
set<int>::iterator iset;
iset = igroups_.begin(); // NOTE inefficient, use some bounds as initial values
lammpsInterface_->group_bounds(*iset,tmpBounds);
xlo = tmpBounds[0];
xhi = tmpBounds[1];
ylo = tmpBounds[2];
yhi = tmpBounds[3];
zlo = tmpBounds[4];
zhi = tmpBounds[5];
for (iset = igroups_.begin(); iset != igroups_.end(); iset++) {
lammpsInterface_->group_bounds(*iset,tmpBounds);
if (xlo>tmpBounds[0]) xlo = tmpBounds[0];
if (xhi<tmpBounds[1]) xhi = tmpBounds[1];
if (ylo>tmpBounds[2]) ylo = tmpBounds[2];
if (yhi<tmpBounds[3]) yhi = tmpBounds[3];
if (zlo>tmpBounds[4]) zlo = tmpBounds[4];
if (zhi<tmpBounds[5]) zhi = tmpBounds[5];
}
bool checkx = true;
bool checky = true;
bool checkz = true;
if (lammpsInterface_->xperiodic() == 1) checkx = false;
if (lammpsInterface_->yperiodic() == 1) checky = false;
if (lammpsInterface_->zperiodic() == 1) checkz = false;
for (int i = 0; i < nNode; ++i) {
if (((coords(0,i)<xlo)||(coords(0,i)>xhi))&&checkx) notInMD = true;
if (((coords(1,i)<ylo)||(coords(1,i)>yhi))&&checky) notInMD = true;
if (((coords(2,i)<zlo)||(coords(2,i)>zhi))&&checkz) notInMD = true;
}
delete [] tmpBounds;
}
else { // use region denoted by regionID_
for (int i = 0; i < nNode; ++i) {
if (!lammpsInterface_->region_match(regionID_,coords(0,i),coords(1,i),coords(2,i))) {
notInMD = true;
}
}
}
return notInMD;
}
//-----------------------------------------------------------------
// checks to see if an element intersects the ghost atoms
//-----------------------------------------------------------------
bool ATC_Transfer::intersect_ghost(int eltIdx)
{
// NOTE we use a check only on the bounding boxes, not the actual shapes
int nNode;
DENS_MAT coords;
bool intersectGhost = false;
feEngine_->element_coordinates(eltIdx,coords);
nNode = coords.nCols();
double xhi,xlo,yhi,ylo,zhi,zlo;
double * tmpBounds;
tmpBounds = new double[6];
set<int>::iterator iset;
iset = igroupsGhost_.begin(); // NOTE inefficient, use some bounds as initial values
lammpsInterface_->group_bounds(*iset,tmpBounds);
xlo = tmpBounds[0];
xhi = tmpBounds[1];
ylo = tmpBounds[2];
yhi = tmpBounds[3];
zlo = tmpBounds[4];
zhi = tmpBounds[5];
for (iset = igroupsGhost_.begin(); iset != igroupsGhost_.end(); iset++) {
lammpsInterface_->group_bounds(*iset,tmpBounds);
if (xlo<tmpBounds[0]) xlo = tmpBounds[0];
if (xhi>tmpBounds[1]) xhi = tmpBounds[1];
if (ylo<tmpBounds[2]) ylo = tmpBounds[2];
if (yhi>tmpBounds[3]) yhi = tmpBounds[3];
if (zlo<tmpBounds[4]) zlo = tmpBounds[4];
if (zhi>tmpBounds[5]) zhi = tmpBounds[5];
}
bool xlower = true;
bool xhigher = true;
bool ylower = true;
bool yhigher = true;
bool zlower = true;
bool zhigher = true;
for (int i = 0; i < nNode; ++i) {
if (coords(0,i)>=xlo)
xlower = false;
if (coords(0,i)<=xhi)
xhigher = false;
if (coords(0,i)>=ylo)
ylower = false;
if (coords(0,i)<=yhi)
yhigher = false;
if (coords(0,i)>=zlo)
zlower = false;
if (coords(0,i)<=zhi)
zhigher = false;
}
if (!(xlower||xhigher||ylower||yhigher||zlower||zhigher))
intersectGhost = true;
delete [] tmpBounds;
return intersectGhost;
}
//=================================================================
// memory management and processor information exchange
//=================================================================
//-----------------------------------------------------------------
// memory usage of local atom-based arrays
//-----------------------------------------------------------------
int ATC_Transfer::memory_usage()
{
int bytes = lammpsInterface_->nmax() * 3 * sizeof(double);
return bytes;
}
//-----------------------------------------------------------------
// allocate local atom-based arrays
//-----------------------------------------------------------------
void ATC_Transfer::grow_arrays(int nmax)
{
xref_ =
lammpsInterface_->grow_2d_double_array(xref_,nmax,3,"fix_atc:xref");
}
//-----------------------------------------------------------------
// copy values within local atom-based arrays
//-----------------------------------------------------------------
void ATC_Transfer::copy_arrays(int i, int j)
{
xref_[j][0] = xref_[i][0];
xref_[j][1] = xref_[i][1];
xref_[j][2] = xref_[i][2];
}
//-----------------------------------------------------------------
// pack values in local atom-based arrays for exchange with another proc
//-----------------------------------------------------------------
int ATC_Transfer::pack_exchange(int i, double *buf)
{
buf[0] = xref_[i][0];
buf[1] = xref_[i][1];
buf[2] = xref_[i][2];
return 3;
}
//-----------------------------------------------------------------
// unpack values in local atom-based arrays from exchange with another proc
//-----------------------------------------------------------------
int ATC_Transfer::unpack_exchange(int nlocal, double *buf)
{
xref_[nlocal][0] = buf[0];
xref_[nlocal][1] = buf[1];
xref_[nlocal][2] = buf[2];
return 3;
}
//-----------------------------------------------------------------
// pack values in local atom-based arrays from exchange with another proc
//-----------------------------------------------------------------
int ATC_Transfer::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = xref_[j][0];
buf[m++] = xref_[j][1];
buf[m++] = xref_[j][2];
}
}
else {
if (lammpsInterface_->domain_triclinic() == 0) {
dx = pbc[0]*Xprd_;
dy = pbc[1]*Yprd_;
dz = pbc[2]*Zprd_;
}
else {
dx = pbc[0]*Xprd_ + pbc[5]*XY_ + pbc[4]*XZ_;
dy = pbc[1]*Yprd_ + pbc[3]*YZ_;
dz = pbc[2]*Zprd_;
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = xref_[j][0] + dx;
buf[m++] = xref_[j][1] + dy;
buf[m++] = xref_[j][2] + dz;
}
}
return 3;
}
//-----------------------------------------------------------------
// unpack values in local atom-based arrays from exchange with another proc
//-----------------------------------------------------------------
void ATC_Transfer::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
xref_[i][0] = buf[m++];
xref_[i][1] = buf[m++];
xref_[i][2] = buf[m++];
}
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
bool ATC_Transfer::get_region_bounds(const char * regionName,
double & xmin, double & xmax,
double & ymin, double & ymax,
double & zmin, double & zmax,
double & xscale,
double & yscale,
double & zscale)
{
int iRegion = lammpsInterface_->get_region_id(regionName);
xscale = lammpsInterface_->region_xscale(iRegion);
yscale = lammpsInterface_->region_yscale(iRegion);
zscale = lammpsInterface_->region_zscale(iRegion);
xmin = lammpsInterface_->region_xlo(iRegion);
xmax = lammpsInterface_->region_xhi(iRegion);
ymin = lammpsInterface_->region_ylo(iRegion);
ymax = lammpsInterface_->region_yhi(iRegion);
zmin = lammpsInterface_->region_zlo(iRegion);
zmax = lammpsInterface_->region_zhi(iRegion);
if (strcmp(lammpsInterface_->region_style(iRegion),"block")==0) {
return true;
}
else {
return false;
}
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::reset_nlocal()
{
nLocalTotal_ = lammpsInterface_->nlocal();
int * mask = lammpsInterface_->atom_mask();
nLocal_ = 0;
nLocalGhost_ = 0;
for (int i = 0; i < nLocalTotal_; ++i) {
if (mask[i] & groupbit_) nLocal_++;
if (mask[i] & groupbitGhost_) nLocalGhost_++;
}
int local_data[2] = {nLocal_, nLocalGhost_};
int data[2] = {0, 0};
lammpsInterface_->int_allsum(local_data,data,2);
nInternal_ = data[0];
nGhost_ = data[1];
// set up internal & ghost maps & coordinates
if (nLocal_>0) {
// set map
internalToAtom_.reset(nLocal_);
int j = 0;
// construct internalToAtom map
// : internal index -> local lammps atom index
for (int i = 0; i < nLocalTotal_; ++i) {
if (mask[i] & groupbit_) internalToAtom_(j++) = i;
}
// construct reverse map
for (int i = 0; i < nLocal_; ++i) {
atomToInternal_[internalToAtom_(i)] = i;
}
atomicCoords_.reset(nsd_,nLocal_);
}
if (nLocalGhost_>0) {
// set map
ghostToAtom_.reset(nLocalGhost_);
int j = 0;
for (int i = 0; i < nLocalTotal_; ++i) {
if (mask[i] & groupbitGhost_) ghostToAtom_(j++) = i;
}
ghostAtomCoords_.reset(nsd_,nLocalGhost_);
}
reset_coordinates();
}
void ATC_Transfer::reset_coordinates()
{
// reset atomic coordinates
double ** xx = xref_;
if (atomToElementMapType_ == EULERIAN) xx = lammpsInterface_->xatom();
// fill coordinates
if (nLocal_ > 0) {
for (int i = 0; i < nLocal_; i++) {
for (int j = 0; j < nsd_; j++)
atomicCoords_(j,i) = xx[internalToAtom_(i)][j];
}
}
// fill ghost coordinates
if (nLocalGhost_ > 0) {
for (int i = 0; i < nLocalGhost_; i++) {
for (int j = 0; j < nsd_; j++)
ghostAtomCoords_(j,i) = xx[ghostToAtom_(i)][j];
}
}
}
// shape functions are based on reference positions
// and only change when atoms cross processor boundaries
// (x_ref migrates with atoms but does not change value)
void ATC_Transfer::reset_shape_functions()
{
// FE_Engine allocates all required memory
// assume initial atomic position is the reference position for now
if (!initialized_) {
shpFcnDerivs_ = vector<SPAR_MAT >(nsd_);
shpFcnDerivsGhost_ = vector<SPAR_MAT >(nsd_);
shpFcnDerivsMask_ = vector<SPAR_MAT >(nsd_);
}
// use FE_Engine to compute shape functions at atomic points
if (nLocal_>0) {
feEngine_->evaluate_shape_functions(atomicCoords_,
shpFcn_,
shpFcnDerivs_,
atomToElementMap_);
}
if (nLocalGhost_>0) {
feEngine_->evaluate_shape_functions(ghostAtomCoords_,
shpFcnGhost_,
shpFcnDerivsGhost_,
ghostAtomToElementMap_);
}
// shpFcnMasked are the shape functions over atoms for which we need to do approximate
// atomic quadrature to remove the FE contributions, i.e. all elements with
// internal atoms if atomQuadForInternal is true and all elements containing
// but not fully filled by atoms (e.g. has internal and ghost types) if
// atomQuadForInternal is false (this has *nothing* to do with thermostats!)
// determine full/partial/empty elements
if (nLocal_>0) {
atomMask_.reset(nLocal_);
if (atomQuadForInternal_) {
atomMask_ = true;
nLocalMask_ = nLocal_;
}
else {
Array<bool> shapeFunctionMask(nNodes_);
shapeFunctionMask = false;
atomMask_ = false;
nLocalMask_ = 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);
atomMask_(i) = elementMask_(eltID);
if (atomMask_(i)) nLocalMask_++;
}
}
if( nLocalMask_ > 0 ) {
atomicCoordsMask_.reset(nsd_,nLocalMask_);
int nmask = 0;
for (int i = 0; i < nLocal_; ++i) {
if (atomMask_(i)) {
for (int j = 0; j < nsd_; ++j) {
atomicCoordsMask_(j,nmask) = atomicCoords_(j,i);
}
nmask++;
}
}
// NOTE inefficient but easy to implement
Array<int> tempArray;
feEngine_->evaluate_shape_functions(atomicCoordsMask_,
shpFcnMasked_,
shpFcnDerivsMask_,
tempArray);
}
}
shpWeight_.reset(nNodes_,nNodes_);
DENS_VEC tmpWeights(nNodes_);
DENS_VEC weights(nNodes_);
// Computed shape function weights
// If weight is function of shape function index only, must be (sum_alpha
// N_I^alpha)^-1
if (nLocal_>0) tmpWeights = shpFcn_.col_sum();
lammpsInterface_->allsum(tmpWeights.get_ptr(),weights.get_ptr(),nNodes_);
for (int i = 0; i < nNodes_; i++) {
if (weights(i) > 0.) {
shpWeight_(i,i) = 1./weights(i);
}
}
return;
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::reset_NhatOverlap()
{
// determine which atoms are being thermostatted, if we use a global lambda
// it is all of them or if we used a local atom it is only those in the support
// of shape functions on the boundary, i.e. those containing atoms in their support
// but not having their support completely in the MD region (e.g. they contain both
// internal and ghost atoms). (this has *nothing* to do with FE integration!)
if (nLocal_ > 0) {
// form mask
Array<bool> atomMask(nLocal_);
if (!useLocalizedLambda_) {
atomMask = true;
nLocalLambda_ = nLocal_;
}
else {
atomMask = false;
nLocalLambda_ = 0;
DENS_VEC coords(nsd_);
for (int i = 0; i < nLocal_; ++i) {
for (int j = 0; j < nsd_; ++j) {
coords(j) = atomicCoords_(j,i);
}
int eltID = feEngine_->get_feMesh()->map_to_element(coords);
Array<int> nodes;
feEngine_->get_feMesh()->element_connectivity_unique(eltID,nodes);
for (int j = 0; j < nodes.size(); ++j) {
if (maskMat_(nodes(j),nodes(j)) > 0.) {
atomMask(i) = true;
nLocalLambda_++;
break;
}
}
}
}
// compute shape functions and normalize
NhatOverlap_.reset(nLocalLambda_,nNodeOverlap_);
if( nLocalLambda_ > 0 ) {
DENS_MAT atomicCoords(nsd_,nLocalLambda_);
int count = 0;
for (int i = 0; i < nLocal_; ++i) {
if (atomMask(i)) {
for (int j = 0; j < nsd_; ++j) {
atomicCoords(j,count) = atomicCoords_(j,i);
}
count++;
}
}
// from shpFcnLambda get nz and map sparsity pattern
SPAR_MAT shpFcnLambda;
Array<int> tempArray;
feEngine_->evaluate_shape_functions(atomicCoords,shpFcnLambda,tempArray);
// NhatOverlap__ = shpFcnLambda*maskMat_
DIAG_MAT scaling(maskMat_*shpWeight_);
CLON_VEC scaleVector(scaling);
NhatOverlapWeights_.reset(NhatOverlap_.nCols(),NhatOverlap_.nCols());
for (int i = 0; i < shpFcnLambda.size(); ++i) {
TRIPLET<double> myTriplet = shpFcnLambda.get_triplet(i);
int myCol = nodeToOverlapMap_(myTriplet.j);
if (myCol > -1) {
NhatOverlap_.set(myTriplet.i,myCol,myTriplet.v);
NhatOverlapWeights_(myCol) = scaleVector(myTriplet.j);
}
}
NhatOverlap_.compress();
// set up matrix that grabs only atoms to be used in the thermostat
Trestrict_.reset(nLocalLambda_,nLocal_);
int j = 0;
for (int i = 0; i < nLocal_; ++i) {
if (atomMask(i)) Trestrict_.set(j++,i,1.0);
}
Trestrict_.compress();
}
}
else {
nLocalLambda_ = 0;
}
// compute overlap ghost shape function
// NOTE if all ghosts do not lie in the support of boundary nodes
// this process will not be correct
if (nLocalGhost_>0) {
shpFcnGhostOverlap_ = shpFcnGhost_*maskMat_;
DENS_VEC activeNodes(nNodes_);
for (int i = 0; i < nNodes_; ++i) {
activeNodes(i) = maskMat_(i,i);
}
DENS_VEC weights = shpFcnGhost_*activeNodes;
DIAG_MAT weightMatrix(nLocalGhost_,nLocalGhost_);
for (int i = 0; i < nLocalGhost_; ++i) {
weightMatrix(i,i) = 1./weights(i);
}
shpFcnGhostOverlap_ = weightMatrix*shpFcnGhostOverlap_;
}
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::reset_atomicWeightsGroup()
{
if (nLocal_>0) {
atomicWeights_.reset(nLocal_,nLocal_);
int *mask = lammpsInterface_->atom_mask();
map<int, double>::const_iterator igroup;
for (igroup = Valpha_.begin(); igroup != Valpha_.end(); igroup++) {
int gid = igroup->first;
double weight = igroup->second;
for (int i = 0; i < nLocal_; ++i) {
if (mask[internalToAtom_(i)] & gid) {
atomicWeights_(i,i) = weight;
}
}
}
}
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::reset_atomicWeightsRegion()
{
if (nLocal_>0) {
atomicWeights_.reset(nLocal_,nLocal_);
map<int, double>::const_iterator iregion;
for (iregion = Valpha_.begin(); iregion != Valpha_.end(); iregion++) {
int rid = iregion->first;
double weight = iregion->second;
for (int i = 0; i < nLocal_; ++i) {
if (lammpsInterface_->region_match(rid,
atomicCoords_(0,i),
atomicCoords_(1,i),
atomicCoords_(2,i)))
atomicWeights_(i,i) = weight;
}
}
}
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::reset_atomicWeightsLattice()
{
// setup quadrature weights uniformly
if (nLocal_>0) {
atomicWeights_.reset(nLocal_,nLocal_);
double volume_per_atom = lammpsInterface_->volume_per_atom();
atomicWeights_ = volume_per_atom;
}
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::reset_atomicWeightsElement()
{
const FE_Mesh * feMesh = feEngine_->get_feMesh();
int nElts = feMesh->get_nElements();
DenseVector<int> eltAtomCtLocal(nElts);
DenseVector<int> eltAtomCt(nElts);
DenseVector<int> atomEltID;
DENS_VEC eltAtomVol(nElts);
for (int i = 0; i < nElts; ++i)
eltAtomCtLocal(i) = 0;
if (nLocal_>0) {
atomicWeights_.reset(nLocal_,nLocal_);
atomEltID.reset(nLocal_);
// determine number of atoms in each element, partial sum
const FE_Mesh * feMesh = feEngine_->get_feMesh();
int nElts = feMesh->get_nElements();
for (int i = 0; i < nLocal_; ++i) {
DENS_VEC coords(3);
coords(0) = xref_[internalToAtom_(i)][0];
coords(1) = xref_[internalToAtom_(i)][1];
coords(2) = xref_[internalToAtom_(i)][2];
int thisElt = feMesh->map_to_element(coords);
eltAtomCtLocal(thisElt) += 1;
atomEltID(i) = thisElt;
}
}
// mpi to determine total atoms
lammpsInterface_->int_allsum(eltAtomCtLocal.get_ptr(),eltAtomCt.get_ptr(),nElts);
// divide element volume by total atoms to get atomic volume
if (nLocal_>0) {
for (int i = 0; i < nElts; ++i) {
double minx, maxx, miny, maxy, minz, maxz;
DENS_MAT nodalCoords;
feMesh->element_coordinates(i,nodalCoords);
minx = nodalCoords(0,0); maxx = nodalCoords(0,0);
miny = nodalCoords(1,0); maxy = nodalCoords(1,0);
minz = nodalCoords(2,0); maxz = nodalCoords(2,0);
for (int j = 1; j < nodalCoords.nCols(); ++j) {
if (nodalCoords(0,j)<minx) minx = nodalCoords(0,j);
if (nodalCoords(0,j)>maxx) maxx = nodalCoords(0,j);
if (nodalCoords(1,j)<miny) miny = nodalCoords(1,j);
if (nodalCoords(1,j)>maxy) maxy = nodalCoords(1,j);
if (nodalCoords(2,j)<minz) minz = nodalCoords(2,j);
if (nodalCoords(2,j)>maxz) maxz = nodalCoords(2,j);
}
double eltVol = (maxx-minx)*(maxy-miny)*(maxz-minz);
if (eltVol<0) eltVol *= -1.;
eltAtomVol(i) = eltVol/eltAtomCt(i);
}
for (int i = 0; i < nLocal_; ++i)
atomicWeights_(i,i) = eltAtomVol(atomEltID(i));
}
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::reset_atomicWeightsMultiscale(const SPAR_MAT & shapeFunctionMatrix,
DIAG_MAT & atomicVolumeMatrix)
{
// solve equation \sum_a N_Ia \sum_J N_Ja dV_J = \int_Omega N_I dV
// form left-hand side
SPAR_MAT lhs;
compute_consistent_md_mass_matrix(shapeFunctionMatrix,lhs);
// form right-hand side
DIAG_MAT rhsMatrix;
feEngine_->compute_lumped_mass_matrix(rhsMatrix);
DENS_VEC rhs;
rhs.copy(rhsMatrix.get_ptr(),rhsMatrix.size(),1);
// change entries for all entries if no atoms in shape function support
double totalVolume = rhs.sum();
double averageVolume = totalVolume/nInternal_;
DENS_VEC scale(nNodes_);
for (int i = 0; i < nNodes_; i++) {
if ((abs(lhs(i,i)) > 0.))
scale(i) = 1.;
else
scale(i) = 0.;
}
lhs.row_scale(scale);
for (int i = 0; i < nNodes_; i++) {
if (scale(i) < 0.5) {
lhs.set(i,i,1.);
rhs(i) = averageVolume;
}
}
lhs.compress();
// solve equation
DIAG_MAT preConditioner = lhs.get_diag();
int myMaxIt = lhs.nRows(); // NOTE could also use a fixed parameter
double myTol = 1.e-10; // NOTE could also use a fixed parameter
DENS_VEC nodalAtomicVolume(nNodes_);
int convergence = CG(lhs, nodalAtomicVolume, rhs, preConditioner, myMaxIt, myTol);
if (convergence>0) // error if didn't converge
throw ATC_Error(0,"CG solver did not converge in ATC_Transfer::reset_atomicWeightsMultiscale");
// compute each atom's volume
DENS_VEC atomicVolume(nLocal_);
prolong(nodalAtomicVolume,atomicVolume);
atomicVolumeMatrix.reset(atomicVolume);
}
//-----------------------------------------------------------------
//
//-----------------------------------------------------------------
void ATC_Transfer::compute_consistent_md_mass_matrix(const SPAR_MAT & shapeFunctionMatrix,
SPAR_MAT & mdMassMatrix)
{
// NOTE could use sparse matrix and do communication broadcasting to exchange triplets
int nRows = shapeFunctionMatrix.nRows();
int nCols = shapeFunctionMatrix.nCols();
DENS_MAT massMatrixLocal(nCols,nCols);
DENS_MAT denseMdMassMatrix(nCols,nCols);
if (nLocal_>0)
massMatrixLocal = shapeFunctionMatrix.transMat(shapeFunctionMatrix);
lammpsInterface_->allsum(massMatrixLocal.get_ptr(),
denseMdMassMatrix.get_ptr(),
denseMdMassMatrix.size());
mdMassMatrix.reset(denseMdMassMatrix,1.e-10);
}
//=================================================================
// Interscale operators
//=================================================================
//--------------------------------------------------------
void ATC_Transfer::restrict(const MATRIX & atomData,
MATRIX & nodeData)
{
// computes nodeData = N*atomData where N are the normalized shape functions
DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols());
if (nLocal_>0) {
workNodeArray = shpWeight_*shpFcn_.transMat(atomData);
}
int count = nodeData.nRows()*nodeData.nCols();
lammpsInterface_->allsum(workNodeArray.get_ptr(),nodeData.get_ptr(),count);
return;
}
//--------------------------------------------------------
void ATC_Transfer::restrict_unscaled(const MATRIX & atomData,
MATRIX & nodeData)
{
// computes nodeData = N*DeltaVAtom*atomData where N are the shape functions
DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols());
if (nLocal_>0) {
workNodeArray = shpFcn_.transMat(atomicWeights_*atomData);
}
int count = nodeData.nRows()*nodeData.nCols();
lammpsInterface_->allsum(workNodeArray.get_ptr(),nodeData.get_ptr(),count);
return;
}
//--------------------------------------------------------
void ATC_Transfer::restrict_volumetric_quantity(const MATRIX & atomData,
MATRIX & nodeData,
const SPAR_MAT & shpFcn)
{
// computes nodeData = N*DeltaVAtom*atomData where N are the shape functions
DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols());
if (nLocal_>0) {
workNodeArray = shpFcn.transMat(atomData);
}
int count = nodeData.nRows()*nodeData.nCols();
lammpsInterface_->allsum(workNodeArray.get_ptr(),nodeData.get_ptr(),count);
return;
}
//--------------------------------------------------------
void ATC_Transfer::restrict_volumetric_quantity(const MATRIX & atomData,
MATRIX & nodeData)
{
restrict_volumetric_quantity(atomData,nodeData,shpFcn_);
return;
}
//--------------------------------------------------------
void ATC_Transfer::project(const MATRIX & atomData,
MATRIX & nodeData,
FieldName thisField)
{
// computes M*nodeData = N*DeltaVAtom*atomData where N are the shape functions
restrict_unscaled(atomData,nodeData);
apply_inverse_mass_matrix(nodeData,thisField);
return;
}
//--------------------------------------------------------
void ATC_Transfer::project_volumetric_quantity(const MATRIX & atomData,
MATRIX & nodeData,
FieldName thisField)
{
project_volumetric_quantity(atomData,nodeData,shpFcn_,thisField);
return;
}
//--------------------------------------------------------
void ATC_Transfer::project_volumetric_quantity(const MATRIX & atomData,
MATRIX & nodeData,
const SPAR_MAT & shpFcn,
FieldName thisField)
{
// computes nodeData = N*atomData where N are the shape functions
restrict_volumetric_quantity(atomData,nodeData,shpFcn);
apply_inverse_mass_matrix(nodeData,thisField);
return;
}
//--------------------------------------------------------
void ATC_Transfer::project_md(const MATRIX & atomData,
MATRIX & nodeData,
FieldName thisField)
{
// computes M*nodeData = N*DeltaVAtom*atomData where N are the shape functions
restrict_unscaled(atomData,nodeData);
apply_inverse_md_mass_matrix(nodeData,thisField);
return;
}
//--------------------------------------------------------
void ATC_Transfer::project_md_volumetric_quantity(const MATRIX & atomData,
MATRIX & nodeData,
FieldName thisField)
{
project_md_volumetric_quantity(atomData,nodeData,shpFcn_,thisField);
return;
}
//--------------------------------------------------------
void ATC_Transfer::project_md_volumetric_quantity(const MATRIX & atomData,
MATRIX & nodeData,
const SPAR_MAT & shpFcn,
FieldName thisField)
{
// computes nodeData = N*atomData where N are the shape functions
restrict_volumetric_quantity(atomData,nodeData,shpFcn);
apply_inverse_md_mass_matrix(nodeData,thisField);
return;
}
//--------------------------------------------------------
void ATC_Transfer::prolong(const MATRIX & nodeData,
MATRIX & atomData)
{
// computes the finite element interpolation at atoms atomData = N*nodeData
if (nLocal_>0) {
atomData = shpFcn_*nodeData;
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::prolong_scaled(const MATRIX & nodeData,
MATRIX & atomData)
{
// computes the finite element interpolation at atoms atomData = N*nodeData
if (nLocal_>0) {
atomData = shpFcn_*(shpWeight_*nodeData);
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::prolong_ghost(const MATRIX & nodeData,
MATRIX & atomData)
{
// computes the finite element interpolation at atoms atomData = N*nodeData
if (nLocalGhost_>0) {
atomData = shpFcnGhost_*nodeData;
}
return;
}
//=================================================================
// Interscale physics constructors
//=================================================================
void ATC_Transfer::compute_atomic_mass(MATRIX & atomicMasses)
{
atomicMasses.reset(nLocal_, 1);
if (nLocal_>0) {
int * type = LammpsInterface::instance()->atom_type();
double * mass = LammpsInterface::instance()->atom_mass();
double * rmass = LammpsInterface::instance()->atom_rmass();
int atomIndex;
double atomMass;
for (int i = 0; i < nLocal_; i++) {
atomIndex = internalToAtom_(i);
if (mass) {
atomMass = mass[type[atomIndex]];
}
else {
atomMass = rmass[atomIndex];
}
atomicMasses[i] = atomMass;
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_charge(MATRIX & atomicCharges)
{
atomicCharges.reset(nLocal_, 1);
if (nLocal_>0) {
double * charge = LammpsInterface::instance()->atom_charge();
//double coef = LammpsInterface::instance()->elemental_charge();
int atomIndex;
for (int i = 0; i < nLocal_; i++) {
atomIndex = internalToAtom_(i);
atomicCharges[i] = charge[atomIndex];
//atomicCharges[i] = coef*charge[atomIndex];
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_temperature(MATRIX & T,
const double * const* v,
double ** v2)
{
// v2 is for fractional step
if (v2==NULL) v2 = const_cast<double**>(v);
// T_alpha = 2/(3kB) * 1/2*m_alpha*(v_alpha \dot v2_alpha)
T.reset(nLocal_, 1);
if (nLocal_>0) {
double ma;
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
double kBoltzmann = lammpsInterface_->kBoltzmann();
int atomIdx;
double coef = 1./(nsd_*kBoltzmann);
for (int i = 0; i < nLocal_; i++) {
atomIdx = internalToAtom_(i);
if (mass) {
ma = mass[type[atomIdx]];
}
else {
ma = rmass[atomIdx];
}
for (int j = 0; j < nsd_; j++) {
T[i] += coef*ma*(v[atomIdx][j]*v2[atomIdx][j]);
}
}
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_kinetic_energy(MATRIX & kineticEnergy,
const double * const* v,
double ** v2)
{
// v2 is for fractional step
if (v2==NULL) v2 = const_cast<double**>(v);
// KE_alpha = 1/2*m_alpha*(v_alpha \dot v2_alpha)
if (nLocal_>0) {
kineticEnergy.reset(nLocal_, 1);
double ma;
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
int atomIdx;
for (int i = 0; i < nLocal_; i++) {
atomIdx = internalToAtom_(i);
if (mass) {
ma = mass[type[atomIdx]];
}
else {
ma = rmass[atomIdx];
}
for (int j = 0; j < nsd_; j++) {
kineticEnergy[i] += ma*(v[atomIdx][j]*v2[atomIdx][j]);
}
}
kineticEnergy *= 0.5;
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_temperature_roc(MATRIX & atomicPower,
const double * const* v,
const double * const* f)
{
// atom_power = 2/(3kB) * v_alpha * f
atomicPower.reset(nLocal_,1);
if (nLocal_>0) {
double coef = 2./(nsd_*lammpsInterface_->boltz());
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++) {
atomicPower[i] += coef*(v[atomIdx][j]*f[atomIdx][j]);
}
}
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_power(MATRIX & atomicPower,
const double * const* v,
const double * const* f)
{
// atom_power = v_alpha * f
atomicPower.reset(nLocal_,1);
if (nLocal_>0) {
double coef = lammpsInterface_->ftm2v();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++) {
atomicPower[i] += coef*(v[atomIdx][j]*f[atomIdx][j]);
}
}
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_temperature_roc(MATRIX & atomicPower,
const double * const* v,
const MATRIX & f)
{
// atom_power = 2/(3kB) * v_alpha * f
atomicPower.reset(nLocal_,1);
if (nLocal_>0) {
double coef = 2./(nsd_*(lammpsInterface_->kBoltzmann()));
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++) {
atomicPower[i] += coef*(v[atomIdx][j]*f(i,j));
}
}
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_power(MATRIX & atomicPower,
const double * const* v,
const MATRIX & f)
{
// atom_power = v_alpha * f
atomicPower.reset(nLocal_,1);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++) {
atomicPower[i] += v[atomIdx][j]*f(i,j);
}
}
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_force_strength(MATRIX & atomicForceStrength,
const double * const* f)
{
// atom_force_strength = 2/(3kB) * f * f / ma
atomicForceStrength.reset(nLocal_,1);
if (nLocal_>0) {
double mvv2e = lammpsInterface_->mvv2e();
double coef = 2./(nsd_*(lammpsInterface_->kBoltzmann())) / mvv2e / mvv2e;
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
double ma;
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
if (mass) {
ma = mass[type[atomIdx]];
}
else {
ma = rmass[atomIdx];
}
for (int j = 0; j < nsd_; j++) {
atomicForceStrength[i] += coef*(f[atomIdx][j]*f[atomIdx][j])/ma;
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_force_strength(MATRIX & atomicForceStrength,
const MATRIX & f)
{
// atom_force_strength = 1/(3kB) * f * f / ma
atomicForceStrength.reset(nLocal_,1);
if (nLocal_>0) {
double coef = 1./(nsd_*(lammpsInterface_->kBoltzmann()));
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
double ma;
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
if (mass) {
ma = mass[type[atomIdx]];
}
else {
ma = rmass[atomIdx];
}
for (int j = 0; j < nsd_; j++) {
atomicForceStrength[i] += coef*(f(i,j)*f(i,j))/ma;
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_force_dot(MATRIX & atomicForceDot,
const double * const* f1,
const MATRIX & f2)
{
// atom_force_strength = 2/(3kB) * f * f / ma
atomicForceDot.reset(nLocal_,1);
if (nLocal_>0) {
double mvv2e = lammpsInterface_->mvv2e();
double coef = 1./(nsd_*(lammpsInterface_->kBoltzmann())) / mvv2e;
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
double ma;
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
if (mass) {
ma = mass[type[atomIdx]];
}
else {
ma = rmass[atomIdx];
}
for (int j = 0; j < nsd_; j++) {
atomicForceDot[i] += coef*(f1[atomIdx][j]*f2(i,j))/ma;
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_lambda_power(MATRIX & lambdaPower,
const MATRIX & lambda,
const double * const* v)
{
// atom_lambda_power = -(1/2) m_a * v_a^2 * lambda
lambdaPower.reset(nLocal_,1);
if (nLocal_>0) {
DENS_VEC lambdaAtom(nLocal_);
prolong_scaled(lambda,lambdaAtom);
double ma;
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
int atomIdx;
double coef = 0.5;
for (int i = 0; i < nLocal_; i++) {
atomIdx = internalToAtom_(i);
if (mass) {
ma = mass[type[atomIdx]];
}
else {
ma = rmass[atomIdx];
}
for (int j = 0; j < nsd_; j++) {
lambdaPower[i] -= coef*ma*v[atomIdx][j]*v[atomIdx][j]*lambdaAtom(i);
}
}
}
return;
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_temperature_roc(MATRIX & atomicPower,
const MATRIX & force,
const double * const* v)
{
// atom_lambda_power = v_{\alpha} \dot force
// NOTE account for missing mass
atomicPower.reset(nLocal_,1);
if (nLocal_>0) {
double ma;
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
int atomIdx;
double coef = 2./(nsd_*(lammpsInterface_->kBoltzmann()));
for (int i = 0; i < nLocal_; i++) {
atomIdx = internalToAtom_(i);
if (mass) {
ma = mass[type[atomIdx]];
}
else {
ma = rmass[atomIdx];
}
for (int j = 0; j < nsd_; j++) {
atomicPower[i] += coef*ma*v[atomIdx][j]*force(i,j);
}
}
}
return;
}
//--------------------------------------------------------
// atomic kinematic quantities
void ATC_Transfer::compute_atomic_displacement(DENS_MAT & atomicDisplacement,
const double * const* x)
{
atomicDisplacement.reset(nLocal_, nsd_);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; ++j) {
atomicDisplacement(i,j) = x[atomIdx][j] - xref_[atomIdx][j];
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_centerOfMass_displacement(DENS_MAT & atomicComD,
const double * const* x)
{
atomicComD.reset(nLocal_, nsd_);
if (nLocal_>0) {
double *mass = lammpsInterface_->atom_mass();
if (mass) {
int *type = lammpsInterface_->atom_type();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicComD(i,j) = mass[type[atomIdx]]*(x[atomIdx][j]- xref_[atomIdx][j]);
}
}
else {
double *rmass = lammpsInterface_->atom_rmass();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicComD(i,j) = rmass[atomIdx]*(x[atomIdx][j]- xref_[atomIdx][j]);
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_position(DENS_MAT & atomicDisplacement,
const double * const* x)
{
atomicDisplacement.reset(nLocal_, nsd_);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; ++j)
atomicDisplacement(i,j) = x[atomIdx][j];
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_centerOfMass(DENS_MAT & atomicCom,
const double * const* x)
{
atomicCom.reset(nLocal_, nsd_);
if (nLocal_>0) {
double *mass = lammpsInterface_->atom_mass();
if (mass) {
int *type = lammpsInterface_->atom_type();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicCom(i,j) = mass[type[atomIdx]]*x[atomIdx][j];
}
}
else {
double *rmass = lammpsInterface_->atom_rmass();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicCom(i,j) = rmass[atomIdx]*x[atomIdx][j];
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_velocity(DENS_MAT & atomicVelocity,
const double * const* v)
{
atomicVelocity.reset(nLocal_, nsd_);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; ++j)
atomicVelocity(i,j) = v[atomIdx][j];
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_momentum(DENS_MAT & atomicMomentum,
const double * const* v)
{
atomicMomentum.reset(nLocal_, nsd_);
if (nLocal_>0) {
double *mass = lammpsInterface_->atom_mass();
if (mass) {
int *type = lammpsInterface_->atom_type();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicMomentum(i,j) = mass[type[atomIdx]]*v[atomIdx][j];
}
}
else {
double *rmass = lammpsInterface_->atom_rmass();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicMomentum(i,j) = rmass[atomIdx]*v[atomIdx][j];
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_acceleration(DENS_MAT & atomicAcceleration,
const double * const* f)
{
atomicAcceleration.reset(nLocal_, nsd_);
if (nLocal_>0) {
double coef = lammpsInterface_->ftm2v();
double *mass = lammpsInterface_->atom_mass();
if (mass) {
int *type = lammpsInterface_->atom_type();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicAcceleration(i,j) = coef*f[atomIdx][j]/mass[type[atomIdx]];
}
}
else {
double *rmass = lammpsInterface_->atom_rmass();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicAcceleration(i,j) = coef*f[atomIdx][j]/rmass[atomIdx];
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::compute_atomic_force(DENS_MAT & atomicForce,
const double * const* f)
{
atomicForce.reset(nLocal_, nsd_);
if (nLocal_>0) {
double coef = lammpsInterface_->ftm2v();
for (int i = 0; i < nLocal_; i++) {
int atomIdx = internalToAtom_(i);
for (int j = 0; j < nsd_; j++)
atomicForce(i,j) = coef*f[atomIdx][j];
}
}
}
//=================================================================
// FE mapping operations
//=================================================================
// Mappings between overlap nodes and unique nodes
void ATC_Transfer::map_unique_to_overlap(const MATRIX & uniqueData,
MATRIX & overlapData)
{
for (int i = 0; i < nNodes_; i++) {
if (nodeToOverlapMap_(i)!=-1) {
for (int j = 0; j < uniqueData.nCols(); j++) {
overlapData(nodeToOverlapMap_(i),j) = uniqueData(i,j);
}
}
}
}
void ATC_Transfer::map_overlap_to_unique(const MATRIX & overlapData,
MATRIX & uniqueData)
{
uniqueData.reset(nNodes_,overlapData.nCols());
for (int i = 0; i < nNodeOverlap_; i++) {
for (int j = 0; j < overlapData.nCols(); j++) {
uniqueData(overlapToNodeMap_(i),j) = overlapData(i,j);
}
}
}
//--------------------------------------------------------
//
//--------------------------------------------------------
void ATC_Transfer::construct_prescribed_data_manager (void) {
prescribedDataMgr_ = new PrescribedDataManager(feEngine_,fieldSizes_);
}
//========================================================
// FE functions
//========================================================
//--------------------------------------------------------
// handle nodesets
//--------------------------------------------------------
void ATC_Transfer::set_fixed_nodes()
{
// set fields
prescribedDataMgr_->set_fixed_fields(simTime_,
fields_,dot_fields_,ddot_fields_,dddot_fields_);
// set related data
map<FieldName,int>::const_iterator field;
for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
FieldName thisField = field->first;
int thisSize = field->second;
for (int inode = 0; inode < nNodes_ ; ++inode) {
for (int thisIndex = 0; thisIndex < thisSize ; ++thisIndex) {
if (prescribedDataMgr_->is_fixed(inode,thisField,thisIndex)) {
dot_fieldsMD_ [thisField](inode,thisIndex) = 0.;
ddot_fieldsMD_ [thisField](inode,thisIndex) = 0.;
dot_dot_fieldsMD_[thisField](inode,thisIndex) = 0.;
rhs_[thisField](inode,thisIndex) = 0.;
}
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::set_initial_conditions()
{
// set fields
prescribedDataMgr_->set_initial_conditions(simTime_,
fields_,dot_fields_,ddot_fields_,dddot_fields_);
// set (all) related data
map<FieldName,int>::const_iterator field;
for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
FieldName thisField = field->first;
int thisSize = field->second;
for (int inode = 0; inode < nNodes_ ; ++inode) {
for (int thisIndex = 0; thisIndex < thisSize ; ++thisIndex) {
dot_fieldsMD_ [thisField](inode,thisIndex) = 0.;
ddot_fieldsMD_ [thisField](inode,thisIndex) = 0.;
dot_dot_fieldsMD_[thisField](inode,thisIndex) = 0.;
rhs_[thisField](inode,thisIndex) = 0.;
}
}
}
}
//--------------------------------------------------------
void ATC_Transfer::set_sources()
{
// set fields
prescribedDataMgr_->set_sources(simTime_,sources_);
}
//--------------------------------------------------------
void ATC_Transfer::output()
{
if (lammpsInterface_->comm_rank() == 0) {
map< pair <string, FieldName>, double >::iterator iter;
for (iter = nsetCompute_.begin(); iter != nsetCompute_.end();iter++){
pair <string, FieldName> id = iter->first;
string nsetName = id.first;
FieldName field = id.second;
double sum = 0.0;
const set<int> nset = feEngine_->get_feMesh()->get_nodeset(nsetName);
set< int >::const_iterator itr;
for (itr = nset.begin(); itr != nset.end();itr++){
int node = *itr;
sum += fields_[field](node,0);
}
iter->second = sum;
string name = nsetName + "_" + field_to_string(field);
feEngine_->add_global(name, sum);
}
}
}
//--------------------------------------------------------
// use part, return to OutputManager maps of data & geometry
// base class atomic output or just data creation
void ATC_Transfer::atomic_output()
{
int me = lammpsInterface_->comm_rank();
int nTotal = (int) lammpsInterface_->natoms();
double ** x = lammpsInterface_->xatom();
double ** v = lammpsInterface_->vatom();
double ** f = lammpsInterface_->fatom();
int * type = lammpsInterface_->atom_type();
double * mass = lammpsInterface_->atom_mass();
double * rmass = lammpsInterface_->atom_rmass();
// map for output data
OUTPUT_LIST atom_data;
// geometry
// NOTE this does not properly handle the output of subset of the entire system
DENS_MAT atomicCoordinatesGlobal(nsd_,nTotal);
DENS_MAT atomicCoordinatesLocal(nsd_,nTotal);
DENS_MAT atomicTypeLocal(nTotal,1);
DENS_MAT atomicType(nTotal,1);
int * tag = lammpsInterface_->atom_tag();
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
//cout << "procID : " << me << ", iLocal_ : " << i << flush;
int atom_index = internalToAtom_(i);
//cout << ", internalToAtom_ : " << atom_index << flush;
int global_index = tag[atom_index]-1;
//cout << ", global_index : " << global_index << flush;
for (int j = 0; j < nsd_; j++) {
atomicCoordinatesLocal(j,global_index) = x[atom_index][j];
//cout << ", x : " << x[atom_index][j] << flush;
}
atomicTypeLocal(global_index,0) = type[atom_index];
//cout << "\n";
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
atomicTypeLocal(global_index,0) = 1.0;
for (int j = 0; j < nsd_; j++) {
atomicCoordinatesLocal(j,global_index) = x[atom_index][j];
}
}
}
lammpsInterface_->allsum(atomicCoordinatesLocal.get_ptr(),
atomicCoordinatesGlobal.get_ptr(), nTotal*nsd_);
if (me==0) mdOutputManager_.write_geometry(atomicCoordinatesGlobal);
lammpsInterface_->allsum(atomicTypeLocal.get_ptr(),
atomicType.get_ptr(), nTotal);
DENS_MAT atomicPositionGlobal(nTotal,nsd_);
if (atomicOutputMask_.count("position") > 0) {
DENS_MAT atomicPositionLocal(nTotal,nsd_);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicPositionLocal(global_index,j) = x[atom_index][j];
}
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicPositionLocal(global_index,j) = x[atom_index][j];
}
}
}
lammpsInterface_->allsum(atomicPositionLocal.get_ptr(),
atomicPositionGlobal.get_ptr(), nTotal*nsd_);
if (me==0) atom_data["atomicPosition"] = & atomicPositionGlobal;
}
// point data, output: displacement, velocity, temperature
if (me==0) atom_data["atomicType"] = & atomicType;
// displacement
// NOTE make this a stand alone function
#define CHEAP_UNWRAP
#ifdef CHEAP_UNWRAP
double box_bounds[3][3];
lammpsInterface_->get_box_bounds(box_bounds[0][0],box_bounds[1][0],
box_bounds[0][1],box_bounds[1][1],
box_bounds[0][2],box_bounds[1][2]);
double box_length[3];
for (int k = 0; k < 3; k++) {
box_length[k] = box_bounds[1][k] - box_bounds[0][k];
}
#endif
bool latticePeriodicity[3];
latticePeriodicity[0] = (bool) lammpsInterface_->xperiodic();
latticePeriodicity[1] = (bool) lammpsInterface_->yperiodic();
latticePeriodicity[2] = (bool) lammpsInterface_->zperiodic();
DENS_MAT atomicDisplacementGlobal(nTotal,nsd_);
if (atomicOutputMask_.count("displacement") > 0) {
DENS_MAT atomicDisplacementLocal(nTotal,nsd_);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicDisplacementLocal(global_index,j)
= x[atom_index][j] -xref_[atom_index][j];
if (latticePeriodicity[j]) {
#ifdef CHEAP_UNWRAP
// the cheap version of unwrapping atomic positions
double u = atomicDisplacementLocal(global_index,j);
if (u >= 0.5*box_length[j]) { u -= box_length[j]; }
if (u <= -0.5*box_length[j]) { u += box_length[j]; }
atomicDisplacementLocal(global_index,j) = u;
#else
double x_a[3];
lammpsInterface_->unwrap_coordinates(atom_index,x_a);
atomicDisplacementLocal(global_index,j)
= x_a[j] -xref_[atom_index][j];
#endif
}
}
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicDisplacementLocal(global_index,j)
= x[atom_index][j] -xref_[atom_index][j];
if (latticePeriodicity[j]) {
#ifdef CHEAP_UNWRAP
// the cheap version of unwrapping atomic positions
double u = atomicDisplacementLocal(global_index,j);
if (u >= 0.5*box_length[j]) { u -= box_length[j]; }
if (u <= -0.5*box_length[j]) { u += box_length[j]; }
atomicDisplacementLocal(global_index,j) = u;
#else
double x_a[3];
lammpsInterface_->unwrap_coordinates(atom_index,x_a);
atomicDisplacementLocal(global_index,j)
= x_a[j] -xref_[atom_index][j];
#endif
}
}
}
}
lammpsInterface_->allsum(atomicDisplacementLocal.get_ptr(),
atomicDisplacementGlobal.get_ptr(), nTotal*nsd_);
if (me==0) atom_data["atomicDisplacement"] = & atomicDisplacementGlobal;
}
// velocity
DENS_MAT atomicVelocityGlobal(nTotal,nsd_);
if (atomicOutputMask_.count("velocity") > 0) {
DENS_MAT atomicVelocityLocal(nTotal,nsd_);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicVelocityLocal(global_index,j) = v[atom_index][j];
}
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicVelocityLocal(global_index,j) = v[atom_index][j];
}
}
}
lammpsInterface_->allsum(atomicVelocityLocal.get_ptr(),
atomicVelocityGlobal.get_ptr(), nTotal*nsd_);
if (me==0) atom_data["atomicVelocity"] = & atomicVelocityGlobal;
}
// temperature
DENS_MAT atomicTemperatureGlobal(nTotal,1);
if (atomicOutputMask_.count("temperature") > 0) {
DENS_MAT atomicTemperatureLocal(nTotal,1);
atomicTemperatureLocal = 0.0;
double ma;
double coef = 1./(nsd_*(lammpsInterface_->kBoltzmann()));
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
if (mass) ma = mass[type[atom_index]];
else ma = rmass[atom_index];
for (int j = 0; j < nsd_; j++) {
double vj = v[atom_index][j];
atomicTemperatureLocal(global_index,0) += coef*ma*vj*vj;
}
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
if (mass) ma = mass[type[atom_index]];
else ma = rmass[atom_index];
for (int j = 0; j < nsd_; j++) {
double vj = v[atom_index][j];
atomicTemperatureLocal(global_index,0) += coef*ma*vj*vj;
}
}
}
lammpsInterface_->allsum(atomicTemperatureLocal.get_ptr(),
atomicTemperatureGlobal.get_ptr(), nTotal);
if (me==0) atom_data["atomicTemperature"] = & atomicTemperatureGlobal;
}
// force
DENS_MAT atomicForceGlobal(nTotal,nsd_);
if (atomicOutputMask_.count("force") > 0) {
DENS_MAT atomicForceLocal(nTotal,nsd_);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicForceLocal(global_index,j) = f[atom_index][j];
}
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicForceLocal(global_index,j) = f[atom_index][j];
}
}
}
lammpsInterface_->allsum(atomicForceLocal.get_ptr(),
atomicForceGlobal.get_ptr(), nTotal*nsd_);
if (me==0) atom_data["atomicForce"] = & atomicVelocityGlobal;
}
// atomic stress
double ** atomStress = NULL;
try{
atomStress = lammpsInterface_->compute_vector_data("stress/atom");
}
catch(ATC::ATC_Error& atcError) {
atomStress = NULL;
}
// if compute stress/atom isn't create, this code is skipped
DENS_MAT atomicVirialGlobal(nTotal,6);
if (atomicOutputMask_.count("virial") > 0 && atomStress) {
DENS_MAT atomicVirialLocal(nTotal,6);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
for (int k = 0; k < 6; k++) {
atomicVirialLocal(global_index,k) = atomStress[atom_index][k];
}
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
for (int k = 0; k < 6; k++) {
atomicVirialLocal(global_index,k) = atomStress[atom_index][k];
}
}
}
lammpsInterface_->allsum(atomicVirialLocal.get_ptr(),
atomicVirialGlobal.get_ptr(), nTotal*6);
//atomicVirialGlobal *= nktv2p;
if (me==0) atom_data["atomicVirial"] = & atomicVirialGlobal;
}
// potential energy
// if compute pe/atom isn't created, this code is skipped
double * atomPE = NULL;
try{
atomPE = lammpsInterface_->compute_scalar_data("pe/atom");
}
catch(ATC::ATC_Error& atcError) {
atomPE = NULL;
}
DENS_MAT atomicEnergyGlobal(nTotal,1);
if (atomicOutputMask_.count("potential_energy") > 0 && atomPE) {
DENS_MAT atomicEnergyLocal(nTotal,1);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
atomicEnergyLocal(global_index,0) = atomPE[atom_index];
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
atomicEnergyLocal(global_index,0) = atomPE[atom_index];
}
}
lammpsInterface_->allsum(atomicEnergyLocal.get_ptr(),
atomicEnergyGlobal.get_ptr(), nTotal);
if (me==0) atom_data["atomicPotentialEnergy"] = & atomicEnergyGlobal;
}
// centrosymmetry
// NOTE if hardy doesn't create compute centro/atom this is skipped
// we should move this to base class with a request from derived
//double * atomCentro = lammpsInterface_->atomCentro_compute();
double * atomCentro = NULL;
try{
atomCentro = lammpsInterface_->compute_scalar_data("centro/atom");
}
catch(ATC::ATC_Error& atcError) {
atomCentro = NULL;
}
// if compute centro/atom isn't create, this code is skipped
DENS_MAT atomicCentroGlobal(nTotal,1);
if (atomicOutputMask_.count("centrosymmetry") > 0 && atomCentro) {
DENS_MAT atomicCentroLocal(nTotal,1);
if (nLocal_>0) {
for (int i = 0; i < nLocal_; ++i) {
int atom_index = internalToAtom_(i);
int global_index = tag[atom_index]-1;
for (int j = 0; j < nsd_; j++) {
atomicCentroLocal(global_index,0) = atomCentro[atom_index];
}
}
}
if (nLocalGhost_>0) {
for (int i = 0; i < nLocalGhost_; ++i) {
int atom_index = ghostToAtom_(i);
int global_index = tag[atom_index]-1;
atomicCentroLocal(global_index,0) = atomCentro[atom_index];
}
}
lammpsInterface_->allsum(atomicCentroLocal.get_ptr(),
atomicCentroGlobal.get_ptr(), nTotal);
if (me==0) atom_data["atomicCentrosymmetry"] = & atomicCentroGlobal;
}
if (me==0) mdOutputManager_.write_data(simTime_, & atom_data);
}
//=================================================================
// FE interfaces
//=================================================================
void ATC_Transfer::compute_boundary_flux(const Array2D<bool> & rhsMask,
const FIELDS & fields,
FIELDS & rhs)
{
if (bndyIntType_ == FE_QUADRATURE) {
//cout << "FE QUAD " << lammpsInterface_->comm_rank() << endl;
feEngine_->compute_boundary_flux(rhsMask,
fields,
physicsModel_,
elementToMaterialMap_,
(* bndyFaceSet_),
rhs);
}
else if (bndyIntType_ == FE_INTERPOLATION) {
//cout << "FE INTERP " << lammpsInterface_->comm_rank() << endl;
feEngine_->compute_boundary_flux(rhsMask,
fields,
physicsModel_,
elementToMaterialMap_,
atomMaterialGroups_,
atomicWeights_,
shpFcn_,
shpFcnDerivs_,
fluxMask_,
rhs);
}
else if (bndyIntType_ == NO_QUADRATURE) {
//cout << "RESET BOUNDARY FLUX " << lammpsInterface_->comm_rank() << endl;
FIELDS::const_iterator field;
for (field = fields.begin(); field != fields.end(); field++) {
FieldName thisFieldName = field->first;
if (rhsMask(thisFieldName,FLUX)) {
int nrows = (field->second).nRows();
int ncols = (field->second).nCols();
rhs[thisFieldName].reset(nrows,ncols);
}
}
}
}
//-----------------------------------------------------------------
void ATC_Transfer::compute_flux(const Array2D<bool> & rhsMask,
const FIELDS & fields,
GRAD_FIELDS & flux,
const PhysicsModel * physicsModel)
{
if (! physicsModel) { physicsModel = physicsModel_; }
feEngine_->compute_flux(rhsMask,
fields,
physicsModel,
elementToMaterialMap_,
flux);
}
//-----------------------------------------------------------------
void ATC_Transfer::evaluate_rhs_integral(const Array2D<bool> & rhsMask,
const FIELDS & fields, FIELDS & rhs,
const IntegrationDomainType domain,
const PhysicsModel * physicsModel)
{
if (!physicsModel) physicsModel = physicsModel_;
// compute B and N flux NOTE linear or not
if (domain == FE_DOMAIN ) {
feEngine_->compute_rhs_vector(rhsMask,
fields,
physicsModel,
elementToMaterialMap_,
rhs,
&elementMask_);
if (nLocalMask_ > 0) {
feEngine_->compute_rhs_vector(rhsMask,
fields,
physicsModel,
atomMaterialGroups_,
atomicWeightsMask_,
shpFcnMasked_,
shpFcnDerivsMask_,
rhsAtomDomain_);
for (FIELDS::const_iterator field = fields.begin();
field != fields.end(); field++) {
FieldName thisFieldName = field->first;
rhs[thisFieldName] -= rhsAtomDomain_[thisFieldName];
}
}
}
else if (domain == ATOM_DOMAIN) {
// NOTE should we let this actually do a full atom quadrature integration?
if (nLocalMask_ > 0) {
feEngine_->compute_rhs_vector(rhsMask,
fields,
physicsModel,
atomMaterialGroups_,
atomicWeightsMask_,
shpFcnMasked_,
shpFcnDerivsMask_,
rhs);
}
}
else { // domain == FULL_DOMAIN
feEngine_->compute_rhs_vector(rhsMask,
fields,
physicsModel,
elementToMaterialMap_,
rhs);
}
}
//-----------------------------------------------------------------
void ATC_Transfer::compute_rhs_vector(const Array2D<bool> & rhsMask,
const FIELDS & fields, FIELDS & rhs,
const IntegrationDomainType domain,
const PhysicsModel * physicsModel)
{
if (!physicsModel) physicsModel = physicsModel_;
// NOTE what about FE_DOMAIN & Extrinsic???
int index;
// compute FE contributions
evaluate_rhs_integral(rhsMask,fields,rhs,domain,physicsModel);
// NOTE change massMask to rhsMask
Array<FieldName> massMask;
for (FIELDS::const_iterator field = fields.begin();
field != fields.end(); field++)
{
FieldName thisFieldName = field->first;
if (rhsMask(thisFieldName,PRESCRIBED_SOURCE)) {
if (is_intrinsic(thisFieldName)) {
rhs[thisFieldName] += fluxMaskComplement_*sources_[thisFieldName];
} // NOTE perhaps replace fluxMask in future
else {
rhs[thisFieldName] += sources_[thisFieldName];
}
}
// add in sources from extrinsic models
if (rhsMask(thisFieldName,EXTRINSIC_SOURCE)) {
if (is_intrinsic(thisFieldName)) {
rhs[thisFieldName] += fluxMaskComplement_*extrinsicSources_[thisFieldName];
}
else {
rhs[thisFieldName] += extrinsicSources_[thisFieldName];
}
}
}
}
//-----------------------------------------------------------------
void ATC_Transfer::compute_atomic_sources(const Array2D<bool> & fieldMask,
const FIELDS & fields,
FIELDS & atomicSources)
{
// NOTE change massMask to rhsMask
Array<FieldName> massMask;
for (FIELDS::const_iterator field = fields.begin();
field != fields.end(); field++) {
FieldName thisFieldName = field->first;
if (is_intrinsic(thisFieldName)) {
atomicSources[thisFieldName] = 0.;
if (fieldMask(thisFieldName,FLUX)) {
atomicSources[thisFieldName] = boundaryFlux_[thisFieldName];
}
if (fieldMask(thisFieldName,PRESCRIBED_SOURCE)) {
atomicSources[thisFieldName] -= fluxMask_*sources_[thisFieldName];
} // NOTE perhaps replace fluxMask in future
// add in sources from extrinsic models
if (fieldMask(thisFieldName,EXTRINSIC_SOURCE))
atomicSources[thisFieldName] -= fluxMask_*extrinsicSources_[thisFieldName];
}
}
}
void ATC_Transfer::update_filter(MATRIX & filteredQuantity,
const MATRIX & unfilteredQuantity,
MATRIX & unfilteredQuantityOld,
const double dt)
{
double tau = timeFilterManager_.get_filter_scale();
filteredQuantity = 1./(1./dt+1./(2*tau))*( 1./(2*tau)*
(unfilteredQuantity+unfilteredQuantityOld) +
(1./dt-1./(2*tau))*filteredQuantity);
unfilteredQuantityOld = unfilteredQuantity;
return;
}
//--------------------------------------------------------
void ATC_Transfer::update_filter_implicit(MATRIX & filteredQuantity,
const MATRIX & unfilteredQuantity,
const double dt)
{
double tau = timeFilterManager_.get_filter_scale();
filteredQuantity = (1./(1.+dt/tau))*((dt/tau)*unfilteredQuantity + filteredQuantity);
return;
}
};
Event Timeline
Log In to Comment