trilinos_stokes_initialization_process.h

Go to the documentation of this file.

#ifndef KRATOS_TRILINOS_STOKES_INITIALIZATION_PROCESS_H
#define KRATOS_TRILINOS_STOKES_INITIALIZATION_PROCESS_H
 
// System includes
 
// External includes
#include "Epetra_MpiComm.h"
 
// Project includes
 
// Application includes
#include "containers/model.h"
#include "solving_strategies/schemes/residualbased_incrementalupdate_static_scheme.h"
#include "custom_strategies/builder_and_solvers/trilinos_block_builder_and_solver_periodic.h"
#include "mpi/includes/mpi_communicator.h"
#include "mpi/utilities/parallel_fill_communicator.h"
 
#include "../FluidDynamicsApplication/custom_processes/stokes_initialization_process.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
template<class TSparseSpace,
         class TDenseSpace,
         class TLinearSolver
         >
class  TrilinosStokesInitializationProcess : public StokesInitializationProcess<TSparseSpace,TDenseSpace,TLinearSolver>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( TrilinosStokesInitializationProcess);
 
    typedef StokesInitializationProcess<TSparseSpace,TDenseSpace,TLinearSolver> BaseType;
 
 
    TrilinosStokesInitializationProcess(Epetra_MpiComm& rComm,
                                        ModelPart& rModelPart,
                                        typename TLinearSolver::Pointer pLinearSolver,
                                        unsigned int DomainSize,
                                        const Variable<int>& PeriodicPairIndicesVar):
        BaseType(rModelPart,pLinearSolver,DomainSize,this),
        mrComm(rComm),
        mrPeriodicVar(PeriodicPairIndicesVar)
    {
        KRATOS_TRY;
 
        ModelPart& rReferenceModelPart = BaseType::mrReferenceModelPart;
        typename TLinearSolver::Pointer& pBaseLinearSolver = BaseType::mpLinearSolver;
        unsigned int BaseDomainSize = BaseType::mDomainSize;
        typename ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver>::Pointer& pSolutionStrategy = BaseType::mpSolutionStrategy;
 
        // Initialize new model part (same nodes, new elements, no conditions)
        if(rModelPart.GetModel().HasModelPart("StokesModelPart"))
            rModelPart.GetModel().DeleteModelPart("StokesModelPart");
        ModelPart& rStokesModelPart = rModelPart.GetModel().CreateModelPart("StokesModelPart");
        rStokesModelPart.GetNodalSolutionStepVariablesList() = rReferenceModelPart.GetNodalSolutionStepVariablesList();
        rStokesModelPart.SetBufferSize(1);
        rStokesModelPart.SetNodes( rReferenceModelPart.pNodes() );
        rStokesModelPart.SetProcessInfo( rReferenceModelPart.pGetProcessInfo() );
        rStokesModelPart.SetProperties( rReferenceModelPart.pProperties() );
 
        // Create a communicator for the new model part and copy the partition information about nodes.
        Communicator& rReferenceComm = rReferenceModelPart.GetCommunicator();
        typename Communicator::Pointer pStokesMPIComm = Kratos::make_shared<MPICommunicator>( &(rReferenceModelPart.GetNodalSolutionStepVariablesList()) );
        pStokesMPIComm->SetNumberOfColors( rReferenceComm.GetNumberOfColors() ) ;
        pStokesMPIComm->NeighbourIndices() = rReferenceComm.NeighbourIndices();
        pStokesMPIComm->LocalMesh().SetNodes( rReferenceComm.LocalMesh().pNodes() );
        pStokesMPIComm->InterfaceMesh().SetNodes( rReferenceComm.InterfaceMesh().pNodes() );
        pStokesMPIComm->GhostMesh().SetNodes( rReferenceComm.GhostMesh().pNodes() );
        for (unsigned int i = 0; i < rReferenceComm.GetNumberOfColors(); i++)
        {
            pStokesMPIComm->pInterfaceMesh(i)->SetNodes( rReferenceComm.pInterfaceMesh(i)->pNodes() );
            pStokesMPIComm->pLocalMesh(i)->SetNodes( rReferenceComm.pLocalMesh(i)->pNodes() );
            pStokesMPIComm->pGhostMesh(i)->SetNodes( rReferenceComm.pGhostMesh(i)->pNodes() );
        }
        rStokesModelPart.SetCommunicator( pStokesMPIComm );
 
        // Retrieve Stokes element model
        std::string ElementName;
        if (BaseDomainSize == 2)
            ElementName = std::string("StationaryStokes2D");
        else
            ElementName = std::string("StationaryStokes3D");
 
        const Element& rReferenceElement = KratosComponents<Element>::Get(ElementName);
 
        // Generate Stokes elements
        for (ModelPart::ElementsContainerType::iterator itElem = rReferenceModelPart.ElementsBegin(); itElem != rReferenceModelPart.ElementsEnd(); itElem++)
        {
            Element::Pointer pElem = rReferenceElement.Create(itElem->Id(), itElem->GetGeometry(), itElem->pGetProperties() );
            rStokesModelPart.Elements().push_back(pElem);
            rStokesModelPart.GetCommunicator().LocalMesh().Elements().push_back(pElem);
        }
 
        // Solution scheme: Linear static scheme
        auto pScheme = Kratos::make_shared< ResidualBasedIncrementalUpdateStaticScheme< TSparseSpace, TDenseSpace > >();
 
        // Builder and solver
        int guess_row_size;
        if(BaseDomainSize == 2)
            guess_row_size = 15;
        else
            guess_row_size = 40;
 
        auto pBuildAndSolver = Kratos::make_shared<TrilinosBlockBuilderAndSolverPeriodic<TSparseSpace,TDenseSpace,TLinearSolver> >(
            mrComm,
            guess_row_size,
            pBaseLinearSolver,
            mrPeriodicVar
        );
 
        // Strategy
        bool ReactionFlag = false;
        bool ReformDofSetFlag = false;
        bool CalculateNormDxFlag = false;
        bool MoveMeshFlag = false;
 
        pSolutionStrategy = Kratos::make_shared< ResidualBasedLinearStrategy<TSparseSpace, TDenseSpace, TLinearSolver> >(
            rStokesModelPart,
            pScheme,
            pBuildAndSolver,
            ReactionFlag,
            ReformDofSetFlag,
            CalculateNormDxFlag,
            MoveMeshFlag
        );
 
        pSolutionStrategy->SetEchoLevel(0);
        pSolutionStrategy->Check();
 
        BaseType::mIsCleared = false;
 
        KRATOS_CATCH("");
    }
 
 
    virtual ~TrilinosStokesInitializationProcess()
    {
        BaseType::mpSolutionStrategy->Clear();
    }
 
 
 
 
 
 
 
 
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << " TrilinosStokesInitializationProcess";
        return buffer.str();
    }
 
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << " TrilinosStokesInitializationProcess";
    }
 
 
    void PrintData(std::ostream& rOStream) const override
    {
    }
 
 
 
 
 
protected:
 
 
 
    Epetra_MpiComm& mrComm;
 
    const Variable<int>& mrPeriodicVar;
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
    TrilinosStokesInitializationProcess & operator=( TrilinosStokesInitializationProcess const& rOther)
    {
        return *this;
    }
 
 
    TrilinosStokesInitializationProcess( TrilinosStokesInitializationProcess const& rOther)
    { }
 
 
 
}; // Class  TrilinosStokesInitializationProcess
 
 
 
 
 
 
 
template<class TSparseSpace,
         class TDenseSpace,
         class TLinearSolver
         >
inline std::istream & operator >>(std::istream& rIStream,
                                  TrilinosStokesInitializationProcess<TSparseSpace, TDenseSpace, TLinearSolver>& rThis)
{
    return rIStream;
}
 
 
template<class TSparseSpace,
         class TDenseSpace,
         class TLinearSolver
         >
inline std::ostream & operator <<(std::ostream& rOStream,
                                  const  TrilinosStokesInitializationProcess<TSparseSpace, TDenseSpace, TLinearSolver>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
} // namespace Kratos.
 
#endif // KRATOS_TRILINOS_STOKES_INITIALIZATION_PROCESS_H