sensitivity_builder_scheme.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Suneth Warnakulasuriya
//
 
#if !defined(KRATOS_SENSITIVITY_BUILDER_SCHEME_H_INCLUDED)
#define KRATOS_SENSITIVITY_BUILDER_SCHEME_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "containers/global_pointers_vector.h"
#include "includes/define.h"
#include "includes/model_part.h"
#include "includes/parallel_environment.h"
#include "response_functions/adjoint_response_function.h"
#include "utilities/global_pointer_utilities.h"
#include "utilities/openmp_utils.h"
#include "utilities/parallel_utilities.h"
 
namespace Kratos
{
 
class SensitivityBuilderScheme
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(SensitivityBuilderScheme);
 
    using NodeType = ModelPart::NodeType;
 
    using ConditionType = ModelPart::ConditionType;
 
    using ElementType = ModelPart::ElementType;
 
 
    SensitivityBuilderScheme()
        : mRank(ParallelEnvironment::GetDefaultDataCommunicator().Rank())
    {
        // Allocate auxiliary memory.
        // This needs to be done in the constructor because, this scheme
        // is used to calculate sensitivities w.r.t. element quantities
        // and in there we don't usually pass the model part. Hence, we
        // can not call SensitivityBuilderScheme::Initialize
        // method.
        const int number_of_threads = ParallelUtilities::GetNumThreads();
        mSensitivityMatrices.resize(number_of_threads);
        mAdjointVectors.resize(number_of_threads);
        mPartialSensitivity.resize(number_of_threads);
    }
 
    virtual ~SensitivityBuilderScheme() = default;
 
 
    virtual int Check(
        const ModelPart& rModelPart,
        const ModelPart& rSensitivityModelPart) const
    {
        return 0;
    }
 
    virtual void Initialize(
        ModelPart& rModelPart,
        ModelPart& rSensitivityModelPart,
        AdjointResponseFunction& rResponseFunction)
    {
        const auto& r_communicator = rModelPart.GetCommunicator();
        const auto& r_nodes = rModelPart.Nodes();
        const int number_of_nodes = r_nodes.size();
 
        std::vector<int> indices;
        indices.resize(number_of_nodes);
        IndexPartition<int>(number_of_nodes).for_each([&](const int Index) {
            indices[Index] = (r_nodes.begin() + Index)->Id();
        });
 
        const auto& r_data_communicator = r_communicator.GetDataCommunicator();
        mGlobalPointerNodalMap = GlobalPointerUtilities::RetrieveGlobalIndexedPointersMap(
            r_nodes, indices, r_data_communicator);
    }
 
    virtual void InitializeSolutionStep(
        ModelPart& rModelPart,
        ModelPart& rSensitivityModelPart,
        AdjointResponseFunction& rResponseFunction)
    {}
 
    virtual void FinalizeSolutionStep(
        ModelPart& rModelPart,
        ModelPart& rSensitivityModelPart,
        AdjointResponseFunction& rResponseFunction)
    {}
 
    virtual void Finalize(
        ModelPart& rModelPart,
        ModelPart& rSensitivityModelPart,
        AdjointResponseFunction& rResponseFunction)
    {}
 
    virtual void Update(
        ModelPart& rModelPart,
        ModelPart& rSensitivityModelPart,
        AdjointResponseFunction& rResponseFunction)
    {}
 
    virtual void Clear()
    {
        mSensitivityMatrices.clear();
        mAdjointVectors.clear();
        mPartialSensitivity.clear();
        mGlobalPointerNodalMap.clear();
    }
 
