Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/coupling/kobayashi.cc
b0a1ddf8eccclammps_2020

Last Change
Raw File

kobayashi.cc
View Options

/**
* @file kobayashi.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Fri Jul 11 15:47:44 2014
*
* @brief Kobayashi's bridging method
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "kobayashi.hh"
#include "arlequin_template.hh"
#include "lib_bridging.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
Kobayashi::Kobayashi(const std::string &name)
: LMObject(name), ArlequinTemplate(name) {
// naming the bridging zones
std::stringstream fname;
fname << name << "-bridging";
}
/* -------------------------------------------------------------------------- */
Kobayashi::~Kobayashi() {}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Kobayashi::init(DomainA &domA, DomainC &domC) {
const UInt Dim = spatial_dimension;
if (lm_world_size > 1)
LM_FATAL("this coupler work only in sequential");
std::stringstream fname;
// initialize the bridging_zone object
bridging_zone.init(domA, domC);
// allocate the vectors necessary to continue
UInt size = bridging_zone.getNumberPoints();
this->allocate(size);
// build weights
this->computeContinuumWeights(domC.getContainer());
this->computeAtomWeights(domA.getContainer());
UInt _size = bridging_zone.getNumberElems();
DUMP("Allocating vectors", DBG_INFO);
rhs.resize(size, Dim);
CG_vel.resize(_size, Dim);
CG_disp.resize(_size, Dim);
DUMP("Computing initial least square constraint", DBG_INFO);
buildCGMatrix();
// treat now MD boundary zone
MDboundary_zone.init(domA, domC);
DUMP("Init ... done", DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void Kobayashi::coupling(CouplingStage stage, DomainA &domA, DomainC &domC) {
if (stage != COUPLING_STEP4)
return;
MDboundary_zone.projectAtomicFieldOnMesh(_velocity);
MDboundary_zone.projectAtomicFieldOnMesh(_displacement);
buildLeastSquareRHS(domA.getContainer());
correctContinuum(domC.getContainer());
correctAtoms(domA.getContainer());
}
/* -------------------------------------------------------------------------- */
void Kobayashi::buildCGMatrix() {
UInt nb_nodes = bridging_zone.getNumberNodes();
UInt nb_atoms [[gnu::unused]] = bridging_zone.getNumberPoints();
Eigen::MatrixXf ls_mat(nb_nodes, nb_nodes);
bridging_zone.buildLeastSquareMatrix(ls_mat);
LM_TOIMPLEMENT;
// ls_mat.factoLU();
// ls_mat.inverse();
// math::Matrix shapes(nb_atoms, nb_nodes);
// bridging_zone.extractShapeMatrix(shapes);
// shapes.transpose();
// A = new math::Matrix(nb_nodes, nb_atoms);
// A->matProduct(ls_mat, shapes);
}
/* -------------------------------------------------------------------------- */
template <typename ContA>
void Kobayashi::buildLeastSquareRHS(ContA &pointList) {
rhs.setZero();
UInt at_index = 0;
for (auto &&at : pointList) {
rhs(at_index) = at.velocity()[0];
++at_index;
}
CG_vel = A * rhs.matrix();
at_index = 0;
for (auto &&at : pointList) {
rhs(at_index) = at.displacement()[0];
++at_index;
}
CG_disp = A * rhs.matrix();
}
/* -------------------------------------------------------------------------- */
template <typename ContB> void Kobayashi::correctContinuum(ContB &meshList) {
auto &weightFE = this->weightFE.getArray();
UInt j = 0;
for (auto &&n : meshList) {
if (weightFE[j] <= .5) {
Real speed_projection = this->quality;
if (weightFE[j] <= 1e-2)
speed_projection = 1;
Real &vel = n.velocity()[0];
Real &disp = n.displacement()[0];
Real corrected = vel + speed_projection * (CG_vel(j) - vel);
n.velocity()[0] = corrected;
corrected = disp + speed_projection * (CG_disp(j) - disp);
n.displacement()[0] = corrected;
}
++j;
}
}
/* -------------------------------------------------------------------------- */
template <typename ContA> void Kobayashi::correctAtoms(ContA &pointList) {
UInt natom = bridging_zone.getNumberLocalMatchedPoints();
constexpr UInt Dim = ContA::Dim;
ContainerArray<Real> atomic_vel;
ContainerArray<Real> atomic_disp;
atomic_vel.assign(natom * Dim, 0);
atomic_disp.assign(natom * Dim, 0);
bridging_zone.projectAtomicFieldOnMesh(_displacement, atomic_disp);
bridging_zone.projectAtomicFieldOnMesh(_velocity, atomic_vel);
auto &weightMD = this->weightMD.getArray();
UInt at_index = 0;
for (auto &&at : pointList) {
if (weightMD[at_index] <= .5) {
Real speed_projection = this->quality;
if (weightMD[at_index] <= 1e-2)
speed_projection = 1;
Real &vel = at.velocity()[0];
Real disp = at.displacement()[0];
vel = vel + speed_projection * (atomic_vel[at_index] - vel);
Vector<Dim> corrected;
corrected[0] = disp + this->quality * (atomic_disp[at_index] - disp);
at.displacement() = corrected;
}
++at_index;
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC KOBAYASHI
This keyword is used to set a coupling method
which will implement the method of Kobayashi and co-workers \\ \\
\textit{A simple dynamical scale-coupling method for concurrent simulation of
hybridized atomistic/coarse-grained-particle system,} \\
\textbf{Ryo Kobayashi, Takahide Nakamura, Shuji Ogata,} \\
\textit{International Journal for Numerical Methods in Engineering Volume 83,
Issue 2, pages 249-268, 9 July 2010.} \\ \\
This coupling component works only in 1D and in sequential.
*/
/* LMHERITANCE arlequin_template */
/* LMEXAMPLE
COUPLING_CODE md fe KOBAYASHI GEOMETRY 3 BOUNDARY 4 GRID_DIVISIONX 10 QUALITY
1e-3 VERBOSE
*/
void Kobayashi::declareParams() { ArlequinTemplate::declareParams(); }
/* -------------------------------------------------------------------------- */
DECLARE_COUPLER_MAKE_CALL(Kobayashi, domA, domC)
__END_LIBMULTISCALE__
c4science · Help