trilinos_monotonicity_preserving_solver.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Daniel Diez
//
 
#if !defined(KRATOS_TRILINOS_MONOTONICITY_PRESERVING_SOLVER_H_INCLUDED )
#define  KRATOS_TRILINOS_MONOTONICITY_PRESERVING_SOLVER_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "factories/linear_solver_factory.h"
#include "linear_solvers/linear_solver.h"
#include "utilities/parallel_utilities.h"
 
namespace Kratos
{
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TReordererType = Reorderer<TSparseSpaceType, TDenseSpaceType> >
class TrilinosMonotonicityPreservingSolver
    : public LinearSolver<TSparseSpaceType, TDenseSpaceType,  TReordererType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(TrilinosMonotonicityPreservingSolver);
 
    typedef LinearSolver<TSparseSpaceType, TDenseSpaceType, TReordererType> BaseType;
 
    typedef typename TSparseSpaceType::MatrixType SparseMatrixType;
 
    typedef typename TSparseSpaceType::VectorType VectorType;
 
    typedef typename TDenseSpaceType::MatrixType DenseMatrixType;
 
    typedef LinearSolverFactory<TSparseSpaceType,TDenseSpaceType> LinearSolverFactoryType;
 
 
    TrilinosMonotonicityPreservingSolver()
    {
    }
 
    TrilinosMonotonicityPreservingSolver(
        typename BaseType::Pointer pLinearSolver
        ) : BaseType (),
            mpLinearSolver(pLinearSolver)
    {
    }
 
    TrilinosMonotonicityPreservingSolver(Parameters ThisParameters)
        : BaseType ()
    {
        KRATOS_TRY
 
        Parameters default_parameters( R"(
        {
            "solver_type"                    : "monotonicity_preserving",
            "inner_solver_settings"          : {
                "preconditioner_type"            : "amg",
                "solver_type"                    : "amgcl",
                "smoother_type"                  : "ilu0",
                "krylov_type"                    : "lgmres",
                "coarsening_type"                : "aggregation",
                "max_iteration"                  : 100,
                "provide_coordinates"            : false,
                "gmres_krylov_space_dimension"   : 100,
                "verbosity"                      : 1,
                "tolerance"                      : 1e-6,
                "scaling"                        : false,
                "block_size"                     : 1,
                "use_block_matrices_if_possible" : true,
                "coarse_enough"                  : 1000,
                "max_levels"                     : -1,
                "pre_sweeps"                     : 1,
                "post_sweeps"                    : 1,
                "use_gpgpu"                      : false
            }
 
        }  )" );
 
        // Now validate agains defaults -- this also ensures no type mismatch
        ThisParameters.ValidateAndAssignDefaults(default_parameters);
 
        mpLinearSolver = LinearSolverFactoryType().Create(ThisParameters["inner_solver_settings"]);
 
        KRATOS_CATCH("")
    }
 
    TrilinosMonotonicityPreservingSolver(const TrilinosMonotonicityPreservingSolver& Other) : BaseType(Other) {}
 
 
    ~TrilinosMonotonicityPreservingSolver() override {}
 
 
 
    TrilinosMonotonicityPreservingSolver& operator=(const TrilinosMonotonicityPreservingSolver& Other)
    {
        BaseType::operator=(Other);
        return *this;
    }
 
 
    bool AdditionalPhysicalDataIsNeeded() override
    {
        return true;
    }
 
