Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/coupling/xiao.cc
5c77d804fe71features/coupler_compon…

Last Change
Raw File

xiao.cc
View Options

/**
* @file xiao.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Jan 15 17:00:43 2014
*
* @brief Bridging Domain/Arlequin method
*
* @section LICENSE
*
* Copyright INRIA and CEA
*
* The LibMultiScale is a C++ parallel framework for the multiscale
* coupling methods dedicated to material simulations. This framework
* provides an API which makes it possible to program coupled simulations
* and integration of already existing codes.
*
* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
* within ScAlApplix team and CEA/DPTA Ile de France.
* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
* within the LSMS/ENAC laboratory.
*
* This software is governed by the CeCILL-C license under French law and
* abiding by the rules of distribution of free software. You can use,
* modify and/ or redistribute the software under the terms of the CeCILL-C
* license as circulated by CEA, CNRS and INRIA at the following URL
* "http://www.cecill.info".
*
* As a counterpart to the access to the source code and rights to copy,
* modify and redistribute granted by the license, users are provided only
* with a limited warranty and the software's author, the holder of the
* economic rights, and the successive licensors have only limited
* liability.
*
* In this respect, the user's attention is drawn to the risks associated
* with loading, using, modifying and/or developing or reproducing the
* software by the user in light of its specific status of free software,
* that may mean that it is complicated to manipulate, and that also
* therefore means that it is reserved for developers and experienced
* professionals having in-depth computer knowledge. Users are therefore
* encouraged to load and test the software's suitability as regards their
* requirements in conditions enabling the security of their systems and/or
* data to be ensured and, more generally, to use and operate it in the
* same conditions as regards security.
*
* The fact that you are presently reading this means that you have had
* knowledge of the CeCILL-C license and that you accept its terms.
*
*/
#define TIMER
/* -------------------------------------------------------------------------- */
#include "xiao.hh"
#include "lib_bridging.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename DomainA, typename DomainC, UInt Dim>
void Xiao<DomainA, DomainC, Dim>::coupling(CouplingStage stage) {
if (!this->is_in_continuum && !this->is_in_atomic)
return;
if (stage != COUPLING_STEP3)
return;
STARTTIMER("BuildRHS");
DUMP("Clean RHS", DBG_INFO);
Parent::cleanRHS();
DUMP("Build RHS", DBG_INFO);
if (this->is_in_continuum)
Parent::buildContinuRHS();
STOPTIMER("BuildRHS");
DUMP("Synch with migration bridging", DBG_INFO);
STARTTIMER("updatesForMigrations");
Parent::bridging_zone.updateForMigration();
Parent::size_constraint = Parent::bridging_zone.getLocalPoints();
DUMP("Synch with migration boundary", DBG_INFO);
Parent::MDboundary_zone.updateForMigration();
STOPTIMER("updatesForMigrations");
DUMP("Correcting surface effect", DBG_INFO);
STARTTIMER("projectAtomicVelocitiesOnMe.hh");
// Parent::MDboundary_zone.projectAtomicDOFsOnMesh(_displacement);
Parent::MDboundary_zone.projectAtomicDOFsOnMesh(_velocity);
STOPTIMER("projectAtomicVelocitiesOnMe.hh");
STARTTIMER("BuildRHS");
DUMP("Build RHS", DBG_INFO);
if (this->is_in_atomic)
Parent::buildAtomsRHS();
STOPTIMER("BuildRHS");
// synchronize contributions by summing
DUMP("RHS redistrib/reduction", DBG_INFO);
STARTTIMER("SynchRHS");
Parent::bridging_zone.synchronizeVectorBySum("rhs", Parent::rhs, Dim);
STOPTIMER("SynchRHS");
DUMP("solving constraint", DBG_INFO);
STARTTIMER("Solving constraint");
Parent::solveConstraint();
STOPTIMER("Solving constraint");
DUMP("Applying constraint", DBG_INFO);
STARTTIMER("Correcting");
if (this->is_in_atomic)
Parent::applyAtomsCorrection();
if (this->is_in_continuum)
Parent::applyContinuCorrection();
STOPTIMER("Correcting");
DUMP(this->getID() << " : coupling stage all done", DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void Xiao<DomainA, DomainC, Dim>::init() {
// allocate the bridging zone (parallel version)
Parent::bridging_zone.setParam("GEOMETRY", this->bridging_geom);
Parent::bridging_zone.setParam("CHECK_COHERENCY", this->check_coherency);
Parent::bridging_zone.setParam("CENTROID", this->centroid_flag);
// set the pbc pairs
if (this->is_in_continuum)
Parent::bridging_zone.setPBCPairs(Parent::domC.getPBCpairs());
// initialize the bridging_zone object
Parent::bridging_zone.init();
// allocate the vectors necessary to continue
if (this->is_in_atomic || Parent::bridging_zone.getLocalPoints())
Parent::allocate(Parent::bridging_zone.getLocalPoints());
// build weights & dump vtu files with the weights attached & correct weights
if (this->is_in_continuum) {
Parent::computeContinuumWeights();
Parent::correctContinuumWeights();
}
if (this->is_in_atomic) {
Parent::computeAtomWeights();
Parent::correctAtomWeights();
}
/* build constraint matrix */
if (this->is_in_continuum)
Parent::buildContinuConstraintMatrix();
if (this->is_in_atomic)
Parent::buildAtomsConstraintMatrix();
/* synch constraint matrix */
Parent::bridging_zone.synchronizeVectorBySum("A", Parent::A);
/* now treat the boundary zone */
Parent::MDboundary_zone.setParam("GEOMETRY", this->boundary_geom);
Parent::MDboundary_zone.setParam("CHECK_COHERENCY", this->check_coherency);
Parent::MDboundary_zone.setParam("CENTROID", this->centroid_flag);
// initialize the bridging_zone object
Parent::MDboundary_zone.init();
// Parent::MDboundary_zone.projectAtomicDOFsOnMesh(_displacement);
Parent::MDboundary_zone.projectAtomicDOFsOnMesh(_velocity);
}
/* -------------------------------------------------------------------------- */
/* LMDESC XIAO
This class implements the BridgingDomain/Arlequin method.
Two zones are declared, one for the coupling (light blue atoms)
one for providing a stiff boundary condition to the atoms (purple atoms).
It defines a coupling area as depicted below: \\
\includegraphics{images/bridging-zone-BM} \\
In the bridging part a linear weight function is built the weight the
energies in the Arlequin method flavor.
The detailed algorithm is:\\
\begin{equation*}
\label{initial-BM}
\left\{
\begin{array}{l}
\displaystyle\hat{M}_I \ddot{\mathbf{u}}_I = - \mathbf{f}^R_I + \sum_{k=1}^L
\mathbf{\lambda}_k
\frac{\partial \mathbf{g}_k}{\partial \mathbf{u}_I} \\
\displaystyle \hat{m}_i \ddot{\mathbf{d}}_i = \mathbf{f}^R_i + \sum_{k=1}^L
\mathbf{\lambda}_k
\frac{\partial \mathbf{g}_k}{\partial \mathbf{d}_i}
\end{array}
\right.
\end{equation*}\\
where $\alpha(X_I)$ and $\alpha(X_i)$ are the weight associated with
nodes and atoms respectively, where
$\hat{M}_I=\alpha(X_I)M_I$ and $\hat{m}_i=(1-\alpha(X_i))m_i$
and where the constraint multipliers are computed by solving
$H \Lambda = \mathbf{g}^\star$ with:\\
\begin{align*}
&H_{ik} = \Delta t \left( \sum_J \varphi_J(\mathbf{X}_I) \hat{M_J}^{-1}
\frac{\partial \mathbf{g}_k}{\partial \mathbf{u}_J} -
\hat{m}_I^{-1} \frac{\partial \mathbf{g}_k}{\partial \mathbf{d}_i}
\right),\\
&\mathbf{g}_i^\star = \sum_J \varphi_J(\mathbf{X}_I) \dot{\mathbf{u}}_J^\star
- \dot{\mathbf{d}}_i^\star,
\end{align*}
where $u^\star, d^\star$ are the displacement obtained after one time step
by the integration scheme when the constraints are not applied. \\ \\
For debugging puposes several computes are registered to the central
system:\\
\begin{itemize}
\item weight-fe-\${couplerID} : it is a compute containing the weights
associated with nodes
inside of the bridging zone (the $\alpha(X_I)$).
\item weight-md-\${couplerID} : it is a compute containing the weights
associated with atoms
inside of the bridging zone (the $\alpha(X_i)$).
\item lambdas-fe-\${couplerID} : it is a compute containing the $\hat{M}_I$.
\item lambdas-md-\${couplerID} : it is a compute containing the $\hat{m}_i$.
\end{itemize}
*/
/* LMEXAMPLE COUPLING_CODE bdmID md fe XIAO GEOMETRY 3 BOUNDARY 4 GRID_DIVISIONX
* 10 GRID_DIVISIONY 10 QUALITY 1e-3 */
/* LMHERITANCE arlequin_template */
template <typename DomainA, typename DomainC, UInt Dim>
void Xiao<DomainA, DomainC, Dim>::declareParams() {
ArlequinTemplate<DomainA, DomainC, Dim>::declareParams();
/* LMKEYWORD CENTROID
Set the selection of the bridging zones based on centroid.
*/
this->parseTag("CENTROID", centroid_flag, false);
}
/* -------------------------------------------------------------------------- */
DECLARE_ATOMIC_CONTINUUM_TEMPLATE(Xiao)
__END_LIBMULTISCALE__
c4science · Help