Homec4science
Diffusion lammps (lammps-icms) Browse

lammps/lib/atc/DislocationExtractor.h
db1e17197fcflammps-icms

Last Change
Raw File

DislocationExtractor.h
View Options

///////////////////////////////////////////////////////////////////////////////
//
// Copyright 2010, Alexander Stukowski
// All rights reserved. See README.txt for more information.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef __DXA_DISLOCATION_EXTRACTOR_H
#define __DXA_DISLOCATION_EXTRACTOR_H
// Include the LAMMPS headers that we need.
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "domain.h"
#include "update.h"
#include "neigh_list.h"
using namespace LAMMPS_NS;
// This is must be defined when we build the LAMMPS module (and not the DXA standalone tool).
#define DXA_LAMMPS_MODULE
// Enables additional sanity checks.
#define DEBUG_DISLOCATIONS
// In LAMMPS, we always use double precision floating-point numbers.
#define USE_DOUBLE_PRECISION_FP
#include "../src/dxa/DXATracing.h"
#include "../src/util/FILEStream.h"
/**
*
*/
class DXADislocationExtractor : private DXATracing, Pointers
{
public:
/// Constructor.
///
/// Expects a pointer to the global LAMMPS object.
///
/// Exact mode flag:
///
/// Enables serial processing, i.e. the master processor performs the complete analysis.
/// This allows for proper handling of periodic boundary conditions.
DXADislocationExtractor(LAMMPS* lmp, bool _exactMode = false) : DXATracing(msgLoggerStream, verboseLoggerStream), Pointers(lmp)
{
exactMode = _exactMode;
// Connect internal logging streams to global LAMMPS streams.
if(comm->me == 0) {
msgLoggerStream.rdbuf()->setFiles(screen, logfile);
verboseLoggerStream.rdbuf()->setFiles(screen, logfile);
}
this->processor = comm->me;
}
/// Extracts all dislocations from the current simulation state.
///
/// First copies atoms from LAMMPS array to internal memory.
/// Also adopts the nearest neighbor lists from LAMMPS.
///
/// Parameters:
///
/// neighborList A full neighbor list provided by LAMMPS, which contains neighbors of ghost atoms.
/// cnaCutoff Cutoff radius to be used for the common neighbor analysis.
///
void extractDislocations(NeighList* neighborList, double cnaCutoff, int lineSmoothingLevel = DEFAULT_LINE_SMOOTHING_LEVEL, int lineCoarseningLevel = DEFAULT_LINE_COARSENING_LEVEL)
{
// Check parameters.
if(cnaCutoff <= 0.0) lammps->error->all(FLERR,"Common neighbor analysis cutoff radius must be positive.");
// Release any memory allocated during last call.
cleanup();
// Setup simulation cell.
if(exactMode) {
pbc[0] = domain->periodicity[0] != 0;
pbc[1] = domain->periodicity[1] != 0;
pbc[2] = domain->periodicity[2] != 0;
}
else {
pbc[0] = pbc[1] = pbc[2] = false;
}
simulationCell(0,0) = domain->h[0];
simulationCell(1,1) = domain->h[1];
simulationCell(2,2) = domain->h[2];
simulationCell(1,2) = domain->h[3];
simulationCell(0,2) = domain->h[4];
simulationCell(0,1) = domain->h[5];
simulationCell(1,0) = simulationCell(2,0) = simulationCell(2,1) = 0;
simulationCellOrigin.X = domain->boxlo[0];
simulationCellOrigin.Y = domain->boxlo[1];
simulationCellOrigin.Z = domain->boxlo[2];
setCNACutoff(cnaCutoff);
setupSimulationCell(cnaCutoff);
timestep = update->ntimestep;
// Transfer LAMMPS atoms to internal array.
if(exactMode)
transferAllLAMMPSAtoms();
else
transferLocalLAMMPSAtoms(neighborList);
if(exactMode == false || processor == 0) {
// Perform common neighbor analysis to identify crystalline atoms.
performCNA();
// Order the neighbor lists of crystalline atoms.
orderCrystallineAtoms();
// Cluster crystalline atoms.
clusterAtoms();
#ifdef PROCESSED_ATOMS
// Enable this code block to get the processed atoms.
{
ostringstream ss;
ss << "proc_" << comm->me << ".atoms.dump";
ofstream stream(ss.str().c_str());
writeAtomsDumpFile(stream);
}
#endif
// Create the nodes of the interface mesh.
createInterfaceMeshNodes();
// Create the facets of the interface mesh.
createInterfaceMeshFacets();
#ifdef DEBUG_DISLOCATIONS
// Check the generated mesh.
validateInterfaceMesh();
#endif
// Generate and advance Burgers circuits on the interface mesh.
traceDislocationSegments();
// Smooth dislocation lines.
smoothDislocationSegments(lineSmoothingLevel, lineCoarseningLevel);
}
if(exactMode) {
// Wrap segments at simulation cell boundaries.
wrapDislocationSegments();
// Distribute extracted dislocation segments to back all processors.
broadcastDislocations();
}
// Clip dislocation lines at processor domain.
Point3 subOrigin;
Matrix3 subCell;
if(domain->triclinic) {
subOrigin = simulationCellOrigin + simulationCell * Vector3(domain->sublo_lamda);
subCell.setColumn(0, simulationCell.column(0) * (domain->subhi_lamda[0] - domain->sublo_lamda[0]));
subCell.setColumn(1, simulationCell.column(1) * (domain->subhi_lamda[1] - domain->sublo_lamda[1]));
subCell.setColumn(2, simulationCell.column(2) * (domain->subhi_lamda[2] - domain->sublo_lamda[2]));
} else {
subOrigin = Point3(domain->sublo);
subCell.setColumn(0, Vector3(domain->subhi[0] - domain->sublo[0], 0, 0));
subCell.setColumn(1, Vector3(0, domain->subhi[1] - domain->sublo[1], 0));
subCell.setColumn(2, Vector3(0, 0, domain->subhi[2] - domain->sublo[2]));
}
clipDislocationLines(subOrigin, subCell);
#ifdef OUTPUT_SEGMENTS
// Enable this code block to see the extracted dislocation segments on each processor.
{
ostringstream ss;
ss << "proc_" << comm->me << ".dislocations.vtk";
ofstream stream(ss.str().c_str());
writeDislocationsVTKFile(stream);
}
#endif
}
/// Returns the list of extracted dislocation segments.
const vector<DislocationSegment*>& getSegments() const { return this->segments; }
protected:
/// Copies local atoms and neighbor lists from LAMMPS to internal array.
void transferLocalLAMMPSAtoms(NeighList* neighborList)
{
// Allocate internal atoms array.
int numTotalAtoms = atom->nlocal + atom->nghost;
inputAtoms.resize(numTotalAtoms);
numLocalInputAtoms = numTotalAtoms;
// Transfer local and ghost atoms to internal array.
for(int i = 0; i < numTotalAtoms; i++) {
InputAtom& a = inputAtoms[i];
a.tag = i + 1;
a.pos.X = atom->x[i][0];
a.pos.Y = atom->x[i][1];
a.pos.Z = atom->x[i][2];
a.flags = 0;
a.cluster = NULL;
a.numNeighbors = 0;
a.setFlag(ATOM_IS_LOCAL_ATOM);
}
// Adopt nearest neighbor lists from LAMMPS.
if(neighborList->ghostflag) {
DISLOCATIONS_ASSERT(neighborList->inum + neighborList->gnum == numTotalAtoms);
int tt = 0;
for(int ii = 0; ii < numTotalAtoms; ii++) {
int i = neighborList->ilist[ii];
DISLOCATIONS_ASSERT(i < numTotalAtoms);
InputAtom& ai = inputAtoms[i];
int* jlist = neighborList->firstneigh[i];
int jnum = neighborList->numneigh[i];
if(ii >= neighborList->inum)
tt += jnum;
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
DISLOCATIONS_ASSERT(j >= 0 && j < numTotalAtoms);
if(j < i) continue;
InputAtom& aj = inputAtoms[j];
if(DistanceSquared(aj.pos, ai.pos) >= cnaCutoff * cnaCutoff) continue;
if(ai.numNeighbors == MAX_ATOM_NEIGHBORS)
raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", ai.tag, MAX_ATOM_NEIGHBORS);
ai.addNeighbor(&aj);
if(aj.numNeighbors == MAX_ATOM_NEIGHBORS)
raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", aj.tag, MAX_ATOM_NEIGHBORS);
aj.addNeighbor(&ai);
}
}
}
else {
for(int ii = 0; ii < neighborList->inum; ii++) {
int i = neighborList->ilist[ii];
InputAtom& ai = inputAtoms[i];
int* jlist = neighborList->firstneigh[i];
int jnum = neighborList->numneigh[i];
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
DISLOCATIONS_ASSERT(j >= 0 && j < numTotalAtoms);
if(j < i) continue;
InputAtom& aj = inputAtoms[j];
if(DistanceSquared(aj.pos, ai.pos) < cnaCutoff * cnaCutoff) {
if(ai.numNeighbors == MAX_ATOM_NEIGHBORS)
raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", ai.tag, MAX_ATOM_NEIGHBORS);
ai.addNeighbor(&aj);
if(aj.numNeighbors == MAX_ATOM_NEIGHBORS)
raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", aj.tag, MAX_ATOM_NEIGHBORS);
aj.addNeighbor(&ai);
}
if(j >= atom->nlocal) {
for(int kk = 0; kk < jnum; kk++) {
int k = jlist[kk];
k &= NEIGHMASK;
DISLOCATIONS_ASSERT(k < (int)inputAtoms.size());
if(k == j) continue;
if(k < atom->nlocal) continue;
InputAtom& ak = inputAtoms[k];
if(DistanceSquared(aj.pos, ak.pos) >= cnaCutoff * cnaCutoff) continue;
if(aj.hasNeighbor(&ak)) continue;
if(aj.numNeighbors == MAX_ATOM_NEIGHBORS)
raiseError("Maximum number of nearest neighbors exceeded. Atom %i has more than %i nearest neighbors (built-in maximum number).", aj.tag, MAX_ATOM_NEIGHBORS);
aj.addNeighbor(&ak);
}
}
}
}
}
}
/// Copies atoms from LAMMPS to internal array and gathers all atoms on the master processor.
void transferAllLAMMPSAtoms()
{
// Determine maximum number of local atoms a single processor has.
int nlocalatoms_max;
MPI_Reduce(&atom->nlocal, &nlocalatoms_max, 1, MPI_INT, MPI_MAX, 0, world);
if(processor == 0) {
// Allocate internal atoms array.
inputAtoms.resize(atom->natoms);
numLocalInputAtoms = atom->natoms;
vector<InputAtom>::iterator currentAtom = inputAtoms.begin();
// Copy local atoms of master processor into internal array.
for(int i = 0; i < atom->nlocal; i++, ++currentAtom) {
currentAtom->tag = currentAtom - inputAtoms.begin() + 1;
currentAtom->pos.X = atom->x[i][0];
currentAtom->pos.Y = atom->x[i][1];
currentAtom->pos.Z = atom->x[i][2];
currentAtom->flags = 0;
currentAtom->cluster = NULL;
currentAtom->numNeighbors = 0;
currentAtom->setFlag(ATOM_IS_LOCAL_ATOM);
}
if(comm->nprocs > 1) {
// Allocate receive buffer.
vector<double> buffer(nlocalatoms_max * 3);
// Receive atoms from other processors.
for(int iproc = 1; iproc < comm->nprocs; iproc++) {
MPI_Status status;
MPI_Recv(buffer.empty() ? NULL : &buffer.front(), nlocalatoms_max * 3, MPI_DOUBLE, iproc, 0, world, &status);
int ndoubles;
MPI_Get_count(&status, MPI_DOUBLE, &ndoubles);
int nReceived = ndoubles / 3;
vector<double>::const_iterator data = buffer.begin();
for(int i = 0; i < nReceived; i++, ++currentAtom) {
currentAtom->tag = currentAtom - inputAtoms.begin() + 1;
currentAtom->pos.X = *data++;
currentAtom->pos.Y = *data++;
currentAtom->pos.Z = *data++;
currentAtom->flags = 0;
currentAtom->cluster = NULL;
currentAtom->numNeighbors = 0;
currentAtom->setFlag(ATOM_IS_LOCAL_ATOM);
}
}
}
DISLOCATIONS_ASSERT(currentAtom == inputAtoms.end());
}
else {
// Allocate and fill send buffer.
vector<double> buffer(atom->nlocal * 3);
vector<double>::iterator data = buffer.begin();
for(int i = 0; i < atom->nlocal; i++) {
*data++ = atom->x[i][0];
*data++ = atom->x[i][1];
*data++ = atom->x[i][2];
}
// Send local atom coordinates to master proc.
MPI_Send(buffer.empty() ? NULL : &buffer.front(), buffer.size(), MPI_DOUBLE, 0, 0, world);
}
// Make sure all input atoms are wrapped at periodic boundary conditions.
wrapInputAtoms(NULL_VECTOR);
// Build nearest neighbor lists.
buildNearestNeighborLists();
}
/// This structure is used to communicate dislocation segments between processors.
struct DislocationSegmentComm {
LatticeVector burgersVector;
Vector3 burgersVectorWorld;
int numPoints;
};
/// After the master processor has extracted all dislocation segments, this function broadcasts the dislocations
/// back to all other processor.
void broadcastDislocations()
{
if(comm->nprocs == 1) return; // Nothing to do in serial mode.
// Broadcast number of segments.
int numSegments = segments.size();
MPI_Bcast(&numSegments, 1, MPI_INT, 0, world);
// Allocate send/receive buffer for dislocation segments.
vector<DislocationSegmentComm> segmentBuffer(numSegments);
// The total number of line points (sum of all segments).
int numLinePoints = 0;
if(processor == 0) {
// Fill segment send buffer.
vector<DislocationSegmentComm>::iterator sendItem = segmentBuffer.begin();
for(vector<DislocationSegment*>::const_iterator segment = segments.begin(); segment != segments.end(); ++segment, ++sendItem) {
sendItem->burgersVector = (*segment)->burgersVector;
sendItem->burgersVectorWorld = (*segment)->burgersVectorWorld;
sendItem->numPoints = (*segment)->line.size();
numLinePoints += (*segment)->line.size();
}
}
// Broadcast segments.
MPI_Bcast(segmentBuffer.empty() ? NULL : &segmentBuffer.front(), segmentBuffer.size() * sizeof(segmentBuffer[0]), MPI_CHAR, 0, world);
if(processor != 0) {
// Extract segments from receive buffer.
segments.reserve(segmentBuffer.size());
for(vector<DislocationSegmentComm>::const_iterator receiveItem = segmentBuffer.begin(); receiveItem != segmentBuffer.end(); ++receiveItem) {
DislocationSegment* newSegment = segmentPool.construct(receiveItem->burgersVector, receiveItem->burgersVectorWorld);
newSegment->index = segments.size() + 1;
segments.push_back(newSegment);
numLinePoints += receiveItem->numPoints;
}
}
// Allocate send/receive buffer for dislocation points.
vector<Point3> pointBuffer(numLinePoints);
if(processor == 0) {
// Fill point send buffer.
vector<Point3>::iterator sendItem = pointBuffer.begin();
for(vector<DislocationSegment*>::const_iterator segment = segments.begin(); segment != segments.end(); ++segment) {
for(deque<Point3>::const_iterator p = (*segment)->line.begin(); p != (*segment)->line.end(); ++p, ++sendItem)
*sendItem = *p;
}
DISLOCATIONS_ASSERT(sendItem == pointBuffer.end());
}
// Broadcast segments.
MPI_Bcast(pointBuffer.empty() ? NULL : &pointBuffer.front(), pointBuffer.size() * sizeof(pointBuffer[0]), MPI_CHAR, 0, world);
if(processor != 0) {
// Extract points from receive buffer.
vector<DislocationSegment*>::const_iterator segment = segments.begin();
vector<Point3>::const_iterator p = pointBuffer.begin();
for(vector<DislocationSegmentComm>::const_iterator receiveItem = segmentBuffer.begin(); receiveItem != segmentBuffer.end(); ++receiveItem, ++segment) {
(*segment)->line.assign(p, p + receiveItem->numPoints);
p += receiveItem->numPoints;
}
DISLOCATIONS_ASSERT(p == pointBuffer.end());
}
}
protected:
/// Pointer to the main LAMMPS object.
LAMMPS* lammps;
/// Enables serial processing, which provides exact handling of periodic boundary conditions.
bool exactMode;
/// The output stream to which log messages are sent.
stdio_osyncstream msgLoggerStream;
/// The output stream to which verbose log messages are sent.
stdio_osyncstream verboseLoggerStream;;
};
#endif // __DXA_DISLOCATION_EXTRACTOR_H
c4science · Help