    virtual void CalculateSensitivity(
        ElementType& rCurrentElement,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<NodeType>& rGPSensitivityVector,
        const Variable<double>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentElement, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ElementType& rCurrentElement,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ConditionType>& rGPSensitivityVector,
        const Variable<double>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {}
 
    virtual void CalculateSensitivity(
        ElementType& rCurrentElement,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ElementType>& rGPSensitivityVector,
        const Variable<double>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentElement, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ConditionType& rCurrentCondition,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<NodeType>& rGPSensitivityVector,
        const Variable<double>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentCondition, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ConditionType& rCurrentCondition,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ConditionType>& rGPSensitivityVector,
        const Variable<double>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentCondition, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ConditionType& rCurrentCondition,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ElementType>& rGPSensitivityVector,
        const Variable<double>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {}
 
    virtual void CalculateSensitivity(
        ElementType& rCurrentElement,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<NodeType>& rGPSensitivityVector,
        const Variable<array_1d<double, 3>>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentElement, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ElementType& rCurrentElement,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ConditionType>& rGPSensitivityVector,
        const Variable<array_1d<double, 3>>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {}
 
    virtual void CalculateSensitivity(
        ElementType& rCurrentElement,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ElementType>& rGPSensitivityVector,
        const Variable<array_1d<double, 3>>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentElement, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ConditionType& rCurrentCondition,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<NodeType>& rGPSensitivityVector,
        const Variable<array_1d<double, 3>>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentCondition, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ConditionType& rCurrentCondition,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ConditionType>& rGPSensitivityVector,
        const Variable<array_1d<double, 3>>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        CalculateLocalSensitivityAndGlobalPointersVector(
            rCurrentCondition, rResponseFunction, rSensitivity,
            rGPSensitivityVector, rVariable, rCurrentProcessInfo);
    }
 
    virtual void CalculateSensitivity(
        ConditionType& rCurrentCondition,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        GlobalPointersVector<ElementType>& rGPSensitivityVector,
        const Variable<array_1d<double, 3>>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {}
 
 
    virtual std::string Info() const
    {
        return "SensitivityBuilderScheme";
    }
 
    virtual void PrintInfo(std::ostream& rOStream) const
    {
        rOStream << Info();
    }
 
    virtual void PrintData(std::ostream& rOStream) const
    {
        rOStream << Info();
    }
 
 
protected:
 
    std::vector<Matrix> mSensitivityMatrices;
    std::vector<Vector> mAdjointVectors;
    std::vector<Vector> mPartialSensitivity;
 
    std::unordered_map<int, GlobalPointer<ModelPart::NodeType>> mGlobalPointerNodalMap;
    const int mRank;
 
private:
 
    template <typename TEntityType, typename TDerivativeEntityType, typename TDataType>
    void CalculateLocalSensitivityAndGlobalPointersVector(
        TEntityType& rEntity,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivityVector,
        GlobalPointersVector<TDerivativeEntityType>& rGPSensitivityVector,
        const Variable<TDataType>& rVariable,
        const ProcessInfo& rProcessInfo)
    {
        KRATOS_TRY;
 
        this->CalculateLocalSensitivity(rEntity, rResponseFunction,
                                        rSensitivityVector, rVariable, rProcessInfo);
 
        if (rGPSensitivityVector.size() != 1) {
            rGPSensitivityVector.resize(1);
        }
 
        rGPSensitivityVector(0) = GlobalPointer<TDerivativeEntityType>(&rEntity, mRank);
 
        KRATOS_CATCH("");
    }
 
    template <typename TEntityType, typename TDataType>
    void CalculateLocalSensitivityAndGlobalPointersVector(
        TEntityType& rEntity,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivityVector,
        GlobalPointersVector<NodeType>& rGPSensitivityVector,
        const Variable<TDataType>& rVariable,
        const ProcessInfo& rProcessInfo)
    {
        KRATOS_TRY;
 
        this->CalculateLocalSensitivity(rEntity, rResponseFunction,
                                        rSensitivityVector, rVariable, rProcessInfo);
 
        auto& r_geometry = rEntity.GetGeometry();
        if (rGPSensitivityVector.size() != r_geometry.PointsNumber()) {
            rGPSensitivityVector.resize(r_geometry.PointsNumber());
        }
 
        for (unsigned int i = 0; i < r_geometry.PointsNumber(); ++i) {
            rGPSensitivityVector(i) = mGlobalPointerNodalMap[r_geometry[i].Id()];
        }
 
        KRATOS_CATCH("");
    }
 
    template <class TEntityType, class TDataType>
    void CalculateLocalSensitivity(
        TEntityType& rCurrentEntity,
        AdjointResponseFunction& rResponseFunction,
        Vector& rSensitivity,
        const Variable<TDataType>& rVariable,
        const ProcessInfo& rCurrentProcessInfo)
    {
        KRATOS_TRY;
 
        const auto k = OpenMPUtils::ThisThread();
 
        rCurrentEntity.CalculateSensitivityMatrix(rVariable, mSensitivityMatrices[k], rCurrentProcessInfo);
        rCurrentEntity.GetValuesVector(mAdjointVectors[k]);
 
        KRATOS_ERROR_IF(mAdjointVectors[k].size() != mSensitivityMatrices[k].size2())
            << "mAdjointVectors.size(): " << mAdjointVectors[k].size()
            << " incompatible with mSensitivityMatrices[k].size1(): "
            << mSensitivityMatrices[k].size2() << ". Variable: " << rVariable << std::endl;
 
        rResponseFunction.CalculatePartialSensitivity(
            rCurrentEntity, rVariable, mSensitivityMatrices[k], mPartialSensitivity[k], rCurrentProcessInfo);
 
        KRATOS_ERROR_IF(mPartialSensitivity[k].size() != mSensitivityMatrices[k].size1())
            << "mPartialSensitivity.size(): " << mPartialSensitivity[k].size()
            << " incompatible with mSensitivityMatrices.size1(): "
            << mSensitivityMatrices[k].size1() << ". Variable: " << rVariable << std::endl;
 
        if (rSensitivity.size() != mSensitivityMatrices[k].size1()) {
            rSensitivity.resize(mSensitivityMatrices[k].size1(), false);
        }
 
        noalias(rSensitivity) = prod(mSensitivityMatrices[k], mAdjointVectors[k]) + mPartialSensitivity[k];
 
        KRATOS_CATCH("");
    }
 
 
}; /* Class SensitivityBuilderScheme */
 
 
} /* namespace Kratos.*/
 
#endif /* KRATOS_SENSITIVITY_BUILDER_SCHEME_H_INCLUDED defined */