    void ProvideAdditionalData(
        SparseMatrixType& rA,
        VectorType& rX,
        VectorType& rB,
        typename ModelPart::DofsArrayType& rdof_set,
        ModelPart& r_model_part
    ) override
    {
        std::vector<int> global_ids(rdof_set.size());
        std::vector<double> dofs_values(rdof_set.size());
 
        IndexType i = 0;
        for (const auto& dof : rdof_set) {
            const int global_id = dof.EquationId();
            global_ids[i] = global_id;
            dofs_values[i] = dof.GetSolutionStepValue();
            ++i;
        }
 
        Epetra_Map localmap(-1, global_ids.size(), global_ids.data(), 0, rA.Comm());
        Epetra_Vector dofs_local(Copy, localmap, dofs_values.data());
 
        Epetra_Import dirichlet_importer(rA.ColMap(), dofs_local.Map());
 
        // defining a temporary vector to gather all of the values needed
        Epetra_Vector dofs(rA.ColMap());
 
        // importing in the new temp vector the values
        int ierr = dofs.Import(dofs_local, dirichlet_importer, Insert);
        if (ierr != 0)
            KRATOS_ERROR << "Epetra failure found";
 
        for (int i = 0; i < rA.NumMyRows(); i++) {
            int numEntries; // number of non-zero entries
            double* vals;   // row non-zero values
            int* cols;      // column indices of row non-zero values
            rA.ExtractMyRowView(i, numEntries, vals, cols);
 
            for (int j = 0; j < numEntries; j++) {
                const double value = vals[j];
                if (value > 0.0) {
                    const int row_gid = rA.RowMap().GID(i);
                    const int col_gid = rA.ColMap().GID(cols[j]);
                    if (row_gid > col_gid) {
                        Matrix LHS(2, 2);
                        LHS(0, 0) = value;
                        LHS(0, 1) = -value;
                        LHS(1, 0) = -value;
                        LHS(1, 1) = value;
                        Vector dofs_sol(2);
                        int row_lid = localmap.LID(row_gid);
                        dofs_sol[0] = dofs_local[row_lid];
                        dofs_sol[1] = dofs[cols[j]];
                        Vector RHS = -prod(LHS, dofs_sol);
                        Element::EquationIdVectorType equations_ids = {static_cast<std::size_t>(row_gid), static_cast<std::size_t>(col_gid)};
                        TSparseSpaceType::AssembleLHS(rA, LHS, equations_ids);
                        TSparseSpaceType::AssembleRHS(rB, RHS, equations_ids);
 
                    }
                }
            }
        }
 
        rA.GlobalAssemble();
        rB.GlobalAssemble();
 
        if (mpLinearSolver->AdditionalPhysicalDataIsNeeded()) {
            mpLinearSolver->ProvideAdditionalData(rA,rX,rB,rdof_set,r_model_part);
        }
    }
 
    void InitializeSolutionStep (SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        mpLinearSolver->InitializeSolutionStep(rA,rX,rB);
    }
 
    void FinalizeSolutionStep (SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        mpLinearSolver->FinalizeSolutionStep(rA,rX,rB);
    }
 
    void Clear() override
    {
        mpLinearSolver->Clear();
    }
 
    bool Solve(SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        return mpLinearSolver->Solve(rA,rX,rB);
    }
 
 
 
 
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "Composite Linear Solver. Uses internally the following linear solver " << mpLinearSolver->Info();
        return  buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        BaseType::PrintData(rOStream);
    }
 
 
 
 
 
 
private:
 
 
    typename LinearSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>::Pointer mpLinearSolver;
 
 
 
 
 
 
 
 
 
 
}; // Class TrilinosMonotonicityPreservingSolver
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TPreconditionerType,
         class TReordererType>
inline std::istream& operator >> (std::istream& IStream,
                                  TrilinosMonotonicityPreservingSolver<TSparseSpaceType, TDenseSpaceType,
                                  TReordererType>& rThis)
{
    return IStream;
}
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TPreconditionerType,
         class TReordererType>
inline std::ostream& operator << (std::ostream& OStream,
                                  const TrilinosMonotonicityPreservingSolver<TSparseSpaceType, TDenseSpaceType,
                                  TReordererType>& rThis)
{
    rThis.PrintInfo(OStream);
    OStream << std::endl;
    rThis.PrintData(OStream);
 
    return OStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_TRILINOS_MONOTONICITY_PRESERVING_SOLVER_H_INCLUDED  defined