Homec4science
Diffusion lammps (lammps-icms) Browse

lammps/lib/atc/FieldManager.cpp
e7e5ba4e4688lammps-icms

Last Change
Raw File

FieldManager.cpp
View Options

#include "FieldManager.h"
#include "ATC_Method.h"
#include "LammpsInterface.h"
#include "PerAtomQuantity.h"
#include "TransferOperator.h"
using std::string;
namespace ATC {
typedef PerAtomQuantity<double> PAQ;
//-----------------------------------------------------------------------------
//*
//-----------------------------------------------------------------------------
FieldManager::FieldManager(ATC_Method * atc):
atc_(atc),
interscaleManager_(atc->interscale_manager())
{};
//-----------------------------------------------------------------------------
//* restricted_atom_quantity
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::restricted_atom_quantity(FieldName field, string name, PAQ * atomicQuantity)
{
if (name == "default") { name = field_to_restriction_name(field); }
DENS_MAN * quantity = interscaleManager_.dense_matrix(name);
if (!quantity){
if (field == CHARGE_DENSITY) {
atomicQuantity = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_CHARGE);
}
else if (field == MASS_DENSITY) {
atomicQuantity = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
}
else {
if (!atomicQuantity) {
throw ATC_Error("FieldManager::restricted_atom_quantity - need to supply PAQ if restricted quantity does not already exist");
}
}
quantity = new AtfShapeFunctionRestriction(atc_,atomicQuantity,atc_->accumulant());
interscaleManager_.add_dense_matrix(quantity,name);
}
return quantity;
}
//-----------------------------------------------------------------------------
//* restricted_atom_quantity
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::projected_atom_quantity(FieldName field,string name, PAQ * atomic, DIAG_MAN * normalization)
{
if (atc_->use_md_mass_normalization()) {
if (name == "default") { name = field_to_intrinsic_name(field); }
DENS_MAN * quantity = interscaleManager_.dense_matrix(name);
if (!quantity) {
DENS_MAN * restricted = restricted_atom_quantity(field,field_to_restriction_name(field),atomic);
quantity = new AtfShapeFunctionMdProjection(atc_,restricted,use_mass_matrix(field));
interscaleManager_.add_dense_matrix(quantity,name);
}
return quantity;
}
else {
if (name == "default") { name = field_to_string(field); }
DENS_MAN * quantity = interscaleManager_.dense_matrix(name);
if (quantity) return quantity;
if (atc_->kernel_on_the_fly()) {
if (atc_->kernel_based()) {
quantity = new OnTheFlyKernelAccumulationNormalized(atc_, atomic,
atc_->kernel_function(),
atc_->atom_coarsegraining_positions(),
normalization);
} else {
quantity = new OnTheFlyMeshAccumulationNormalized(atc_, atomic,
atc_->atom_coarsegraining_positions(),
normalization);
}
} else {
quantity = new AtfProjection(atc_, atomic,
atc_->accumulant(),
normalization);
}
interscaleManager_.add_dense_matrix(quantity,name);
return quantity;
}
}
//-----------------------------------------------------------------------------
//* referenced_projected_atom_quantity
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::referenced_projected_atom_quantity(FieldName field,string name, PAQ * atomic, DENS_MAN * reference, DIAG_MAN * normalization)
{
if (name == "default") { name = field_to_string(field); }
DENS_MAN * quantity = interscaleManager_.dense_matrix(name);
if (quantity) return quantity;
if (atc_->use_md_mass_normalization()) {
DENS_MAN * restricted = restricted_atom_quantity(field,field_to_restriction_name(field),atomic);
quantity = new AtfShapeFunctionMdProjectionReferenced(atc_,restricted,reference,use_mass_matrix(field));
}
else if (atc_->kernel_on_the_fly()) {
if (atc_->kernel_based()) {
quantity = new OnTheFlyKernelAccumulationNormalizedReferenced(atc_,
atomic,
atc_->kernel_function(),
atc_->atom_coarsegraining_positions(),
normalization,
reference);
} else {
quantity = new OnTheFlyMeshAccumulationNormalizedReferenced(atc_,
atomic,
atc_->atom_coarsegraining_positions(),
normalization,
reference);
}
} else {
quantity = new AtfProjectionReferenced(atc_, atomic,
atc_->accumulant(),
reference,
normalization);
}
interscaleManager_.add_dense_matrix(quantity,name);
return quantity;
}
//-----------------------------------------------------------------------------
//* scaled_projected_atom_quantity
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::scaled_projected_atom_quantity(FieldName field,string name, PAQ * atomic, double scale, DIAG_MAN * normalization)
{
if (name == "default") { name = field_to_string(field); }
DENS_MAN * quantity = interscaleManager_.dense_matrix(name);
if (quantity) return quantity;
if (atc_->use_md_mass_normalization()) {
DENS_MAN * restricted = restricted_atom_quantity(field,field_to_restriction_name(field),atomic);
quantity = new AtfShapeFunctionMdProjectionScaled(atc_,restricted,scale,use_mass_matrix(field));
}
else if (atc_->kernel_on_the_fly()) {
if (atc_->kernel_based()) {
quantity = new OnTheFlyKernelAccumulationNormalizedScaled(atc_, atomic,
atc_->kernel_function(),
atc_->atom_coarsegraining_positions(),
normalization,
scale);
} else {
quantity = new OnTheFlyMeshAccumulationNormalizedScaled(atc_, atomic,
atc_->atom_coarsegraining_positions(),
normalization,
scale);
}
} else {
quantity = new AtfProjectionScaled(atc_, atomic,
atc_->accumulant(),
scale,
normalization);
}
interscaleManager_.add_dense_matrix(quantity,name);
return quantity;
}
//-----------------------------------------------------------------------------
//* CHARGE_DENSITY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::charge_density(string name)
{
FundamentalAtomQuantity * atomic = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_CHARGE);
return projected_atom_quantity(CHARGE_DENSITY,name,atomic,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* MASS_DENSITY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::mass_density(string name)
{
FundamentalAtomQuantity * atomic = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
return projected_atom_quantity(MASS_DENSITY,name,atomic,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* SPECIES_CONCENTRATION
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::species_concentration(string name)
{
PerAtomQuantity<double> * atomTypeVector = interscaleManager_.per_atom_quantity("atom_species_vector");
if (!atomTypeVector) {
atomTypeVector = new AtomTypeVector(atc_,atc_->type_list(),atc_->group_list());
interscaleManager_.add_per_atom_quantity(atomTypeVector,"atom_species_vector");
}
return projected_atom_quantity(SPECIES_CONCENTRATION,name,atomTypeVector,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* NUMBER_DENSITY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::number_density(string name)
{
PAQ * atomic = interscaleManager_.per_atom_quantity("atomNumber");
if (!atomic) {
atomic = new AtomNumber(atc_);
interscaleManager_.add_per_atom_quantity(atomic, "atomNumber");
}
return projected_atom_quantity(NUMBER_DENSITY,name,atomic,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* MOMENTUM
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::momentum(string name)
{
PAQ * atomic = interscaleManager_.per_atom_quantity("atomMomentum");
if (!atomic) {
atomic = new AtomicMomentum(atc_);
interscaleManager_.add_per_atom_quantity(atomic, "atomMomentum");
}
return projected_atom_quantity(MOMENTUM,name,atomic,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* VELOCITY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::velocity(string name)
{
if (name == "default") { name = field_to_string(VELOCITY); }
DENS_MAN * v = interscaleManager_.dense_matrix(name);
if (v) return v;
if (atc_->use_md_mass_normalization()) {
PAQ * atomic = interscaleManager_.per_atom_quantity("atomMomentum");
if (!atomic) {
atomic = new AtomicMomentum(atc_);
interscaleManager_.add_per_atom_quantity(atomic, "atomMomentum");
}
DENS_MAN * restricted = restricted_atom_quantity(VELOCITY,field_to_restriction_name(VELOCITY),atomic);
v = new AtfShapeFunctionMdProjection(atc_,restricted,VELOCITY);
}
else {
DENS_MAN* p = interscaleManager_.dense_matrix(field_to_string(MOMENTUM));
if (!p) p = nodal_atomic_field(MOMENTUM);
DENS_MAN* m = interscaleManager_.dense_matrix(field_to_string(MASS_DENSITY));
if (!m) m = nodal_atomic_field(MASS_DENSITY);
v = new DenseMatrixQuotient(p,m);
}
interscaleManager_.add_dense_matrix(v,field_to_string(VELOCITY));
return v;
}
//-----------------------------------------------------------------------------
//* PROJECTED_VELOCITY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::projected_velocity(string name)
{
FundamentalAtomQuantity * atomic = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
return projected_atom_quantity(PROJECTED_VELOCITY,name,atomic,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* DISPLACEMENT
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::displacement(string name)
{
if (name == "default") { name = field_to_string(DISPLACEMENT); }
DENS_MAN * u = interscaleManager_.dense_matrix(name);
if (u) return u;
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicMassWeightedDisplacement");
if (!atomic) {
FundamentalAtomQuantity * atomMasses = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
FundamentalAtomQuantity * atomPositions = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION);
atomic = new AtomicMassWeightedDisplacement(atc_,atomPositions, atomMasses, atc_->atom_reference_positions(), INTERNAL);
interscaleManager_.add_per_atom_quantity(atomic,"AtomicMassWeightedDisplacement");
}
if (atc_->use_md_mass_normalization()) {
DENS_MAN * restricted = restricted_atom_quantity(DISPLACEMENT,field_to_restriction_name(DISPLACEMENT),atomic);
u = new AtfShapeFunctionMdProjection(atc_,restricted,VELOCITY);
}
else {
DENS_MAN * q = NULL;
if (atc_->kernel_on_the_fly()) {
if (atc_->kernel_based()) {
q = new OnTheFlyKernelAccumulationNormalized(atc_, atomic,
atc_->kernel_function(),
atc_->atom_coarsegraining_positions(),
atc_->accumulant_inverse_volumes());
} else {
q = new OnTheFlyMeshAccumulationNormalized(atc_, atomic,
atc_->atom_coarsegraining_positions(),
atc_->accumulant_inverse_volumes());
}
} else {
q = new AtfProjection(atc_, atomic,
atc_->accumulant(),
atc_->accumulant_inverse_volumes());
}
interscaleManager_.add_dense_matrix(q,"CoarseGrainedAMWD");
DENS_MAN* m = interscaleManager_.dense_matrix(field_to_string(MASS_DENSITY));
u = new DenseMatrixQuotient(q,m);
}
interscaleManager_.add_dense_matrix(u,name);
return u;
}
//-----------------------------------------------------------------------------
//* REFERENCE_POTENTIAL_ENERGY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::reference_potential_energy(string name)
{
DENS_MAN * rpe = interscaleManager_.dense_matrix(field_to_string(REFERENCE_POTENTIAL_ENERGY));
if (! rpe ) {
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicReferencePotential");
if (!atomic) {
atomic = new AtcAtomQuantity<double>(atc_);
interscaleManager_.add_per_atom_quantity(atomic, "AtomicReferencePotential");
atomic->set_memory_type(PERSISTENT);
}
SPAR_MAN * referenceAccumulant = interscaleManager_.sparse_matrix("ReferenceAccumulant");
if (!referenceAccumulant) {
referenceAccumulant = new SPAR_MAN();
referenceAccumulant->set_memory_type(PERSISTENT);
interscaleManager_.add_sparse_matrix(referenceAccumulant,
"ReferenceAccumulant");
}
DIAG_MAN * referenceAccumulantInverseVolumes = interscaleManager_.diagonal_matrix("ReferenceAccumulantInverseVolumes");
if (!referenceAccumulantInverseVolumes) {
referenceAccumulantInverseVolumes = new DIAG_MAN();
referenceAccumulantInverseVolumes->set_memory_type(PERSISTENT);
interscaleManager_.add_diagonal_matrix(referenceAccumulantInverseVolumes,
"ReferenceAccumulantInverseVolumes");
}
rpe = new AtfProjection(atc_, atomic,
referenceAccumulant,
referenceAccumulantInverseVolumes);
interscaleManager_.add_dense_matrix(rpe,field_to_string(REFERENCE_POTENTIAL_ENERGY));
rpe->set_memory_type(PERSISTENT);
}
return rpe;
}
//-----------------------------------------------------------------------------
//* POTENTIAL_ENERGY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::potential_energy(string name)
{
PerAtomQuantity<double> * atomic = interscaleManager_.per_atom_quantity("AtomicPotentialEnergy");
if (!atomic) {
atomic = new ComputedAtomQuantity(atc_,(atc_->lammps_interface())->compute_pe_name(), atc_->pe_scale());
interscaleManager_.add_per_atom_quantity(atomic,"AtomicPotentialEnergy");
}
DENS_MAN * reference = interscaleManager_.dense_matrix(field_to_string(REFERENCE_POTENTIAL_ENERGY));
if (reference) {
return referenced_projected_atom_quantity(POTENTIAL_ENERGY,name,atomic,reference,atc_->accumulant_inverse_volumes());
}
else {
return projected_atom_quantity(POTENTIAL_ENERGY,name,atomic,atc_->accumulant_inverse_volumes());
}
}
//-----------------------------------------------------------------------------
//* TEMPERATURE
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::temperature(string name)
{
double Tcoef = 1./((atc_->nsd())*(atc_->lammps_interface())->kBoltzmann());
PAQ * atomic = per_atom_quantity("AtomicTwiceFluctuatingKineticEnergy");
return scaled_projected_atom_quantity(TEMPERATURE,name,atomic,Tcoef,atc_->accumulant_weights());
}
//-----------------------------------------------------------------------------
//* KINETIC_TEMPERATURE
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::kinetic_temperature(string name)
{
double Tcoef = 1./((atc_->nsd())*(atc_->lammps_interface())->kBoltzmann());
PAQ * atomic = per_atom_quantity("AtomicTwiceKineticEnergy");
return scaled_projected_atom_quantity(KINETIC_TEMPERATURE,name,atomic,Tcoef,atc_->accumulant_weights());
}
//-----------------------------------------------------------------------------
//* THERMAL_ENERGY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::thermal_energy(string name)
{
double Ecoef = 0.5*atc_->ke_scale();
PAQ * atomic = per_atom_quantity("AtomicTwiceFluctuatingKineticEnergy");
return scaled_projected_atom_quantity(THERMAL_ENERGY,name,atomic,Ecoef,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* KINETIC_ENERGY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::kinetic_energy(string name)
{
double Ecoef = 0.5*atc_->ke_scale();
PAQ * atomic = per_atom_quantity("AtomicTwiceKineticEnergy");
return scaled_projected_atom_quantity(KINETIC_ENERGY,name,atomic,Ecoef,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* CHARGE_FLUX
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::charge_flux(string name)
{
PAQ * atomic = per_atom_quantity("AtomicChargeVelocity");
return projected_atom_quantity(CHARGE_FLUX,name,atomic,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* SPECIES_FLUX
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::species_flux(string name)
{
PAQ * atomic = per_atom_quantity("AtomicSpeciesVelocity");
return projected_atom_quantity(SPECIES_FLUX,name,atomic,atc_->accumulant_inverse_volumes());
}
//-----------------------------------------------------------------------------
//* INTERNAL_ENERGY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::internal_energy(string name)
{
if (name == "default") { name = field_to_string(INTERNAL_ENERGY); }
DENS_MAN * te = thermal_energy(field_to_string(THERMAL_ENERGY));
DENS_MAN * pe = potential_energy(field_to_string(POTENTIAL_ENERGY));
DenseMatrixSum * ie = new DenseMatrixSum(te,pe);
interscaleManager_.add_dense_matrix(ie,name);
return ie;
}
//-----------------------------------------------------------------------------
//* ENERGY
//-----------------------------------------------------------------------------
DENS_MAN * FieldManager::energy(string name)
{
if (name == "default") { name = field_to_string(ENERGY); }
DENS_MAN * ke = kinetic_energy(field_to_string(KINETIC_ENERGY));
DENS_MAN * pe = potential_energy(field_to_string(POTENTIAL_ENERGY));
DenseMatrixSum * e = new DenseMatrixSum(ke,pe);
interscaleManager_.add_dense_matrix(e,name);
return e;
}
//=============================================================================
//* PER ATOM QUANTITIES
//=============================================================================
//-----------------------------------------------------------------------------
//* 2 KE '
//-----------------------------------------------------------------------------
PAQ * FieldManager::atomic_twice_fluctuating_kinetic_energy()
{
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicTwiceFluctuatingKineticEnergy");
if (!atomic) {
FundamentalAtomQuantity * atomMass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
FundamentalAtomQuantity * atomVelocity = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
PAQ * vbar = per_atom_quantity(field_to_prolongation_name(VELOCITY));
atomic = new TwiceFluctuatingKineticEnergy(atc_,atomVelocity,atomMass,vbar);
interscaleManager_.add_per_atom_quantity(atomic, "AtomicTwiceFluctuatingKineticEnergy");
}
return atomic;
}
//-----------------------------------------------------------------------------
//* 2 KE
//-----------------------------------------------------------------------------
PAQ * FieldManager::atomic_twice_kinetic_energy()
{
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicTwiceKineticEnergy");
if (!atomic) {
FundamentalAtomQuantity * atomMass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
FundamentalAtomQuantity * atomVelocity = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
atomic = new TwiceKineticEnergy(atc_,atomVelocity,atomMass);
interscaleManager_.add_per_atom_quantity(atomic, "AtomicTwiceKineticEnergy");
}
return atomic;
}
//-----------------------------------------------------------------------------
//* v'
//-----------------------------------------------------------------------------
PAQ * FieldManager::atomic_fluctuating_velocity()
{
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicFluctuatingVelocity");
if (!atomic) {
FundamentalAtomQuantity * atomVelocity = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
PAQ * atomMeanVelocity = per_atom_quantity(field_to_prolongation_name(VELOCITY));
atomic = new FluctuatingVelocity(atc_,atomVelocity,atomMeanVelocity);
interscaleManager_.add_per_atom_quantity(atomic, "AtomicFluctuatingVelocity");
}
return atomic;
}
//-----------------------------------------------------------------------------
//* q v'
//-----------------------------------------------------------------------------
PAQ * FieldManager::atomic_charge_velocity()
{
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicChargeVelocity");
if (!atomic) {
PAQ * atomVelocity = atomic_fluctuating_velocity();
FundamentalAtomQuantity * atomCharge = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_CHARGE);
atomic = new ChargeVelocity(atc_,atomVelocity,atomCharge);
interscaleManager_.add_per_atom_quantity(atomic, "AtomicChargeVelocity");
}
return atomic;
}
//-----------------------------------------------------------------------------
//* m^a v'
//-----------------------------------------------------------------------------
PAQ * FieldManager::atomic_species_velocity()
{
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicSpeciesVelocity");
if (!atomic) {
PAQ * atomVelocity = atomic_fluctuating_velocity();
PAQ * atomSpecies = atomic_species_vector();
atomic = new SpeciesVelocity(atc_,atomVelocity,atomSpecies);
interscaleManager_.add_per_atom_quantity(atomic, "AtomicSpeciesVelocity");
}
return atomic;
}
//-----------------------------------------------------------------------------
//* [0 1 0 0 ] for type 2 atom
//-----------------------------------------------------------------------------
PAQ * FieldManager::atomic_species_vector()
{
PAQ * atomic = interscaleManager_.per_atom_quantity("AtomicSpeciesVector");
if (!atomic) {
atomic = new AtomTypeVector(atc_,atc_->type_list(),atc_->group_list());
interscaleManager_.add_per_atom_quantity(atomic,"AtomicSpeciesVector");
}
return atomic;
}
//-----------------------------------------------------------------------------
//* Prolonged coarse scale field
//-----------------------------------------------------------------------------
PAQ * FieldManager::prolonged_field(FieldName field)
{
PAQ * quantity = interscaleManager_.per_atom_quantity(field_to_prolongation_name(field));
if (!quantity) {
DENS_MAN * coarseQuantity = interscaleManager_.dense_matrix(field_to_string(field));
if (!coarseQuantity) coarseQuantity = nodal_atomic_field(field);
if (!coarseQuantity) throw ATC_Error("can not prolong quantity: " + field_to_string(field) + " no field registered");
if (atc_->kernel_on_the_fly()) {
quantity = new OnTheFlyShapeFunctionProlongation(atc_,
coarseQuantity,atc_->atom_coarsegraining_positions());
} else {
quantity = new FtaShapeFunctionProlongation(atc_,
coarseQuantity,atc_->interpolant());
}
interscaleManager_.add_per_atom_quantity(quantity,
field_to_prolongation_name(field));
}
return quantity;
}
}
c4science · Help