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_MONOTONICITY_PRESERVING_SOLVER_H_INCLUDED )
#define  KRATOS_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 MonotonicityPreservingSolver
    : public LinearSolver<TSparseSpaceType, TDenseSpaceType,  TReordererType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(MonotonicityPreservingSolver);
 
    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;
 
 
    MonotonicityPreservingSolver()
    {
    }
 
    MonotonicityPreservingSolver(
        typename BaseType::Pointer pLinearSolver
        ) : BaseType (),
            mpLinearSolver(pLinearSolver)
    {
    }
 
    MonotonicityPreservingSolver(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("")
    }
 
    MonotonicityPreservingSolver(const MonotonicityPreservingSolver& Other) : BaseType(Other) {}
 
 
    ~MonotonicityPreservingSolver() override {}
 
 
 
    MonotonicityPreservingSolver& operator=(const MonotonicityPreservingSolver& 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
    {
        Vector dofs_values = ZeroVector(rdof_set.size());
 
        block_for_each(rdof_set, [&](Dof<double>& rDof){
            const std::size_t id = rDof.EquationId();
            dofs_values[id] = rDof.GetSolutionStepValue();
        });
        double *values_vector = rA.value_data().begin();
        std::size_t *index1_vector = rA.index1_data().begin();
        std::size_t *index2_vector = rA.index2_data().begin();
 
        IndexPartition<std::size_t>(rA.size1()).for_each(
            [&](std::size_t i)
            {
                for (std::size_t k = index1_vector[i]; k < index1_vector[i + 1]; k++) {
                    const double value = values_vector[k];
                    if (value > 0.0) {
                        const auto j = index2_vector[k];
                        if (j > i) {
                            rA(i,j) -= value;
                            rA(j,i) -= value;
                            // Values conflicting with other threads
                            auto& r_aii = rA(i,i).ref();
                            AtomicAdd(r_aii, value);
                            auto& r_ajj = rA(j,j).ref();
                            AtomicAdd(r_ajj, value);
                            auto& r_bi = rB[i];
                            AtomicAdd(r_bi, value*dofs_values[j] - value*dofs_values[i]);
                            auto& r_bj = rB[j];
                            AtomicAdd(r_bj, value*dofs_values[i] - value*dofs_values[j]);
                        }
                    }
                }
            }
        );
 
        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 MonotonicityPreservingSolver
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TPreconditionerType,
         class TReordererType>
inline std::istream& operator >> (std::istream& IStream,
                                  MonotonicityPreservingSolver<TSparseSpaceType, TDenseSpaceType,
                                  TReordererType>& rThis)
{
    return IStream;
}
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TPreconditionerType,
         class TReordererType>
inline std::ostream& operator << (std::ostream& OStream,
                                  const MonotonicityPreservingSolver<TSparseSpaceType, TDenseSpaceType,
                                  TReordererType>& rThis)
{
    rThis.PrintInfo(OStream);
    OStream << std::endl;
    rThis.PrintData(OStream);
 
    return OStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_MONOTONICITY_PRESERVING_SOLVER_H_INCLUDED  defined