diff --git a/applications/StructuralMechanicsApplication/CMakeLists.txt b/applications/StructuralMechanicsApplication/CMakeLists.txt index e905d6b51fa6..7184d6f2e7dc 100644 --- a/applications/StructuralMechanicsApplication/CMakeLists.txt +++ b/applications/StructuralMechanicsApplication/CMakeLists.txt @@ -12,6 +12,7 @@ file(GLOB_RECURSE KRATOS_STRUCTURAL_MECHANICS_APPLICATION_CORE ${CMAKE_CURRENT_SOURCE_DIR}/structural_mechanics_application.cpp ${CMAKE_CURRENT_SOURCE_DIR}/structural_mechanics_application_variables.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_conditions/*.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/custom_constraints/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_elements/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/custom_io/*.cpp diff --git a/applications/StructuralMechanicsApplication/custom_constraints/link_constraint.cpp b/applications/StructuralMechanicsApplication/custom_constraints/link_constraint.cpp new file mode 100644 index 000000000000..c6f486163b0e --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_constraints/link_constraint.cpp @@ -0,0 +1,323 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "custom_constraints/link_constraint.hpp" // LinkConstraint +#include "includes/define.h" +#include "includes/process_info.h" +#include "includes/variables.h" // DISPLACEMENT_X, DISPLACEMENT_Y, DISPLACEMENT_Z +#include "includes/serializer.h" // Serializer +#include "includes/variables.h" // CONSTITUTIVE_MATRIX, INTERNAL_FORCES_VECTOR + +// STL includes +#include // std::max_element, std::equal, std::transform + + +namespace Kratos { + + +struct LinkConstraint::Impl +{ + using ValueType = double; + + static void ComputeRelationMatrix(LinkConstraint::MatrixType& rRelationMatrix, + LinkConstraint::MatrixType& rHessian, + LinkConstraint::VectorType& rConstraintGaps, + const std::array& rLastPositions, + const std::array& rLastDisplacements, + const unsigned Dimensions) + { + // Set output sizes. + rRelationMatrix = ZeroMatrix(1, 2 * Dimensions); + rHessian = ZeroMatrix(2 * Dimensions, 2 * Dimensions); + rConstraintGaps = ZeroVector(1); + + // Compute constraint violation. + double current_norm_square = 0.0; + for (std::size_t i_component=0ul; i_component*,3> displacement_components { + &DISPLACEMENT_X, + &DISPLACEMENT_Y, + &DISPLACEMENT_Z + }; + + std::array nodes {&rFirst, &rSecond}; + dofs.resize(2 * Dimensions); + + // Collect dofs from nodes in the following order: + // {u^0_x, ..., u^0_w, u^1_x, ..., u^1_w} + // Where w is the axis with the highest index for the provided dimensions + // (e.g.: z-axis in 3 dimensions). + for (std::size_t i_dimension=0ul; i_dimensionKey(); + + for (std::size_t i_node=0ul; i_node<2; ++i_node) { + Node& r_node = *nodes[i_node]; + const auto it_dof = std::find_if(r_node.GetDofs().begin(), + r_node.GetDofs().end(), + [component_id](auto& rp_dof){ + return rp_dof->GetVariable().Key() == component_id; + }); + KRATOS_ERROR_IF(it_dof == r_node.GetDofs().end()) + << "cannot find DoF " << displacement_components[i_dimension]->Name() + << " in node " << r_node.Id(); + dofs[i_node * Dimensions + i_dimension] = it_dof->get(); + } // for i_node in range(2) + } // for i_dimension in range(Dimensions) + KRATOS_CATCH("") + + return dofs; + } + + + void FetchNodePositions(std::array& rPositions, + std::array& rDisplacements) const { + KRATOS_TRY + + const std::array*,3> displacement_components { + &DISPLACEMENT_X, + &DISPLACEMENT_Y, + &DISPLACEMENT_Z + }; + + for (std::size_t i_component=0u; i_componentCoordinates()[i_component]; + rDisplacements[i_dof] = mNodePair[i_node]->GetSolutionStepValue(*displacement_components[i_component]); + } // for i_node in range(2) + } // for i_component in range(mDimensions) + + KRATOS_CATCH("") + } + + + std::size_t mDimensions; + + bool mIsMeshMoved; + + /// @details MasterSlaveConstraint::GetSlaveDofsVector and MasterSlaveConstraint::GetMasterDofsVector + /// require arrays of mutable Dof pointers, which are only obtainable from mutable nodes, + /// so the nodes' pointers stored here cannot be immutable. Risky business. + std::array mNodePair; + + std::array mLastPositions; + + std::array mLastDisplacements; +}; // struct LinkConstraint::Impl + + +LinkConstraint::LinkConstraint(const IndexType Id, + Node& rFirst, + Node& rSecond, + const std::size_t Dimensions, + bool IsMeshMoved) + : MultifreedomConstraint(/*Id : */Id, + /*rDofs : */Impl::CollectDofs(rFirst, rSecond, Dimensions), + /*rConstraintLabels: */std::vector {Id}), + mpImpl(new Impl{/*mDimensions : */Dimensions, + /*mIsMeshMoved : */IsMeshMoved, + /*mNodePair : */{&rFirst, &rSecond}, + /*mLastPositions : */{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + /*mLastDisplacements : */{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}) +{ + this->Data().GetValue(CONSTITUTIVE_MATRIX) = ZeroMatrix(0); + this->Data().GetValue(INTERNAL_FORCES_VECTOR) = ZeroVector(0); +} + + +LinkConstraint::~LinkConstraint() = default; + + +MasterSlaveConstraint::Pointer LinkConstraint::Clone(IndexType NewId) const +{ + KRATOS_TRY + auto p_other = LinkConstraint::Pointer( + new LinkConstraint(NewId, + *mpImpl->mNodePair.front(), + *mpImpl->mNodePair.back(), + mpImpl->mDimensions, + mpImpl->mIsMeshMoved) + ); + return p_other; + KRATOS_CATCH("") +} + + +void LinkConstraint::Initialize([[maybe_unused]] const ProcessInfo&) +{ +} + + +void LinkConstraint::InitializeSolutionStep(const ProcessInfo&) +{ + KRATOS_TRY + // Store displacements before nonlinear iterations begin. + // At this point, displacements are at their last converged + // state, when the constraint is supposed to be still satisfied. + mpImpl->FetchNodePositions(mpImpl->mLastPositions, + mpImpl->mLastDisplacements); + + if (!mpImpl->mIsMeshMoved) + for (std::size_t i_component=0ul; i_componentmLastPositions.size(); ++i_component) + mpImpl->mLastPositions[i_component] += mpImpl->mLastDisplacements[i_component]; + + Impl::ComputeRelationMatrix(this->Data().GetValue(CONSTITUTIVE_MATRIX), + this->GetHessian(), + this->Data().GetValue(INTERNAL_FORCES_VECTOR), + mpImpl->mLastPositions, + mpImpl->mLastDisplacements, + mpImpl->mDimensions); + KRATOS_CATCH("") +} + + +void LinkConstraint::InitializeNonLinearIteration(const ProcessInfo&) +{ + KRATOS_TRY + [[maybe_unused]] std::array dummy; + mpImpl->FetchNodePositions(dummy, mpImpl->mLastDisplacements); + + Impl::ComputeRelationMatrix(this->Data().GetValue(CONSTITUTIVE_MATRIX), + this->GetHessian(), + this->Data().GetValue(INTERNAL_FORCES_VECTOR), + mpImpl->mLastPositions, + mpImpl->mLastDisplacements, + mpImpl->mDimensions); + KRATOS_CATCH("") +} + + +void LinkConstraint::FinalizeNonLinearIteration(const ProcessInfo& rProcessInfo) +{ +} + + +void LinkConstraint::FinalizeSolutionStep(const ProcessInfo&) +{ + this->Data().GetValue(CONSTITUTIVE_MATRIX).clear(); + this->Data().GetValue(INTERNAL_FORCES_VECTOR).clear(); +} + + +void LinkConstraint::Finalize(const ProcessInfo&) +{ +} + + +void LinkConstraint::CalculateLocalSystem(MatrixType& rRelationMatrix, + VectorType& rConstraintGaps, + [[maybe_unused]] const ProcessInfo& rProcessInfo) const +{ + rRelationMatrix = this->GetData().GetValue(CONSTITUTIVE_MATRIX); + rConstraintGaps = this->GetData().GetValue(INTERNAL_FORCES_VECTOR); +} + + +int LinkConstraint::Check(const ProcessInfo& rProcessInfo) const +{ + // Check whatever the base class checks. + MasterSlaveConstraint::Check(rProcessInfo); + + // Restrict dimensions to 1, 2, or 3. + KRATOS_ERROR_IF_NOT(0 < mpImpl->mDimensions || mpImpl->mDimensions <= 3) << "unsupported dimensions (" << mpImpl->mDimensions << ")"; + + // Make sure that the nodes have the necessary Dofs. + const std::array*,3> displacement_components { + &DISPLACEMENT_X, + &DISPLACEMENT_Y, + &DISPLACEMENT_Z + }; + + for (unsigned i_component=0u; i_componentmDimensions; ++i_component) { + for (unsigned i_node=0u; i_node<2; ++i_node) { + KRATOS_ERROR_IF_NOT(mpImpl->mNodePair[i_node]->HasDofFor(*displacement_components[i_component])) + << "node " << mpImpl->mNodePair[i_node]->Id() + << " has no Dof for " << displacement_components[i_component]->Name(); + } // for i_node in range(2) + } // for i_component in range(mDimensions) + + // Check for overlapping nodes. + constexpr Impl::ValueType tolerance = 1e-16; + KRATOS_ERROR_IF(norm_2(mpImpl->mNodePair.front()->Coordinates() - mpImpl->mNodePair.back()->Coordinates()) < tolerance) + << "coincident nodes in LinkConstraint"; + + return 0; +} + + +void LinkConstraint::save(Serializer& rSerializer) const +{ + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, MultifreedomConstraint); + rSerializer.save("Dimensions", mpImpl->mDimensions); + rSerializer.save("Nodes", mpImpl->mNodePair); + // The rest of the private members should be computed/cached after restarting. +} + + +void LinkConstraint::load(Serializer& rDeserializer) +{ + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rDeserializer, MultifreedomConstraint); + rDeserializer.load("Dimensions", mpImpl->mDimensions); + rDeserializer.load("Nodes", mpImpl->mNodePair); +} + + +std::ostream& operator<<(std::ostream& rStream, const LinkConstraint& rInstance) +{ + return rStream << "LinkConstraint between nodes " + << rInstance.mpImpl->mNodePair.front()->Id() + << " and " + << rInstance.mpImpl->mNodePair.back()->Id(); +} + + +} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_constraints/link_constraint.hpp b/applications/StructuralMechanicsApplication/custom_constraints/link_constraint.hpp new file mode 100644 index 000000000000..f0283afdf526 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_constraints/link_constraint.hpp @@ -0,0 +1,140 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.hpp" // MultifreedomConstraint +#include "includes/code_location.h" // KRATOS_CODE_LOCATION + +// STL includes +#include // std::unique_ptr +#include // std::ostream + + +namespace Kratos { + + +/** @brief A constraint enforcing the distance between two @ref Node "nodes" to remain constant. + * @details Let @f$ n^0_i @f$ and @f$ n^1_i @f$ denote the coordinates of the initial positions of @ref Node "node" + * @f$ n^0 @f$ and @f$ n^1 @f$ that get displaced by @f$ u^0_i @f$ and @f$ u^1_i @f$ respectively. + * @p LinkConstraint then represents the following constraint equation: + * @f[ + * \sqrt{ + * \sum_i \left( (n^0_i + u^0_i) - (n^1_i + u^1_i) \right)^2 + * } + * - + * \sqrt{ + * \sum_i \left( n^0_i - n^1_i \right)^2 + * } + * = 0 + * @f] + * + * @note @p LinkConstraint is nonlinear and as such, should not be imposed via master-slave elimination (the slave + * DoF's coefficient may vanish). Consider using the @ref AugmentedLagrangeConstraintAssembler "penalty method", + * the method of Lagrange multipliers, or @ref AugmentedLagrangeConstraintAssembler "augmented lagrange" multipliers. + */ +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinkConstraint final + : public MultifreedomConstraint +{ +public: + KRATOS_CLASS_POINTER_DEFINITION(LinkConstraint); + + LinkConstraint() noexcept = default; + + LinkConstraint(const IndexType Id, + Node& rFirst, + Node& rSecond, + const std::size_t Dimensions, + bool IsMeshMoved); + + /// @internal + ~LinkConstraint(); + + MasterSlaveConstraint::Pointer Clone(IndexType NewId) const override; + + void Initialize(const ProcessInfo& rProcessInfo) override; + + void InitializeSolutionStep(const ProcessInfo& rProcessInfo) override; + + void InitializeNonLinearIteration(const ProcessInfo& rProcessInfo) override; + + void FinalizeNonLinearIteration(const ProcessInfo& rProcessInfo) override; + + void FinalizeSolutionStep(const ProcessInfo& rProcessInfo) override; + + void Finalize(const ProcessInfo& rProcessInfo) override; + + void CalculateLocalSystem(MatrixType& rRelationMatrix, + VectorType& rConstraintGaps, + const ProcessInfo& rProcessInfo) const override; + + int Check(const ProcessInfo& rProcessInfo) const override; + + friend std::ostream& operator<<(std::ostream& rStream, const LinkConstraint& rInstance); + + /// @name Unsupported Interface + /// @{ + + MasterSlaveConstraint::Pointer Create(IndexType, + DofPointerVectorType&, + DofPointerVectorType&, + const MatrixType&, + const VectorType&) const override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported by LinkConstraint.";} + + MasterSlaveConstraint::Pointer Create(IndexType, + NodeType&, + const VariableType&, + NodeType&, + const VariableType&, + const double, + const double) const override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported by LinkConstraint.";} + + void SetDofList(const DofPointerVectorType&, + const DofPointerVectorType&, + const ProcessInfo&) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported by LinkConstraint.";} + + void SetSlaveDofsVector(const DofPointerVectorType&) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported by LinkConstraint.";} + + void SetMasterDofsVector(const DofPointerVectorType&) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported by LinkConstraint.";} + + void SetLocalSystem(const MatrixType&, + const VectorType&, + const ProcessInfo&) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported by LinkConstraint.";} + + /// @} + +private: + LinkConstraint(LinkConstraint&& rRhs) = delete; + + LinkConstraint(const LinkConstraint& rRhs) = delete; + + LinkConstraint& operator=(LinkConstraint&& rRhs) = delete; + + LinkConstraint& operator=(const LinkConstraint& rRhs) = delete; + + friend class Serializer; + + void save(Serializer& rSerializer) const override; + + void load(Serializer& rDeserializer) override; + + struct Impl; + std::unique_ptr mpImpl; +}; // class LinkConstraint + + +} // namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_constraints_to_python.cpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_constraints_to_python.cpp new file mode 100644 index 000000000000..9e8d0de64a91 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_constraints_to_python.cpp @@ -0,0 +1,34 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Mate Kelemen +// + +// Project includes +#include "custom_python/add_custom_constraints_to_python.hpp" // AddCustomConstraintsToPython +#include "custom_constraints/link_constraint.hpp" // LinkConstraint + + +namespace Kratos::Python { + + +void AddCustomConstraintsToPython(pybind11::module& rModule) +{ + pybind11::class_(rModule, "LinkConstraint") + .def(pybind11::init(), + pybind11::arg("Id"), + pybind11::arg("FirstNode"), + pybind11::arg("SecondNode"), + pybind11::arg("Dimensions"), + pybind11::arg("IsMeshMoved")) + .doc() = "A constraint enforcing the distance between two nodes to remain constant." + ; +} + + +} // namespace Kratos::Python diff --git a/applications/StructuralMechanicsApplication/custom_python/add_custom_constraints_to_python.hpp b/applications/StructuralMechanicsApplication/custom_python/add_custom_constraints_to_python.hpp new file mode 100644 index 000000000000..6de970b5899e --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_python/add_custom_constraints_to_python.hpp @@ -0,0 +1,20 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Mate Kelemen +// + +#pragma once + +// External includes +#include + +namespace Kratos::Python +{ + void AddCustomConstraintsToPython(pybind11::module& rModule); +} // namespace Kratos::Python. diff --git a/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp b/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp index 3b1bc79ad035..3b6fae63902d 100644 --- a/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp +++ b/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp @@ -23,6 +23,7 @@ #include "custom_python/add_custom_processes_to_python.h" #include "custom_python/add_custom_utilities_to_python.h" #include "custom_python/add_custom_constitutive_laws_to_python.h" +#include "custom_python/add_custom_constraints_to_python.hpp" #include "custom_python/add_custom_response_functions_to_python.h" namespace Kratos::Python { @@ -41,6 +42,7 @@ PYBIND11_MODULE(KratosStructuralMechanicsApplication,m) AddCustomProcessesToPython(m); AddCustomUtilitiesToPython(m); AddCustomConstitutiveLawsToPython(m); + AddCustomConstraintsToPython(m); AddCustomResponseFunctionUtilitiesToPython(m); py::class_,VariableData >(m,"ShellCrossSectionVariable"); diff --git a/applications/StructuralMechanicsApplication/python_scripts/link_constraint_process.py b/applications/StructuralMechanicsApplication/python_scripts/link_constraint_process.py new file mode 100644 index 000000000000..e166821aa819 --- /dev/null +++ b/applications/StructuralMechanicsApplication/python_scripts/link_constraint_process.py @@ -0,0 +1,81 @@ +# --- Kratos Imports --- +import KratosMultiphysics +import KratosMultiphysics.StructuralMechanicsApplication + +# --- STD Imports --- +import typing + + +class LinkConstraintProcess(KratosMultiphysics.Process): + """@brief Create and manage @ref LinkConstraint "constraints" between pairs of @ref Node "nodes" that force them to keep the distance between each other constant.""" + + + def __init__(self, + model: KratosMultiphysics.Model, + parameters: KratosMultiphysics.Parameters): + super().__init__() + parameters.ValidateAndAssignDefaults(self.GetDefaultParameters()) + self.__model_part: KratosMultiphysics.ModelPart = model.GetModelPart(parameters["model_part_name"].GetString()) + self.__node_pairs: "list[tuple[int,int]]" = [] + self.__dimensions: int = parameters["dimensions"].GetInt() + self.__is_mesh_moved: bool = parameters["move_mesh_flag"].GetBool() + self.__interval: KratosMultiphysics.IntervalUtility = KratosMultiphysics.IntervalUtility(parameters) + self.__is_active: bool = False + self.__constraints: "list[KratosMultiphysics.StructuralMechanicsApplication.LinkConstraint]" = [] + + pair: KratosMultiphysics.Parameters + for pair in parameters["node_pairs"].values(): + self.__node_pairs.append((pair[0].GetInt(), pair[1].GetInt())) + + + @classmethod + def GetDefaultParameters(cls: "typing.Type[LinkConstraintProcess]") -> KratosMultiphysics.Parameters: + return KratosMultiphysics.Parameters(R"""{ + "model_part_name" : "", + "node_pairs" : [], + "dimensions" : 3, + "move_mesh_flag" : false, + "interval" : [0, "End"] + }""") + + + def ExecuteBeforeSolutionLoop(self) -> None: + self.__is_active = self.__interval.IsInInterval(self.__model_part.ProcessInfo[KratosMultiphysics.TIME]) + + id: int = self.__GetLastConstraintId() + 1 + for id_first, id_second in self.__node_pairs: + constraint = KratosMultiphysics.StructuralMechanicsApplication.LinkConstraint( + id, + self.__model_part.GetNode(id_first), + self.__model_part.GetNode(id_second), + self.__dimensions, + self.__is_mesh_moved) + + constraint.Set(KratosMultiphysics.ACTIVE, self.__is_active) + + self.__constraints.append(constraint) + self.__model_part.AddMasterSlaveConstraint(constraint) + id += 1 + + + def ExecuteInitializeSolutionStep(self) -> None: + is_active: bool = self.__interval.IsInInterval(self.__model_part.ProcessInfo[KratosMultiphysics.TIME]) + if is_active != self.__is_active: + self.__is_active = is_active + for constraint in self.__constraints: + constraint.Set(KratosMultiphysics.ACTIVE, self.__is_active) + + + def __GetLastConstraintId(self) -> int: + constraint_count: int = len(self.__model_part.MasterSlaveConstraints) + last_constraint_id: int + if constraint_count: + last_constraint_id = self.__model_part.MasterSlaveConstraints[constraint_count - 1] + else: + last_constraint_id = 0 + return self.__model_part.GetCommunicator().GetDataCommunicator().MaxAll(last_constraint_id) + + +def Factory(parameters: KratosMultiphysics.Parameters, + model: KratosMultiphysics.Model) -> LinkConstraintProcess: + return LinkConstraintProcess(model, parameters["Parameters"]) diff --git a/applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_solver.py b/applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_solver.py index cef48f81e83f..ddab582916e8 100755 --- a/applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_solver.py +++ b/applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_solver.py @@ -473,11 +473,15 @@ def _CreateBuilderAndSolver(self) -> KratosMultiphysics.BuilderAndSolver: return KratosMultiphysics.ResidualBasedEliminationBuilderAndSolverWithConstraints(linear_solver) else: return KratosMultiphysics.ResidualBasedEliminationBuilderAndSolver(linear_solver) + elif builder_and_solver_name == "p_multigrid": + parameters = self.settings["builder_and_solver_settings"]["advanced_settings"] + return KratosMultiphysics.PMultigridBuilderAndSolver(linear_solver, parameters) else: message: str = f"Invalid type for builder and solver '{builder_and_solver_name}'. Options are:\n" message += "'block'\n" message += "'block_lagrange'\n" message += "'elimination'" + message += "'p_multigrid'" raise ValueError(message) def _CreateScheme(self): diff --git a/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d.json b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d.json new file mode 100644 index 000000000000..4344107c8707 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d.json @@ -0,0 +1,92 @@ +{ + "analysis_stage": "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis", + "problem_data": { + "problem_name": "test_link_constraint_2d", + "parallel_type": "OpenMP", + "echo_level": 1, + "start_time": 0, + "end_time": 1 + }, + "solver_settings": { + "time_stepping": { + "time_step": 0.5 + }, + "solver_type": "Static", + "model_part_name": "root", + "domain_size": 2, + "echo_level": 0, + "analysis_type": "non_linear", + "reform_dofs_at_each_step" : false, + "move_mesh_flag" : false, + "max_iteration" : 1e2, + "model_import_settings": { + "input_type": "mdpa", + "input_filename" : "link_constraint_2d" + }, + "material_import_settings": { + "materials_filename": "link_constraint_2d_materials.json" + }, + "convergence_criterion": "displacement_criterion", + "displacement_relative_tolerance": 1e-4, + "displacement_absolute_tolerance": 1e-9, + "builder_and_solver_settings" : { + "type" : "p_multigrid", + "advanced_settings" : { + "max_iterations" : 1, + "verbosity" : 1, + "linear_solver_settings" : {"solver_type" : "skyline_lu_factorization"}, + "constraint_imposition_settings" : { + "method" : "augmented_lagrange", + "penalty_factor" : "1e6 * max", + "max_iterations" : 2e1, + "tolerance" : 1e-9, + "verbosity" : 1 + } + } + } + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "model_part_name": "root.dirichlet", + "variable_name": "DISPLACEMENT", + "interval": [0, "End"], + "constrained": [true, true, false], + "value": [0, 0, 0] + } + }, + { + "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", + "python_module" : "link_constraint_process", + "process_name" : "LinkConstraintProcess", + "Parameters" : { + "model_part_name" : "root", + "node_pairs" : [[2, 3]], + "dimensions" : 2, + "move_mesh_flag" : false, + "interval" : [0, 0.6] + } + } + ], + "loads_process_list": [ + { + "python_module": "assign_vector_by_direction_to_condition_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorByDirectionToConditionProcess", + "Parameters": { + "model_part_name": "root.neumann", + "variable_name": "POINT_LOAD", + "interval": [ + 0, "End" + ], + "modulus": "t", + "direction": [-2, 7, 0] + } + } + ] + } +} diff --git a/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d.mdpa b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d.mdpa new file mode 100644 index 000000000000..c476b2eff02a --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d.mdpa @@ -0,0 +1,31 @@ +Begin Properties 0 +End Properties + +Begin Nodes + 1 0 0 0 + 2 1 0 0 + 3 0 1 0 + 4 1 1 0 +End Nodes + +Begin Elements SmallDisplacementElement2D3N + 1 0 1 2 3 + 2 0 3 2 4 +End Elements + +Begin Conditions PointLoadCondition2D1N + 1 0 3 +End Conditions + +Begin SubModelPart neumann + Begin SubModelPartConditions + 1 + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart dirichlet + Begin SubModelPartNodes + 1 + 4 + End SubModelPartNodes +End SubModelPart diff --git a/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d_materials.json b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d_materials.json new file mode 100644 index 000000000000..363be074d2d9 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_2d_materials.json @@ -0,0 +1,19 @@ +{ + "properties": [ + { + "model_part_name": "root", + "properties_id": 1, + "Material": { + "constitutive_law": { + "name": "LinearElasticPlaneStress2DLaw" + }, + "Variables": { + "DENSITY": 1.0, + "YOUNG_MODULUS": 1.0, + "POISSON_RATIO": 0.3 + }, + "Tables": null + } + } + ] +} diff --git a/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d.json b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d.json new file mode 100644 index 000000000000..f34e54d7eebf --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d.json @@ -0,0 +1,93 @@ +{ + "analysis_stage": "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis", + "problem_data": { + "problem_name": "test_link_constraint_3d", + "parallel_type": "OpenMP", + "echo_level": 1, + "start_time": 0.0, + "end_time": 1.0 + }, + "solver_settings": { + "time_stepping": { + "time_step": 0.5 + }, + "solver_type": "Static", + "model_part_name": "root", + "domain_size": 3, + "echo_level": 0, + "analysis_type": "non_linear", + "move_mesh_flag" : true, + "max_iteration" : 1e1, + "model_import_settings": { + "input_type": "mdpa", + "input_filename" : "link_constraint_3d" + }, + "material_import_settings": { + "materials_filename": "link_constraint_3d_materials.json" + }, + "convergence_criterion": "displacement_criterion", + "displacement_relative_tolerance": 1e-4, + "displacement_absolute_tolerance": 1e-9, + "builder_and_solver_settings" : { + "type" : "p_multigrid", + "advanced_settings" : { + "max_iterations" : 1, + "verbosity" : 1, + "linear_solver_settings" : {"solver_type" : "skyline_lu_factorization"}, + "constraint_imposition_settings" : { + "method" : "augmented_lagrange", + "penalty_factor" : "1e6 * max", + "max_iterations" : 5e1, + "tolerance" : 1e-9, + "verbosity" : 1 + } + } + } + }, + "processes": { + "constraints_process_list": [ + { + "python_module": "assign_vector_variable_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorVariableProcess", + "Parameters": { + "model_part_name": "root.dirichlet", + "variable_name": "DISPLACEMENT", + "interval": [0.0, "End"], + "constrained": [true, true, true], + "value": [0.0, 0.0, 0.0] + } + }, + { + "kratos_module" : "KratosMultiphysics.StructuralMechanicsApplication", + "python_module" : "link_constraint_process", + "process_name" : "LinkConstraintProcess", + "Parameters" : { + "model_part_name" : "root", + "node_pairs" : [[2, 3]], + "dimensions" : 3, + "interval" : [0, 2], + "move_mesh_flag" : true + } + } + ], + "loads_process_list": [ + { + "python_module": "assign_vector_by_direction_to_condition_process", + "kratos_module": "KratosMultiphysics", + "process_name": "AssignVectorByDirectionToConditionProcess", + "Parameters": { + "model_part_name": "root.neumann", + "variable_name": "POINT_LOAD", + "interval": [ + 0.0, "End" + ], + "modulus": "3e-1 * t", + "direction": [1, 3, 2] + } + } + ], + "list_other_processes": [] + }, + "output_processes": {} +} diff --git a/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d.mdpa b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d.mdpa new file mode 100644 index 000000000000..fb896830facd --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d.mdpa @@ -0,0 +1,30 @@ +Begin Properties 0 +End Properties + +Begin Nodes + 1 0 0 0 + 2 1 0 0 + 3 0 1 0 + 4 0 0 1 +End Nodes + +Begin Elements SmallDisplacementElement3D4N + 1 0 1 2 3 4 +End Elements + +Begin Conditions PointLoadCondition3D1N + 1 0 2 +End Conditions + +Begin SubModelPart neumann + Begin SubModelPartConditions + 1 + End SubModelPartConditions +End SubModelPart + +Begin SubModelPart dirichlet + Begin SubModelPartNodes + 1 + 4 + End SubModelPartNodes +End SubModelPart diff --git a/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d_materials.json b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d_materials.json new file mode 100644 index 000000000000..ad5f9dc8be13 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/constraints/link_constraint_3d_materials.json @@ -0,0 +1,19 @@ +{ + "properties": [ + { + "model_part_name": "root", + "properties_id": 1, + "Material": { + "constitutive_law": { + "name": "LinearElastic3DLaw" + }, + "Variables": { + "DENSITY": 1.0, + "YOUNG_MODULUS": 1.0, + "POISSON_RATIO": 0.3 + }, + "Tables": null + } + } + ] +} diff --git a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py index df548b1d9217..11f1e758f5ec 100644 --- a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py +++ b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py @@ -15,6 +15,8 @@ ##### SELF-CONTAINED TESTS ##### # CL tests from test_constitutive_law import TestConstitutiveLaw as TTestConstitutiveLaw +# Constraint tests +from test_link_constraint import TestLinkConstraint # Processes test from test_mass_calculation import TestMassCalculation as TTestMassCalculation from test_compute_center_of_gravity import TestComputeCenterOfGravity as TTestComputeCenterOfGravity @@ -275,6 +277,8 @@ def AssembleTestSuites(): smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestConstitutiveLaw])) nightSuite.addTest(TInitialStateElasticityTest('test_execution')) nightSuite.addTest(TInitialStrainShellQ4ThickTest('test_execution')) + # Constraint tests + smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCase(TestLinkConstraint)) # Mass calculation tests smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestMassCalculation])) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TTestComputeCenterOfGravity])) diff --git a/applications/StructuralMechanicsApplication/tests/test_link_constraint.py b/applications/StructuralMechanicsApplication/tests/test_link_constraint.py new file mode 100644 index 000000000000..0222a7bf60f3 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/test_link_constraint.py @@ -0,0 +1,134 @@ +# --- Kratos Imports --- +import KratosMultiphysics +import KratosMultiphysics.KratosUnittest +from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis +from KratosMultiphysics.KratosUnittest import WorkFolderScope + +# --- STD Imports --- +import pathlib +import contextlib +import os +from typing import Generator +import math + + +class TestLinkConstraint(KratosMultiphysics.KratosUnittest.TestCase): + + def test_LinkConstraint2DWithoutMovedMesh(self) -> None: + self.__Run2D(False) + + def test_LinkConstraint2DWithMovedMesh(self) -> None: + self.__Run2D(True) + + def test_LinkConstraint3DWithoutMovedMesh(self) -> None: + self.__Run3D(False) + + def test_LinkConstraint3DWithMovedMesh(self) -> None: + self.__Run3D(True) + + @staticmethod + def __GetNodeDistance(first: KratosMultiphysics.Node, + second: KratosMultiphysics.Node, + initial_configuration: bool = False) -> float: + initial_positions: "list[list[float]]" = [[first.X0, first.Y0, first.Z0], + [second.X0, second.Y0, second.Z0]] + displacements: "list[list[float]]" = [[first.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X), + first.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_Y), + first.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_Z)], + [second.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X), + second.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_Y), + second.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_Z)]] + distance: float = 0.0 + for i_component in range(3): + diff: float = initial_positions[0][i_component] - initial_positions[1][i_component] + if not initial_configuration: + diff += displacements[0][i_component] - displacements[1][i_component] + distance += diff * diff + return math.sqrt(distance) + + + def __Run2D(self, move_mesh_flag: bool) -> None: + dimensions = 2 + + with WorkFolderScope("constraints", pathlib.Path(__file__).absolute()): + # Load config. + with open("link_constraint_2d.json") as project_parameters_file: + parameters = KratosMultiphysics.Parameters(project_parameters_file.read()) + + parameters["solver_settings"]["move_mesh_flag"].SetBool(move_mesh_flag) + parameters["processes"]["constraints_process_list"][1]["Parameters"]["move_mesh_flag"].SetBool(move_mesh_flag) + + # The test is set up such that 2 pseudo time steps are taken. + # + # At the end of the first one, the constraint is supposed to + # be active so distance between the constrained nodes is supposed + # to remain constant. + # + # At the end of the second step, the constraint is supposed to be + # inactive, allowing the distance between the constrained nodes to + # change. + + # Run the analysis until the first step and check whether constraints are satisfied. + parameters["problem_data"]["end_time"].SetDouble(0.5) + model = KratosMultiphysics.Model() + self.__Run(model, parameters, move_mesh_flag = move_mesh_flag) + root_model_part = model.GetModelPart("root") + + constrained_node_ids: "tuple[int,int]" = (2, 3) + first: KratosMultiphysics.Node = root_model_part.GetNode(constrained_node_ids[0]) + second: KratosMultiphysics.Node = root_model_part.GetNode(constrained_node_ids[1]) + initial_distance: float = self.__GetNodeDistance(first, second, initial_configuration = True) + distance: float = self.__GetNodeDistance(first, second) + self.assertAlmostEqual(distance, initial_distance, places = 3) + + # Run the analysis until the first step and check whether constraints are not satisfied. + parameters["problem_data"]["end_time"].SetDouble(1.0) + model = KratosMultiphysics.Model() + self.__Run(model, parameters, move_mesh_flag = move_mesh_flag) + root_model_part = model.GetModelPart("root") + + constrained_node_ids: "tuple[int,int]" = (2, 3) + first: KratosMultiphysics.Node = root_model_part.GetNode(constrained_node_ids[0]) + second: KratosMultiphysics.Node = root_model_part.GetNode(constrained_node_ids[1]) + initial_distance: float = self.__GetNodeDistance(first, second, initial_configuration = True) + distance: float = self.__GetNodeDistance(first, second) + self.assertNotAlmostEqual(distance, initial_distance, places = 0) + + + def __Run3D(self, move_mesh_flag: bool) -> None: + dimensions = 3 + + with WorkFolderScope("constraints", pathlib.Path(__file__).absolute()): + # Load config. + with open("link_constraint_3d.json") as project_parameters_file: + parameters = KratosMultiphysics.Parameters(project_parameters_file.read()) + + parameters["solver_settings"]["move_mesh_flag"].SetBool(move_mesh_flag) + parameters["processes"]["constraints_process_list"][1]["Parameters"]["move_mesh_flag"].SetBool(move_mesh_flag) + + model = KratosMultiphysics.Model() + self.__Run(model, parameters, move_mesh_flag = move_mesh_flag) + root_model_part = model.GetModelPart("root") + + constrained_node_ids: "tuple[int,int]" = (2, 3) + first: KratosMultiphysics.Node = root_model_part.GetNode(constrained_node_ids[0]) + second: KratosMultiphysics.Node = root_model_part.GetNode(constrained_node_ids[1]) + initial_distance: float = self.__GetNodeDistance(first, second, initial_configuration = True) + distance: float = self.__GetNodeDistance(first, second) + self.assertAlmostEqual(distance, initial_distance, places = 3) + + + def __Run(self, + model: KratosMultiphysics.Model, + project_parameters: KratosMultiphysics.Parameters, + move_mesh_flag: bool): + project_parameters["solver_settings"]["move_mesh_flag"].SetBool(move_mesh_flag) + project_parameters["processes"]["constraints_process_list"][1]["Parameters"]["move_mesh_flag"].SetBool(move_mesh_flag) + + KratosMultiphysics.ModelPart = model.CreateModelPart("root") + analysis_stage = StructuralMechanicsAnalysis(model, project_parameters) + analysis_stage.Run() + + +if __name__ == "__main__": + KratosMultiphysics.KratosUnittest.main() diff --git a/kratos/containers/nodal_data.h b/kratos/containers/nodal_data.h index 2f9e98a21abe..420584dd98b2 100644 --- a/kratos/containers/nodal_data.h +++ b/kratos/containers/nodal_data.h @@ -186,7 +186,7 @@ class KRATOS_API(KRATOS_CORE) NodalData final ///@{ /// Copy constructor. - NodalData(NodalData const& rOther); + NodalData(NodalData const& rOther) = default; ///@} diff --git a/kratos/includes/dof.h b/kratos/includes/dof.h index 874a2ca1837c..2fb6081ac7af 100644 --- a/kratos/includes/dof.h +++ b/kratos/includes/dof.h @@ -184,10 +184,18 @@ class Dof{ { } + Dof(Dof&&) noexcept = default; + + Dof(const Dof&) noexcept = default; + ///@} ///@name Operators ///@{ + Dof& operator=(Dof&&) noexcept = default; + + Dof& operator=(const Dof&) noexcept = default; + template typename TVariableType::Type& operator()(const TVariableType& rThisVariable, IndexType SolutionStepIndex = 0) { diff --git a/kratos/includes/kratos_application.h b/kratos/includes/kratos_application.h index aa9b30c7a7d4..4f13d4675fdc 100644 --- a/kratos/includes/kratos_application.h +++ b/kratos/includes/kratos_application.h @@ -27,6 +27,7 @@ #include "input_output/logger.h" #include "utilities/quaternion.h" #include "constraints/linear_master_slave_constraint.h" +#include "solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp" // LinearMultifreedomConstraint // Geometries definition #include "geometries/register_kratos_components_for_geometry.h" @@ -489,6 +490,7 @@ class KRATOS_API(KRATOS_CORE) KratosApplication { // Master-Slave base constraint const MasterSlaveConstraint mMasterSlaveConstraint; const LinearMasterSlaveConstraint mLinearMasterSlaveConstraint; + const LinearMultifreedomConstraint mLinearMultifreedomConstraint; // Periodic Condition const PeriodicCondition mPeriodicCondition; diff --git a/kratos/python/add_constraint_to_python.cpp b/kratos/python/add_constraint_to_python.cpp index ed25cc780302..6d1da19efa95 100644 --- a/kratos/python/add_constraint_to_python.cpp +++ b/kratos/python/add_constraint_to_python.cpp @@ -20,6 +20,7 @@ #include "includes/master_slave_constraint.h" #include "constraints/linear_master_slave_constraint.h" #include "constraints/slip_constraint.h" +#include "solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp" // LinearMultifreedomConstraint namespace Kratos::Python { @@ -133,6 +134,33 @@ void AddConstraintToPython(pybind11::module &m) >()) ; + pybind11::class_(m, "LinearMultifreedomConstraint", "Class representing (part of) a linear multifreedom constraint equation.") + .def(pybind11::init([](const IndexType Id, + const LinearMultifreedomConstraint::DofPointerVectorType& rDofs, + const std::vector& rConstraintLabels, + const LinearMultifreedomConstraint::MatrixType& rRelationMatrix, + const LinearMultifreedomConstraint::VectorType& rConstraintGaps){ + return std::make_shared( + Id, + LinearMultifreedomConstraint::DofPointerVectorType(rDofs), + rConstraintLabels, + rRelationMatrix, + rConstraintGaps); + }), + (std::string {"@brief Class representing (part of) a linear multifreedom constraint equation.\n"} + + "@param Id Identifier of the constraint instance. Unrelated to the identifier of the constrain equation.\n" + + "@param Dofs List of DoFs participating in the constrain equation(s).\n" + + "@param ConstraintLabels Identifiers of the constraint equations.\n" + + "@param RelationMatrix Matrix storing the coefficients of the participating DoFs. Each row of the matrix defines a constraint equation. The coefficients must be in the same order as the provided DoFs. The number of rows must match the size of ConstraintLabels.\n" + + "@param ConstraintGaps Constants of the constraint equations. The constraint gap vector's size must match the number of rows in RelationMatrix.").c_str(), + pybind11::arg("Id"), + pybind11::arg("Dofs"), + pybind11::arg("ConstraintLabels"), + pybind11::arg("RelationMatrix"), + pybind11::arg("ConstraintGaps")) + ; pybind11::class_ (m, "SlipConstraint") diff --git a/kratos/python/add_strategies_to_python.cpp b/kratos/python/add_strategies_to_python.cpp index 410c760f9f02..b7dcc68a1a04 100644 --- a/kratos/python/add_strategies_to_python.cpp +++ b/kratos/python/add_strategies_to_python.cpp @@ -66,6 +66,7 @@ #include "solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver.h" #include "solving_strategies/builder_and_solvers/residualbased_block_builder_and_solver.h" #include "solving_strategies/builder_and_solvers/residualbased_block_builder_and_solver_with_lagrange_multiplier.h" +#include "solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp" // Linear solvers #include "linear_solvers/linear_solver.h" @@ -512,6 +513,17 @@ namespace Kratos:: Python .def(py::init< LinearSolverType::Pointer, Parameters > ()) ; + typedef PMultigridBuilderAndSolver PMultigridBuilderAndSolverType; + py::class_(m, "PMultigridBuilderAndSolver") + .def(py::init<>()) + .def(py::init(), + (std::string {"@param LinearSolver Pointer to a linear solver. This is ignored because PMultigridBuilderAndSolver constructs its own instances internally.\n"} + + "@param Settings\n" + + "@details Default parameters:\n" + PMultigridBuilderAndSolverType().GetDefaultParameters().PrettyPrintJsonString()).c_str(), + py::arg("LinearSolver"), + py::arg("Settings")) + ; + //******************************************************************** //******************************************************************** //******************************************************************** diff --git a/kratos/solving_strategies/builder_and_solvers/builder_and_solver.cpp b/kratos/solving_strategies/builder_and_solvers/builder_and_solver.cpp new file mode 100644 index 000000000000..800f2dd2da88 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/builder_and_solver.cpp @@ -0,0 +1,21 @@ +// Project includes +#include "solving_strategies/builder_and_solvers/builder_and_solver.h" // BuilderAndSolver +#include "spaces/ublas_space.h" // TUblasSparseSpace, TUblasDenseSpace +#include "linear_solvers/linear_solver.h" // LinearSolver + + +namespace Kratos { + + +template class KRATOS_API(KRATOS_CORE) BuilderAndSolver, + TUblasDenseSpace, + LinearSolver, + TUblasDenseSpace>>; + +template class KRATOS_API(KRATOS_CORE) BuilderAndSolver, + TUblasDenseSpace, + LinearSolver, + TUblasDenseSpace>>; + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/builder_and_solver.h b/kratos/solving_strategies/builder_and_solvers/builder_and_solver.h index c5c11accb786..6de40cbbaaaf 100644 --- a/kratos/solving_strategies/builder_and_solvers/builder_and_solver.h +++ b/kratos/solving_strategies/builder_and_solvers/builder_and_solver.h @@ -60,7 +60,7 @@ template, class TLinearSolver //= LinearSolver > -class BuilderAndSolver +class KRATOS_API(KRATOS_CORE) BuilderAndSolver { public: ///@name Type Definitions diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.cpp new file mode 100644 index 000000000000..225bb3b00fe1 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.cpp @@ -0,0 +1,458 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.hpp" // AugmentedLagrangeConstraintAssembler +#include "solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp" // MapRowContribution +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_utilities.hpp" // ProcessMasterSlaveConstraint, ProcessMultifreedomConstraint +#include "solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp" // Scaling +#include "spaces/ublas_space.h" // TUblasSparseSpace, TUblasDenseSpace +#include "utilities/sparse_matrix_multiplication_utility.h" // SparseMatrixMultiplicationUtility + +// STL includes +#include // std::filesystem::current_path + + +namespace Kratos { + + +namespace detail { + + +template +void ApplyDirichletConditions(typename TSparse::MatrixType& rRelationMatrix, + typename TSparse::VectorType& rConstraintGaps, + TItDof itDofBegin, + [[maybe_unused]] TItDof itDofEnd, + const std::string& rConstraintAssemblerName, + int Verbosity) +{ + // Sanity checks. + KRATOS_ERROR_IF_NOT(static_cast(std::distance(itDofBegin, itDofEnd)) == rRelationMatrix.size2()); + KRATOS_ERROR_IF_NOT(rConstraintGaps.size() == rRelationMatrix.size1()); + + KRATOS_TRY + + /// @todo Make applying dirichlet conditions on penalty constraints more robust and efficient. + + std::size_t forced_dirichlet_dof_count, + total_forced_dirichlet_dof_count = 0ul, + iteration_count = 0ul; + std::vector mutexes(rRelationMatrix.size2()); + + do { + forced_dirichlet_dof_count = 0ul; + + IndexPartition(rRelationMatrix.size1()).for_each([&rRelationMatrix, + &rConstraintGaps, + &forced_dirichlet_dof_count, + &mutexes, + itDofBegin](std::size_t i_constraint){ + const auto i_entry_begin = rRelationMatrix.index1_data()[i_constraint]; + const auto i_entry_end = rRelationMatrix.index1_data()[i_constraint + 1]; + auto free_dof_count = i_entry_end - i_entry_begin; + + for (auto i_entry=i_entry_begin; i_entry& r_dof = *(itDofBegin + i_dof); + KRATOS_DEBUG_ERROR_IF_NOT(i_dof == r_dof.EquationId()); + + std::scoped_lock lock(mutexes[i_dof]); + if (r_dof.IsFixed()) { + --free_dof_count; + // A Dirichlet condition is set on the current DoF, which also appears in + // a constraint equation => explicitly impose this condition on the constraint + // equation. + const typename TSparse::DataType constraint_coefficient = rRelationMatrix.value_data()[i_entry]; + const typename TSparse::DataType dirichlet_condition = r_dof.GetSolutionStepValue(); + AtomicAdd(rConstraintGaps[i_constraint], + static_cast(constraint_coefficient * dirichlet_condition)); + rRelationMatrix.value_data()[i_entry] = static_cast(0); + } // if r_dof.IsFixed() + } // for i_entry in range(i_entry_begin, i_entry_end) + + if (free_dof_count == 1) { + AtomicAdd(forced_dirichlet_dof_count, static_cast(1)); + + // If there's only 1 DoF left that does not have a Dirichlet condition + // set on it, the constraint itself becomes equivalent to a Dirichlet + // condition => find and constrain the last remaining free DoF. + for (auto i_entry=i_entry_begin; i_entry& r_dof = *(itDofBegin + i_dof); + + std::scoped_lock lock(mutexes[i_dof]); + if (r_dof.IsFree()) { + const typename TSparse::DataType constraint_coefficient = rRelationMatrix.value_data()[i_entry]; + KRATOS_ERROR_IF_NOT(constraint_coefficient); + const typename TSparse::DataType constraint_gap = rConstraintGaps[i_constraint]; + const auto dirichlet_condition = -constraint_gap / constraint_coefficient; + + r_dof.GetSolutionStepValue() = dirichlet_condition; + r_dof.FixDof(); + rRelationMatrix.value_data()[i_entry] = static_cast(0); + rConstraintGaps[i_constraint] = static_cast(0); + break; + } + } // for i_entry in range(i_entry_begin, i_entry_end) + } // if free_dof_count == 1 + }); // for i_constraint in range(rRelationMatrix.size1()) + + ++iteration_count; + total_forced_dirichlet_dof_count += forced_dirichlet_dof_count; + } while (forced_dirichlet_dof_count); + + if (2 <= Verbosity && total_forced_dirichlet_dof_count) { + std::cout << rConstraintAssemblerName << ": " + << "propagated Dirichlet conditions to " << total_forced_dirichlet_dof_count << " DoFs " + << "in " << iteration_count - 1 << " iterations\n"; + } + KRATOS_CATCH("") +} + + +} // namespace detail + + +template +struct AugmentedLagrangeConstraintAssembler::Impl +{ + using Interface = AugmentedLagrangeConstraintAssembler; + + /// @brief A map associating slave IDs with constraint indices and the number of constraint objects referencing it. + std::unordered_map + > mConstraintIdMap; + + std::optional mMaybeTransposeRelationMatrix; + + std::unique_ptr mpPenaltyFunctor; + + int mVerbosity; +}; // struct AugmentedLagrangeConstraintAssembler::Impl + + +template + +AugmentedLagrangeConstraintAssembler::AugmentedLagrangeConstraintAssembler() noexcept + : AugmentedLagrangeConstraintAssembler(Parameters()) +{} + + +template +AugmentedLagrangeConstraintAssembler::AugmentedLagrangeConstraintAssembler(Parameters Settings) + : AugmentedLagrangeConstraintAssembler(Settings, "unnamed") +{ +} + + +template +AugmentedLagrangeConstraintAssembler::AugmentedLagrangeConstraintAssembler(Parameters Settings, + std::string&& rInstanceName) + : Base(ConstraintImposition::AugmentedLagrange, std::move(rInstanceName)), + mpImpl(new Impl{/*mConstraintIdToIndexMap: */std::unordered_map>(), + /*mMaybeTransposeRelationMatrix: */std::optional(), + /*mpScaling*/ nullptr, + /*mVerbosity: */1}) +{ + KRATOS_INFO_IF(this->Info(), 2 <= mpImpl->mVerbosity && !this->GetName().empty()) + << ": aliased to " << this->GetName(); + + // Parse the penalty factor and validate other settings. + Parameters default_parameters = this->GetDefaultParameters(); + Parameters default_penalty_factor = default_parameters["penalty_factor"].Clone(); + default_parameters.RemoveValue("penalty_factor"); + std::optional maybe_penalty_factor = Settings.Has("penalty_factor") ? Settings["penalty_factor"].Clone() : std::optional(); + if (maybe_penalty_factor.has_value()) Settings.RemoveValue("penalty_factor"); + Settings.ValidateAndAssignDefaults(default_parameters); + Settings.AddValue("penalty_factor", maybe_penalty_factor.has_value() ? maybe_penalty_factor.value() : default_penalty_factor); + mpImpl->mpPenaltyFunctor = std::make_unique(Settings["penalty_factor"]); + + // Parse other algorithmic settings. + Vector algorithmic_parameters(1); + algorithmic_parameters[0] = Settings["tolerance"].Get(); + this->SetValue(this->GetAlgorithmicParametersVariable(), algorithmic_parameters); + this->SetValue(NL_ITERATION_NUMBER, Settings["max_iterations"].Get()); + + // Parse miscellaneous settings. + this->mpImpl->mVerbosity = Settings["verbosity"].Get(); +} + + +template +AugmentedLagrangeConstraintAssembler::~AugmentedLagrangeConstraintAssembler() = default; + + +template +void AugmentedLagrangeConstraintAssembler::Allocate(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + typename Base::DofSet& rDofSet) +{ + KRATOS_TRY + + this->Clear(); + if (rConstraints.empty()) { + this->GetRelationMatrix() = typename TSparse::MatrixType(0, rLhs.size2(), 0); + this->GetRelationMatrix().index1_data()[0] = 0; + this->GetRelationMatrix().set_filled(1, 0); + return; + } + + detail::MakeRelationTopology(rDofSet.size(), + rConstraints, + rProcessInfo, + this->GetRelationMatrix(), + this->GetConstraintGapVector(), + mpImpl->mConstraintIdMap); + + // At this point, the sparsity pattern of the unconstrained LHS matrix (K), as well + // as the relation matrix (A) are constructed. The current policy in Kratos is that + // constraint imposition happens on K in-place, meaning K will have the union of both + // matrices' sparsity patterns (K and A^TA) after imposition. + // In the case of nonlinear constraints, the constraints' Hessian (H) also comes into + // play, which is a subset of A^TA so no new entries need to be added. + { + typename TSparse::MatrixType transpose; + SparseMatrixMultiplicationUtility::TransposeMatrix(transpose, this->GetRelationMatrix()); + SparseMatrixMultiplicationUtility::MatrixMultiplication(transpose, this->GetRelationMatrix(), this->GetHessian()); + } + MergeMatrices(rLhs, this->GetHessian()); + + KRATOS_CATCH("") +} + + +template +void AugmentedLagrangeConstraintAssembler::Assemble(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename Base::DofSet& rDofSet, + const bool AssembleLhs, + const bool AssembleRhs) +{ + KRATOS_TRY + + /// @todo @p AssembleLhs and @p AssembleRhs are currently ignored, + /// and everything gets assembled. Think it through what needs + /// to be done and what can be omitted for each case. - matekelemen + detail::AssembleRelationMatrix(rConstraints, + rProcessInfo, + this->GetRelationMatrix(), + this->GetHessian(), + this->GetConstraintGapVector(), + mpImpl->mConstraintIdMap); + + // Propagate obvious dirichlet conditions. + detail::ApplyDirichletConditions(this->GetRelationMatrix(), + this->GetConstraintGapVector(), + rDofSet.begin(), + rDofSet.end(), + this->GetName(), + mpImpl->mVerbosity); + + KRATOS_CATCH("") +} + + +template +void AugmentedLagrangeConstraintAssembler::Initialize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + typename Base::DofSet::iterator itDofBegin, + typename Base::DofSet::iterator itDofEnd) +{ + KRATOS_TRY + + // Output the relation matrix and constraint gaps if the verbosity is high enough. + if (4 <= mpImpl->mVerbosity) { + KRATOS_INFO(this->GetName() + ": ") + << "write relation matrix to " + << (std::filesystem::current_path() / "relation_matrix.mm") + << "\n"; + TSparse::WriteMatrixMarketMatrix("relation_matrix.mm", this->GetRelationMatrix(), false); + + KRATOS_INFO(this->GetName() + ": ") + << "write constraint gaps to " + << (std::filesystem::current_path() / "constraint_gaps.mm") + << "\n"; + TSparse::WriteMatrixMarketVector("constraint_gaps.mm", this->GetConstraintGapVector()); + + KRATOS_INFO(this->GetName() + ": ") + << "write constraint hessian to " + << (std::filesystem::current_path() / "hessian.mm") + << "\n"; + TSparse::WriteMatrixMarketMatrix("hessian.mm", this->GetHessian(), false); + } // if 4 <= verbosity + + // Reset the relation matrix' transpose to make sure that propagated dirichlet + // conditions are reflected in the transpose as well. + mpImpl->mMaybeTransposeRelationMatrix.reset(); + + // Fetch function-wide constants. + using SparseUtils = SparseMatrixMultiplicationUtility; + mpImpl->mpPenaltyFunctor->template Cache(rLhs); + const typename TSparse::DataType penalty_factor = this->GetPenaltyFactor(); + + // Apply initial penalty- and lagrange terms to the left- and right hand sides. + { + const typename TSparse::MatrixType& r_transpose_relation_matrix = this->GetTransposeRelationMatrix(); + KRATOS_ERROR_IF(r_transpose_relation_matrix.size1() != this->GetRelationMatrix().size2() + || r_transpose_relation_matrix.size2() != this->GetRelationMatrix().size1() + || r_transpose_relation_matrix.nnz() != this->GetRelationMatrix().nnz()) + << "the transpose of the relation matrix is uninitialized or out of date " + << "(" << this->GetRelationMatrix().size1() << "x" << this->GetRelationMatrix().size2() << "/" << this->GetRelationMatrix().nnz() << ") " + << "(" << this->GetTransposeRelationMatrix().size1() << "x" << this->GetTransposeRelationMatrix().size2() << "/" << this->GetTransposeRelationMatrix().nnz() << ") "; + + typename TSparse::MatrixType relation_product; + SparseUtils::MatrixMultiplication(r_transpose_relation_matrix, + this->GetRelationMatrix(), + relation_product); + + // Add terms to the LHS matrix. + InPlaceMatrixAdd(rLhs, + relation_product, + penalty_factor); + + InPlaceMatrixAdd(rLhs, + this->GetHessian(), + penalty_factor); + + // Add terms to the RHS vector. + BalancedProduct(r_transpose_relation_matrix, + this->GetConstraintGapVector(), + rRhs); + } + KRATOS_CATCH("") +} + + +template +void AugmentedLagrangeConstraintAssembler::InitializeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs) +{ + KRATOS_TRY + + // Compute the constraint residuals. + typename TSparse::VectorType constraint_residual = this->GetConstraintGapVector(); + BalancedProduct(this->GetRelationMatrix(), rSolution, constraint_residual); + + // Update the RHS. + BalancedProduct(this->GetTransposeRelationMatrix(), + constraint_residual, + rRhs, + static_cast(-this->GetPenaltyFactor())); + + KRATOS_CATCH("") +} + + +template +bool +AugmentedLagrangeConstraintAssembler::FinalizeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + PMGStatusStream::Report& rReport, + PMGStatusStream& rStream) +{ + KRATOS_TRY + + // Compute the constraint residuals. + typename TSparse::VectorType constraint_residual = this->GetConstraintGapVector(); + BalancedProduct(this->GetRelationMatrix(), rSolution, constraint_residual); + + // Update status. + rReport.maybe_constraint_residual = TSparse::TwoNorm(constraint_residual); + rReport.constraints_converged = rReport.maybe_constraint_residual.value() <= this->GetTolerance(); + + // Decide on whether to keep looping and log state. + if (rReport.constraints_converged) { + rStream.Submit(rReport.Tag(2), mpImpl->mVerbosity); + return true; + } /*if converged*/ else { + const int max_iterations = this->GetValue(NL_ITERATION_NUMBER); + if (static_cast(rReport.constraint_iteration) + 1 < max_iterations) { + rStream.Submit(rReport.Tag(3), mpImpl->mVerbosity); + return false; + } /*if iIteration < max_iterations*/ else { + rStream.Submit(rReport.Tag(2), mpImpl->mVerbosity); + return true; + } + } /*if not converged*/ + + KRATOS_CATCH("") +} + + +template +void AugmentedLagrangeConstraintAssembler::Finalize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + typename Base::DofSet& rDofSet) +{ +} + + +template +void AugmentedLagrangeConstraintAssembler::Clear() +{ + Base::Clear(); + mpImpl->mConstraintIdMap = decltype(mpImpl->mConstraintIdMap)(); + mpImpl->mMaybeTransposeRelationMatrix = decltype(mpImpl->mMaybeTransposeRelationMatrix)(); +} + + +template +Parameters AugmentedLagrangeConstraintAssembler::GetDefaultParameters() +{ + return Parameters(R"({ + "method" : "augmented_lagrange", + "penalty_factor" : "norm", + "tolerance" : 1e-6, + "max_iterations" : 1e1, + "verbosity" : 1 + })"); +} + + +template +typename TSparse::DataType AugmentedLagrangeConstraintAssembler::GetPenaltyFactor() const +{ + return mpImpl->mpPenaltyFunctor->Evaluate(); +} + + +template +typename TSparse::MatrixType& +AugmentedLagrangeConstraintAssembler::GetTransposeRelationMatrix() +{ + KRATOS_TRY + if (!mpImpl->mMaybeTransposeRelationMatrix.has_value()) { + mpImpl->mMaybeTransposeRelationMatrix.emplace(); + SparseMatrixMultiplicationUtility::TransposeMatrix(mpImpl->mMaybeTransposeRelationMatrix.value(), + this->GetRelationMatrix()); + } + return mpImpl->mMaybeTransposeRelationMatrix.value(); + KRATOS_CATCH("") +} + + +template class AugmentedLagrangeConstraintAssembler,TUblasDenseSpace>; + +template class AugmentedLagrangeConstraintAssembler,TUblasDenseSpace>; + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.hpp new file mode 100644 index 000000000000..c917010e69fc --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.hpp @@ -0,0 +1,139 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project Includes +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp" // ConstraintAssembler +#include "includes/kratos_parameters.h" // Parameters + + +namespace Kratos { + + +/** @brief Class implementing constraint imposition with the augmented Lagrange multiplier method. + * @details In a nutshell, the augmented Lagrange method is a combination of the penalty method + * and the method of Lagrange multipliers. The solution potentially involves multiple + * iterations, depending on the choice of penalty factor as well as the properties of + * the system to be solved. + * + * Default settings: + * @code + * { + * "method" : "augmented_lagrange", + * "penalty_factor" : "norm", + * "initial_lagrange_multiplier" : 0.0, + * "tolerance" : 1e-6, + * "max_iterations" : 1e1, + * "verbosity" : 1 + * } + * @endcode + * - @p "method" name of the constrained assembler to refer to in @ref ConstraintAssemblerFactory. + * - @p "penalty_factor" Coefficient to scale the constraint equations by. It can either be a + * numeric literal or a string representing a function of the left hand side + * matrix' diagonal. See @ref Scaling for more information. + * - @p "max" Infinity norm (abs max) of the matrix' diagonal. + * - @p "norm" 2-norm of the matrix' diagonal. + * - @p "tolerance" Target absolute norm of the constraint equations (without penalty scaling). + * - @p "max_iterations" Stop iterating if the desired tolerance was not reached after this many iterations. + * - @p "verbosity" Level of verbosity. Every level includes lower verbosity levels and + * adds new events to report: + * - @p 0 No messages, not even in case of failure. + * - @p 1 Warnings and failure messages. + * - @p 2 Aggregated status reports. + * - @p 3 Per-iteration status reports. + * - @p 4 Output system matrices and vectors. + * - @p 5 Write a VTU output with the solution and residuals at each iteration. + */ +template +class AugmentedLagrangeConstraintAssembler : public ConstraintAssembler +{ +public: + using Base = ConstraintAssembler; + + AugmentedLagrangeConstraintAssembler() noexcept; + + AugmentedLagrangeConstraintAssembler(Parameters Settings); + + AugmentedLagrangeConstraintAssembler(Parameters Settings, + std::string&& rInstanceName); + + ~AugmentedLagrangeConstraintAssembler(); + + /// @copydoc Base::Allocate + void Allocate(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + typename Base::DofSet& rDofSet) override; + + /// @copydoc Base::Assemble + void Assemble(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename Base::DofSet& rDofSet, + const bool AssembleLhs, + const bool AssembleRhs) override; + + /// @copydoc Base::Initialize + void Initialize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + typename Base::DofSet::iterator itDofBegin, + typename Base::DofSet::iterator itDofEnd) override; + + /// @copydoc Base::InitializeSolutionStep + void InitializeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs) override; + + /// @copydoc Base::FinalizeSolutionStep + bool FinalizeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + PMGStatusStream::Report& rReport, + PMGStatusStream& rStream) override; + + /// @copydoc Base::Finalize + void Finalize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + typename Base::DofSet& rDofSet) override; + + /// @copydoc Base::Clear + void Clear() override; + + static Parameters GetDefaultParameters(); + + /// @internal + typename TSparse::DataType GetPenaltyFactor() const; + + /// @internal + typename TSparse::DataType GetTolerance() const + { + return this->GetValue(AugmentedLagrangeConstraintAssembler::GetAlgorithmicParametersVariable())[0]; + } + + /// @internal + static const Variable& GetAlgorithmicParametersVariable() noexcept + { + return SHAPE_FUNCTIONS_VECTOR; + } + +private: + typename TSparse::MatrixType& GetTransposeRelationMatrix(); + + struct Impl; + std::unique_ptr mpImpl; +}; // class AugmentedLagrangeConstraintAssembler + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp new file mode 100644 index 000000000000..c96f0e98e0d3 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp @@ -0,0 +1,282 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/status_stream.hpp" // PMGStatusStream +#include "containers/data_value_container.h" // DataValueContainer +#include "containers/pointer_vector_set.h" // PointerVectorSet +#include "includes/master_slave_constraint.h" // MasterSlaveConstraint +#include "includes/indexed_object.h" // IndexedObject +#include "includes/dof.h" // Dof +#include "includes/variables.h" // IDENTIFIER +#include "includes/smart_pointers.h" // KRATOS_CLASS_POINTER_DEFINITION + + +namespace Kratos { + + +/// @brief Enum class representing imposition methods for multifreedom constraints. +enum class ConstraintImposition +{ + None = 0, + MasterSlave = 1, + Lagrange = 2, + AugmentedLagrange = 3, + Penalty = 4 +}; // enum class ConstraintImposition + + +/// @brief Interface for assembling and imposing multifreedom constraints. +template +class ConstraintAssembler : public DataValueContainer +{ +public: + KRATOS_CLASS_POINTER_DEFINITION(ConstraintAssembler); + + using DofSet = PointerVectorSet>; + + using ConstraintArray = PointerVectorSet; + + ConstraintAssembler() noexcept + : ConstraintAssembler(ConstraintImposition::MasterSlave, "unnamed") + {} + + ConstraintAssembler(ConstraintImposition Method, + std::string&& rInstanceName) + : DataValueContainer(), + mRelationMatrix(), + mConstraintGapVector(), + mName(std::move(rInstanceName)) + { + std::string method_name; + + switch (Method) { + case ConstraintImposition::None: + method_name = "none"; + break; + case ConstraintImposition::MasterSlave: + method_name = "master_slave"; + break; + case ConstraintImposition::AugmentedLagrange: + method_name = "augmented_lagrange"; + break; + default: + // Other imposition methods are not implemented yet. + KRATOS_ERROR << "Unsupported constraint imposition: " << (int)Method; + } // switch Method + + this->SetValue(ConstraintAssembler::GetImpositionVariable(), method_name); + } + + /// @details Define an overriding virtual destructor to ensure compile time errors + // if the base class' destructor turns non-virtual. + virtual ~ConstraintAssembler() override = default; + + /// @brief Allocate memory for the constraint gap vector and relation matrix, and compute its topology. + /// @details This function is responsible for large memory allocations, as well as computing + /// the sparsity pattern of the relation matrix. It must also modify the provided + /// left hand side matrix, solution vector, right hand side vector, and DoF list such that + /// these containers will not require reallocation during later stages of the solution process. + /// @param rConstraints Constraint set of the related @ref ModelPart. + /// @param rProcessInfo Current @ref ProcessInfo of the related @ref ModelPart. + /// @param rLhs Unconstrained left hand side matrix' topology. + /// @param rDofSet Unconstrained set of @ref Dof "DoFs". + /// @note This function should be invoked @b after the unconstrained system is allocated, but @b before + /// it is assembled. + virtual void Allocate(const ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + DofSet& rDofSet) + { + } + + /// @brief Compute and assemble constraint contributions into the preallocated relation matrix and constraint gap vector. + /// @details This function is responsible for computing the entries of the relation matrix + /// as well as the constraint gap vector. + /// @param rConstraints Constraint set of the related @ref ModelPart. + /// @param rProcessInfo @ref ProcessInfo of the related @ref ModelPart. + /// @param rDofSet Unconstrained set of @ref Dof "DoFs". + /// @param AssembleLhs Indicates whether to assemble data structures necessary for updating + /// the left hand side matrix of the unconstrained system. + /// @param AssembleRhs Indicates whether to assemble data structures necessary for updating + /// the right hand side vector of the unconstrained system. + /// @note This function must be preceded by a call to @ref ConstraintAssembler::Allocate, and should not make large scale + /// reallocations. + virtual void Assemble(const ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + DofSet& rDofSet, + const bool AssembleLhs, + const bool AssembleRhs) + { + } + + /// @brief Prepare the linear system for the solution loop. + /// @details This function is supposed to perform tasks on the linear system + /// that are required only once, before calls to the linear solver begin. + /// Constraint imposition methods that do not require a solution loop + /// (for example, master-slave elimination one-shots the constraints) + /// should manipulate the system here. If the set of @ref Dof "DoFs" has + /// to be changed, it should also be carried out here. + /// @param rLhs Unconstrained left hand side matrix with topology to accomodate constraint imposition. + /// @param rRhs Unconstrained right hand side vector with space to accomodate constraint imposition. + /// @param itDofBegin Iterator pointing to the first @ref Dof "DoF" of the unconstrained system. + /// @param itDofEnd Sentinel of the unconstrained system's array of @ref Dof "DoFs". + virtual void Initialize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + typename DofSet::iterator itDofBegin, + typename DofSet::iterator itDofEnd) + { + } + + /// @brief Manipulate the linear system before invoking the linear solver in the solution loop's current iteration. + /// @param rLhs Left hand side matrix. + /// @param rSolution Unconverged solution vector. + /// @param rRhs Right hand side vector. + virtual void InitializeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs) + { + } + + /// @brief Perform constraint-related tasks after invoking the linear solver in the current iteration of the solution loop. + /// @details This function is supposed to evaluate the convergence of constraint imposition, + /// decide whether to request more iterations in the solution loop. + /// @param rLhs Constrained left hand side matrix. + /// @param rSolution Converged solution vector. + /// @param rRhs Constrained right hand side vector. + /// @param rReport Status information on the solution loop. + /// @param rStream Output stream to report status to. + /// @warning The solution loop will continue indefinitely unless this function eventually + /// returns @p true. + virtual bool FinalizeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + PMGStatusStream::Report& rReport, + PMGStatusStream& rStream) = 0; + + /// @brief Perform tasks related to constraint imposition after constraints converged. + /// @param rLhs Constrained left hand side matrix. + /// @param rSolution Converged solution vector. + /// @param rRhs Constrained right hand side vector. + /// @param rDofSet Constrained set of @ref Dof "DoFs". + virtual void Finalize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + DofSet& rDofSet) + { + this->Apply(rSolution); + } + + /// @brief Transform the input vector from the independent space to the constrained space. + /// @param rSolution Input solution vector in the independent space. + virtual void Apply(typename TSparse::VectorType& rSolution) const + { + } + + /// @brief Release memory related to the linear system and constraints. + /// @details Derived classes must call the @ref ConstraintAssembler::Clear "Clear" + /// function of their parents at some point. + virtual void Clear() + { + mRelationMatrix = typename TSparse::MatrixType(); + mConstraintGapVector = typename TSparse::VectorType(); + } + + ConstraintImposition GetImposition() const + { + const std::string method_name = this->GetValue(ConstraintAssembler::GetImpositionVariable()); + if (method_name == "none") { + return ConstraintImposition::None; + } else if (method_name == "master_slave") { + return ConstraintImposition::MasterSlave; + } else if (method_name == "augmented_lagrange") { + return ConstraintImposition::AugmentedLagrange; + } else { + // Other imposition methods are not implemented yet. + KRATOS_ERROR << "Unsupported constraint imposition: \"" << method_name << "\""; + } + } + + const typename TSparse::MatrixType& GetRelationMatrix() const noexcept + { + return mRelationMatrix; + } + + const typename TSparse::MatrixType& GetHessian() const noexcept + { + return mHessian; + } + + const typename TSparse::VectorType& GetConstraintGapVector() const noexcept + { + return mConstraintGapVector; + } + + const std::string& GetName() const noexcept + { + return mName; + } + +protected: + // PGrid inherits an assembled unconstrained system, as well as an assembled relation matrix, + // and constructs a restriction operator it can use to directly compute the system and relation + // matrix at its own level. As a result, it does not need to do allocation and assembly, which + // means it should be able to directly define the relation matrix and constraint gap vector + // of its own constraint assembler. + // => long story short, PGrid needs access to the protected members of this class. + template + friend class PGrid; + + typename TSparse::MatrixType& GetRelationMatrix() noexcept + { + return mRelationMatrix; + } + + typename TSparse::MatrixType& GetHessian() noexcept + { + return mHessian; + } + + typename TSparse::VectorType& GetConstraintGapVector() noexcept + { + return mConstraintGapVector; + } + +private: + ConstraintAssembler(ConstraintAssembler&&) noexcept = default; + + ConstraintAssembler(const ConstraintAssembler&) = delete; + + ConstraintAssembler& operator=(const ConstraintAssembler&) = delete; + + ConstraintAssembler& operator=(ConstraintAssembler&&) noexcept = default; + + static const Variable& GetImpositionVariable() noexcept + { + return IDENTIFIER; + } + + typename TSparse::MatrixType mRelationMatrix; + + typename TSparse::MatrixType mHessian; + + typename TSparse::VectorType mConstraintGapVector; + + std::string mName; +}; // class ConstraintImposition + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.cpp new file mode 100644 index 000000000000..e89b0d39e327 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.cpp @@ -0,0 +1,59 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.hpp" // ConstraintAssemblerFactory +#include "solving_strategies/builder_and_solvers/p_multigrid/noop_constraint_assembler.hpp" // NoOpConstraintAssembler +#include "solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.hpp" // MasterSlaveConstraintAssembler +#include "solving_strategies/builder_and_solvers/p_multigrid/augmented_lagrange_constraint_assembler.hpp" // AugmentedLagrangeConstraintAssembler +#include "spaces/ublas_space.h" // TUblasSparseSpace, TUblasDenseSpace + + +namespace Kratos { + + +template +typename ConstraintAssembler::Pointer +ConstraintAssemblerFactory(Parameters Settings, + std::string&& rInstanceName) +{ + KRATOS_TRY + const std::string imposition_name = Settings["method"].Get(); + if (imposition_name == "none") { + return std::make_shared>(Settings, std::move(rInstanceName)); + } else if (imposition_name == "master_slave") { + return std::make_shared>(Settings, std::move(rInstanceName)); + } else if (imposition_name == "augmented_lagrange") { + return std::make_shared>(Settings, std::move(rInstanceName)); + } else { + std::stringstream message; + message << "Unsupported constraint imposition \"" << imposition_name << "\". Options are:\n" + << "\t\"none\"\n" + << "\t\"master_slave\"\n" + << "\t\"augmented_lagrange\""; + KRATOS_ERROR << message.str(); + } + KRATOS_CATCH("") +} + + +#define KRATOS_INSTANTIATE_CONSTRAINT_ASSEMBLER_FACTORY(TSparse, TDense) \ + template ConstraintAssembler::Pointer ConstraintAssemblerFactory(Parameters, std::string&&) + +KRATOS_INSTANTIATE_CONSTRAINT_ASSEMBLER_FACTORY(TUblasSparseSpace, TUblasDenseSpace); + +KRATOS_INSTANTIATE_CONSTRAINT_ASSEMBLER_FACTORY(TUblasSparseSpace, TUblasDenseSpace); + +#undef KRATOS_INSTANTIATE_CONSTRAINT_ASSEMBLER_FACTORY + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.hpp new file mode 100644 index 000000000000..d110899cb344 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.hpp @@ -0,0 +1,41 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project Includes +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp" // ConstraintAssembler +#include "includes/kratos_parameters.h" // Parameters + + +namespace Kratos { + + +/** @brief Construct a @ref ConstraintAssembler instance. + * @param Settings Parameters configuring the requested constraint assembler. + * @param rInstanceName Name of the constraint assembler instance to be constructed. + * @details This factory operates similary to @ref LinearSolverFactory. + * @p "method" setting controls what type of constraint assembler to construct, + * while the rest of the settings are passed on to the instance's constructor. + * + * Options for @p "method": + * - @p "master_slave" @ref MasterSlaveConstraintAssembler + * - @p "augmented_lagrange" @ref AugmentedLagrangeConstraintAssembler + * - @p "none" @ref NoOpConstraintAssembler + */ +template +typename ConstraintAssembler::Pointer +ConstraintAssemblerFactory(Parameters Settings, + std::string&& rInstanceName); + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_utilities.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_utilities.hpp new file mode 100644 index 000000000000..1b2fadd68749 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/constraint_utilities.hpp @@ -0,0 +1,398 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp" // ConstraintAssembler +#include "solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp" // MakeSparseTopology +#include "includes/master_slave_constraint.h" // MasterSlaveConstraint +#include "includes/variables.h" // CONSTRAINT_LABELS + +// STL includes +#include // std::vector +#include // std::unordered_map + + +namespace Kratos::detail { + + +/// @brief Compute the constraint gradient and gaps of a @ref LinearMasterSlaveConstraint. +/// @details This utility function computes the constraint gradient of constraint gaps of a +/// constraint that does not inherit from @ref MultifreedomConstraint. It also collects +/// the constraint indices of the returned entries, as well as the IDs of the @ref Dof "DoFs" +/// involved. +/// @param rConstraintIndices Output vector containing the indices of each returned constraint equation. +/// @param rDofIds Output vector of @ref Dof IDs that participate in the returned constraint equations. +/// @param rRelationMatrix Output matrix containing the gradients of each constraint equation defined by +/// the input constraint. Each row represents a gradient of the constraint equation +/// whose constraint index is defined by the corresponding entry in @p rConstraintIndices. +/// Columns represent coefficients of DoFs that are identified by their IDs at the +/// corresponding position in @p rDofIds. +/// @param rConstraintGaps Output vector of constraint gaps related to constraint equations defined by +/// @p rConstraintIndices. +/// @param rConstraint Constraint to extract values from. It must not be an instance inheriting from +/// @ref MultifreedomConstraint. +/// @param rSlaveDofIds List of slave @ref Dof "DoFs"' IDs. +/// @param rMasterDofIds List of master @ref Dof "DoFs"' IDs. +/// @param rConstraintIdMap Map relating slave DoF IDs to constraint equation indices, as well as the number +/// of constraint objects defining the constraint equation. +/// @see ProcessMultifreedomConstraint +inline void ProcessMasterSlaveConstraint(std::vector& rConstraintIndices, + std::vector& rDofIds, + MasterSlaveConstraint::MatrixType& rRelationMatrix, + [[maybe_unused]] MasterSlaveConstraint::VectorType& rConstraintGaps, + const MasterSlaveConstraint& rConstraint, + const std::vector& rSlaveDofIds, + const std::vector& rMasterDofIds, + const std::unordered_map>& rConstraintIdMap) +{ + // Constraint identifiers are the slave DoFs' IDs. + rConstraintIndices.resize(rSlaveDofIds.size()); + std::transform(rSlaveDofIds.begin(), + rSlaveDofIds.end(), + rConstraintIndices.begin(), + [&rConstraintIdMap](std::size_t slave_id){ + return rConstraintIdMap.at(slave_id).first; + }); + + // Comb DoFs. + rDofIds.resize(rMasterDofIds.size() + rSlaveDofIds.size()); + std::copy(rMasterDofIds.begin(), rMasterDofIds.end(), rDofIds.begin()); + std::copy(rSlaveDofIds.begin(), rSlaveDofIds.end(), rDofIds.begin() + rMasterDofIds.size()); + + // Reconstruct the constraint equations. + if (rRelationMatrix.size2()) { + rRelationMatrix.resize(rRelationMatrix.size1(), + rRelationMatrix.size2() + rSlaveDofIds.size(), + /*preserve entries=*/true); + for (std::size_t i_slave=0ul; i_slave + // 0 = A [u^m \\ u^s] + b +} + + +/// @brief Collect @ref Dof IDs and constraint equation indices of a @ref MultifreedomConstraint. +/// @param rConstraintIndices Output vector containing the indices of each returned constraint equation. +/// @param rDofIds Output vector of @ref Dof IDs that participate in the returned constraint equations. +/// @param rConstraint Constraint to extract values from. It must not be an instance inheriting from +/// @ref MultifreedomConstraint. +/// @param rMasterDofIds List of master @ref Dof "DoFs"' IDs. +/// @param rConstraintIdMap Map relating slave DoF IDs to constraint equation indices, as well as the number +/// of constraint objects defining the constraint equation. +/// @see ProcessMasterSlaveConstraint +inline void ProcessMultifreedomConstraint(std::vector& rConstraintIndices, + std::vector& rDofIds, + const MasterSlaveConstraint& rConstraint, + const std::vector& rMasterDofIds, + const std::unordered_map>& rConstraintIdMap) +{ + const auto& r_constraint_labels = rConstraint.GetData().GetValue(CONSTRAINT_LABELS); + rConstraintIndices.resize(r_constraint_labels.size()); + std::transform(r_constraint_labels.begin(), + r_constraint_labels.end(), + rConstraintIndices.begin(), + [&rConstraintIdMap](std::size_t constraint_label){ + return rConstraintIdMap.at(constraint_label).first; + }); + + rDofIds = rMasterDofIds; +} + + +template +void MakeRelationTopology(std::size_t SystemSize, + const typename ConstraintAssembler::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename TSparse::MatrixType& rRelationMatrix, + typename TSparse::VectorType& rConstraintGaps, + std::unordered_map>& rConstraintIdMap) +{ + KRATOS_TRY + + // Build constraint ID => index map. + rConstraintIdMap.clear(); + + { + MasterSlaveConstraint::IndexType i_constraint = 0; + MasterSlaveConstraint::EquationIdVectorType constraint_labels, master_ids; + + for (const auto& r_constraint : rConstraints) { + r_constraint.EquationIdVector(constraint_labels, master_ids, rProcessInfo); + + if (constraint_labels.empty()) { + // Constraint is a MultifreedomConstraint. + const auto& r_constraint_labels = r_constraint.GetData().GetValue(CONSTRAINT_LABELS); + constraint_labels.resize(r_constraint_labels.size()); + std::copy(r_constraint_labels.begin(), + r_constraint_labels.end(), + constraint_labels.begin()); + } /*if constraint_labels.empty()*/ + + for (const auto constraint_label : constraint_labels) { + const auto emplace_result = rConstraintIdMap.emplace(static_cast(constraint_label), + std::make_pair(i_constraint, 0ul)); + ++emplace_result.first->second.second; //< increment the number of constraint objects defining the constraint equation. + if (emplace_result.second) ++i_constraint; + } // for constraint_label in constraint_labels + } // for r_constraint in rModelPart.MasterSlaveConstraints + } + + { + std::vector> indices(rConstraintIdMap.size()); + std::vector mutexes(rConstraintIdMap.size()); + + struct TLS { + MasterSlaveConstraint::EquationIdVectorType slaves, masters; + std::vector constraint_labels; + }; + + block_for_each(rConstraints, + TLS(), + [&mutexes, &indices, &rProcessInfo, &rConstraintIdMap](const auto& r_constraint, TLS& r_tls) { + r_constraint.EquationIdVector(r_tls.slaves, r_tls.masters, rProcessInfo); + + if (r_tls.slaves.empty()) { + // Constraint is a MultifreedomConstraint. + const auto& r_constraint_labels = r_constraint.GetData().GetValue(CONSTRAINT_LABELS); + r_tls.constraint_labels.resize(r_constraint_labels.size()); + std::copy(r_constraint_labels.begin(), + r_constraint_labels.end(), + r_tls.constraint_labels.begin()); + } /*if r_tls.slaves.empty()*/ else { + // Constraint is a MasterSlaveConstraint. + r_tls.constraint_labels = r_tls.slaves; + r_tls.masters.insert(r_tls.masters.end(), + r_tls.slaves.begin(), + r_tls.slaves.end()); + } + + for (const auto i_slave : r_tls.constraint_labels) { + const auto i_constraint = rConstraintIdMap[static_cast(i_slave)].first; + std::scoped_lock lock(mutexes[i_constraint]); + indices[i_constraint].insert(r_tls.masters.begin(), r_tls.masters.end()); + } // for i_slave in slave_ids + }); // for r_constraint in rModelPart.MasterSlaveConstraints() + + MakeSparseTopology(indices, + SystemSize, + rRelationMatrix, + /*EnsureDiagonal=*/false); + TSparse::SetToZero(rRelationMatrix); + rConstraintGaps.resize(rRelationMatrix.size1(), false); + //TSparse::SetToZero(rConstraintGaps); //< unnecessary + } + + KRATOS_CATCH("") +} + + +template +void AssembleRelationMatrix(const typename ConstraintAssembler::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename TSparse::MatrixType& rRelationMatrix, + typename TSparse::MatrixType& rHessian, + typename TSparse::VectorType& rConstraintGaps, + std::unordered_map>& rConstraintIdMap) +{ + KRATOS_TRY + + // Function-wide variables. + std::vector constraint_mutexes(rConstraintIdMap.size()); + std::vector dof_mutexes(rHessian.size1()); + + // Init. + TSparse::SetToZero(rRelationMatrix); + TSparse::SetToZero(rHessian); + TSparse::SetToZero(rConstraintGaps); + + struct TLS { + MasterSlaveConstraint::EquationIdVectorType slave_ids, master_ids; + std::vector constraint_indices, dof_equation_ids, dof_index_array, reordered_dof_equation_ids; + std::vector matrix_row; + MasterSlaveConstraint::MatrixType relation_matrix, hessian; + MasterSlaveConstraint::VectorType constraint_gaps; + }; + + // Constraint assembly. + block_for_each(rConstraints, + TLS(), + [&constraint_mutexes, + &dof_mutexes, + &rProcessInfo, + &rConstraintIdMap, + &rRelationMatrix, + &rHessian, + &rConstraintGaps](const MasterSlaveConstraint& r_constraint, + TLS& r_tls){ + if (r_constraint.IsActive()) { + r_constraint.EquationIdVector(r_tls.slave_ids, + r_tls.master_ids, + rProcessInfo); + r_constraint.CalculateLocalSystem(r_tls.relation_matrix, + r_tls.constraint_gaps, + rProcessInfo); + + // MasterSlaveConstraints and MultifreedomConstraints have to be handled separately. + // MasterSlaveConstraint has both slave and master DoFs, and interprets the relation + // matrix as a map from master DoFs to slave DoFs. The constraint equation it belongs + // to is defined by the ID of the slave DoFs (MasterSlaveConstraints sharing slave DoFs + // belong to the same constraint equation), + // On the other hand, MultifreedomConstraint only fills the master DoFs, and provides a + // relation matrix whose rows directly represent coefficients in the constraint equation. + // Constraint equations are identified by the CONSTRAINT_LABELS variable stored in the + // DataValueConstainer of the constraint object (constraints with identical CONSTRAINT_LABELS + // belong to the same constraint equation). + if (r_tls.slave_ids.empty()) { + // The constraint is a MultifreedomConstraint. + detail::ProcessMultifreedomConstraint(r_tls.constraint_indices, + r_tls.dof_equation_ids, + r_constraint, + r_tls.master_ids, + rConstraintIdMap); + const auto& r_hessian = r_constraint.GetData().GetValue(GEOMETRIC_STIFFNESS_MATRIX); + r_tls.hessian.resize(r_hessian.size1(), r_hessian.size2()); + std::copy(r_hessian.data().begin(), + r_hessian.data().end(), + r_tls.hessian.data().begin()); + } /*if r_tls.slave_ids.empty()*/ else { + // The constraint is a MasterSlaveConstraint. + detail::ProcessMasterSlaveConstraint(r_tls.constraint_indices, + r_tls.dof_equation_ids, + r_tls.relation_matrix, + r_tls.constraint_gaps, + r_constraint, + r_tls.slave_ids, + r_tls.master_ids, + rConstraintIdMap); + + // The standard MasterSlaveConstraint can only represent linear constraints, + // whose Hessian vanish, so there's no need assemble them into the global + // constraint Hessian. + r_tls.hessian.resize(0, 0); + } // not r_tls.slave_ids.empty() + + KRATOS_ERROR_IF_NOT(r_tls.constraint_indices.size() == r_tls.relation_matrix.size1()) + << "constraint " << r_constraint.Id() << " is ill-formed: " + << "has " << r_tls.constraint_indices.size() << " constraint equations, but its relation matrix is " + << r_tls.relation_matrix.size1() << "x" << r_tls.relation_matrix.size2(); + + KRATOS_ERROR_IF_NOT(r_tls.dof_equation_ids.size() == r_tls.relation_matrix.size2()) + << "constraint " << r_constraint.Id() << " is ill-formed: " + << "defined on " << r_tls.dof_equation_ids.size() << " DoFs, but its relation matrix is " + << r_tls.relation_matrix.size1() << "x" << r_tls.relation_matrix.size2(); + + KRATOS_ERROR_IF_NOT(r_tls.constraint_gaps.size() == r_tls.relation_matrix.size1()) + << "constraint " << r_constraint.Id() << " is ill formed: " + << "relation matrix is " << r_tls.relation_matrix.size1() << "x" << r_tls.relation_matrix.size2() + << " but the constraint gap vector is of size " << r_tls.constraint_gaps.size(); + + r_tls.matrix_row.resize(r_tls.relation_matrix.size2()); + + // Sort DoFs based on their equation IDs because the CSR format expects rows to be sorted. + r_tls.dof_index_array.resize(r_tls.dof_equation_ids.size()); + std::iota(r_tls.dof_index_array.begin(), r_tls.dof_index_array.end(), 0ul); + std::sort(r_tls.dof_index_array.begin(), + r_tls.dof_index_array.end(), + [&r_tls](const std::size_t i_left, const std::size_t i_right) -> bool { + return r_tls.dof_equation_ids[i_left] < r_tls.dof_equation_ids[i_right]; + }); + + r_tls.reordered_dof_equation_ids.resize(r_tls.dof_equation_ids.size()); + std::transform(r_tls.dof_index_array.begin(), + r_tls.dof_index_array.end(), + r_tls.reordered_dof_equation_ids.begin(), + [&r_tls](const auto i_column){ + return r_tls.dof_equation_ids[i_column]; + }); + + // Assemble local rows into the global relation matrix. + for (std::size_t i_row=0ul; i_row lock(constraint_mutexes[i_constraint]); + MapRowContribution(rRelationMatrix, + r_tls.relation_matrix, + i_constraint, + i_row, + r_tls.reordered_dof_equation_ids); + } + + AtomicAdd(rConstraintGaps[i_constraint], + static_cast(r_tls.constraint_gaps[i_row])); + } // for i_row in range(r_tls.slave_ids.size) + + // Map contributions to the Hessian. + if (r_tls.hessian.size1() && r_tls.hessian.size2()) { + for (std::size_t i_row=0ul; i_row lock(dof_mutexes[i_row_dof]); + MapRowContribution(rHessian, + r_tls.hessian, + i_row_dof, + i_reordered_row, + r_tls.reordered_dof_equation_ids); + } // for i_row in range(r_tls.hessian.size1()) + } // if the hessian does not vanish + } // if r_constraint.IsActive + }); // for r_constraint in rModelPart.MasterSlaveConstraints + + KRATOS_CATCH("") +} + + +} // namespace Kratos::detail diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.cpp new file mode 100644 index 000000000000..b125993694d0 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.cpp @@ -0,0 +1,174 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// External includes +#include "tinyexpr/tinyexpr.h" // te_expr, te_compile, te_eval, te_free + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp" // DiagonalScaling, ParseDiagonalScaling, GetDiagonalScaleFactor +#include "spaces/ublas_space.h" // TUblasSparseSpace +#include "utilities/reduction_utilities.h" // MaxReduction, AbsMaxReduction +#include "utilities/profiler.h" // KRATOS_PROFILE_SCOPE + +// STL includes +#include // std::optional + + +namespace Kratos { + + +struct Scaling::Impl { + using ExpressionPtr = std::unique_ptr>; + + std::optional mMaybeAbsMax; + + std::optional mMaybeNorm; + + mutable ExpressionPtr mpExpression; +}; // struct Scaling::Impl + + +Scaling::Scaling() + : Scaling(Parameters()) +{ +} + + +Scaling::Scaling(Parameters Settings) + : mpImpl(new Impl {/*mMaybeAbsMax=*/0.0, + /*mMaybeNorm=*/ 0.0, + /*mpExpression*/ Impl::ExpressionPtr(nullptr, [](te_expr*){})}) +{ + // Sanity checks. + KRATOS_ERROR_IF_NOT(Settings.IsNumber() || Settings.Is()) + << "Expecting a numeric literal or a string, but got " << Settings; + + KRATOS_TRY + + std::string expression; + if (Settings.IsNumber()){ + expression = std::to_string(Settings.Get()); + } /*if Settings.IsNumber*/ else { + expression = Settings.Get(); + } /*if Settings.Is*/ + + te_variable variables[] = { + {"max", &mpImpl->mMaybeAbsMax.value()}, + {"norm", &mpImpl->mMaybeNorm.value()} + }; + + int error = 0; + mpImpl->mpExpression = Impl::ExpressionPtr(/*function ptr: */te_compile(expression.c_str(), variables, 2, &error), + /*deleter: */ [](te_expr* p){if (p) te_free(p);}); + KRATOS_ERROR_IF_NOT(mpImpl->mpExpression.get()) + << "Failed to compile expression (error at position " << error - 1 << " in \"" << expression << "\"). " + << "Make sure the expression is compatible with TinyExpr and is a function of " + << "only \"max\" and \"norm\"."; + + this->ClearCache(); + KRATOS_CATCH("") +} + + +Scaling::~Scaling() = default; + + +void Scaling::ClearCache() +{ + mpImpl->mMaybeAbsMax.reset(); + mpImpl->mMaybeNorm.reset(); +} + + +double Scaling::Evaluate() const +{ + KRATOS_ERROR_IF_NOT(mpImpl->mMaybeAbsMax.has_value() && mpImpl->mMaybeNorm.has_value()) + << "Matrix diagonal norms are not cached. Make sure to invoke Scaling::Cache before Scaling::Evaluate."; + + return te_eval(mpImpl->mpExpression.get()); +} + + +void Scaling::Cache(double DiagonalMax, double DiagonalNorm) +{ + mpImpl->mMaybeAbsMax = DiagonalMax; + mpImpl->mMaybeNorm = DiagonalNorm; +} + + +template +void NormalizeRows(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs) +{ + KRATOS_TRY + IndexPartition(rLhs.size1()).for_each([&rLhs, &rRhs](const std::size_t i_row){ + const auto maybe_diagonal_entry = detail::FindDiagonal(rLhs, i_row); + KRATOS_ERROR_IF_NOT(maybe_diagonal_entry.has_value()) << "row " << i_row << " has no diagonal entry"; + const auto diagonal_entry = maybe_diagonal_entry.value(); + + if (diagonal_entry) { + const auto scale = static_cast(1) / diagonal_entry; + + const auto i_entry_begin = rLhs.index1_data()[i_row]; + const auto i_entry_end = rLhs.index1_data()[i_row + 1]; + for (auto i_entry=i_entry_begin; i_entry(1); + } // for i_entry in range(i_entry_begin, i_entry_end) + //KRATOS_ERROR_IF(rRhs[i_row]) << "row " << i_row << " is empty, but the corresponding RHS component is " << rRhs[i_row]; + } // if not diagonal_entry + }); + KRATOS_CATCH("") +} + + +template +void NormalizeSystem(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + typename TSparse::DataType Coefficient) +{ + block_for_each(&*rLhs.value_data().begin(), + (&*rLhs.value_data().begin()) + rLhs.value_data().size(), + [Coefficient](typename TSparse::DataType& r_entry){r_entry *= Coefficient;}); + block_for_each(&*rRhs.begin(), + (&*rRhs.begin()) + rRhs.size(), + [Coefficient](typename TSparse::DataType& r_entry){r_entry *= Coefficient;}); +} + + +// Explicit instantiations for UBLAS. +#define KRATOS_DEFINE_DIAGONAL_SCALE_FUNCTIONS(TSparseSpace) \ + template void NormalizeRows(typename TSparseSpace::MatrixType&, \ + typename TSparseSpace::VectorType&); \ + template void NormalizeSystem(typename TSparseSpace::MatrixType&, \ + typename TSparseSpace::VectorType&, \ + typename TSparseSpace::DataType) + + +KRATOS_DEFINE_DIAGONAL_SCALE_FUNCTIONS(TUblasSparseSpace); + +KRATOS_DEFINE_DIAGONAL_SCALE_FUNCTIONS(TUblasSparseSpace); + +#undef KRATOS_DEFINE_DIAGONAL_SCALE_FUNCTIONS + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp new file mode 100644 index 000000000000..4738a663ada3 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp @@ -0,0 +1,149 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "includes/kratos_parameters.h" // Parameters +#include "utilities/reduction_utilities.h" // CombinedReduction, MaxReduction, SumReduction + +// STL includes +#include // std::unique_ptr +#include // std::optional + + +namespace Kratos { + + +/** @brief Class handling expressions involving the diagonal of a sparse matrix. + * @details + * @subsection Configuration + * This class can either be configured with a numeric literal that does not + * depend on the matrix's diagonal, or a string that is a function of: + * - @p "max" The infinitiy norm (abs max) of the matrix' diagonal. + * - @p "norm" 2-norm of the matrix' diagonal. + * + * @subsection Usage + * Invoke @ref Scaling::Cache "Cache" with a matrix to evaluate the stored + * expression for later use. The computed value can be retrived from + * @ref Scaling::Evaluate "Evaluate". You should call + * @ref Scaling::ClearCache "ClearCache" whenever the input matrix changes. + * + * @note The input expression is parsed and evaluated with @p tinyexpr. + */ +class KRATOS_API(KRATOS_CORE) Scaling final +{ +public: + /// @brief Construct using the default settings. + Scaling(); + + /// @brief Construct from @ref Parameters representing a numeric literal or a string. + /// @details + Scaling(Parameters Settings); + + /// @internal + ~Scaling(); + + /// @brief Compute and store the provided expression using the diagonal of an input matrix. + /// @note Invoking this function is mandatory before calling @ref Scaling::Evaluate "Evaluate". + template + void Cache(const typename TSparse::MatrixType& rMatrix); + + /// @brief Reset the expression's cached value. + /// @note Subsequent calls to @ref Scaling::Evaluate will fail unless @ref Scaling::Cache is called beforehand. + void ClearCache(); + + /// @brief Get the cached value of the stored expression. + /// @warning A call to @ref Scaling::Cache must precede any call to @p Evaluate. + double Evaluate() const; + +private: + void Cache(double DiagonalMax, double DiagonalNorm); + + Scaling(Scaling&&) = delete; + + Scaling(const Scaling&) = delete; + + Scaling& operator=(Scaling&&) = delete; + + Scaling& operator=(const Scaling&) = delete; + + struct Impl; + std::unique_ptr mpImpl; +}; // class Scaling + + +/// @internal +template +void NormalizeRows(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs); + + +/// @internal +template +void NormalizeSystem(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + typename TSparse::DataType Coefficient); + + +namespace detail { + + +/// @internal +template +std::optional FindDiagonal(const typename TSparse::MatrixType& rMatrix, + const std::size_t iRow) noexcept +{ + const auto i_entry_begin = rMatrix.index1_data()[iRow]; + const auto i_entry_end = rMatrix.index1_data()[iRow + 1]; + + const auto it_column_begin = rMatrix.index2_data().begin() + i_entry_begin; + const auto it_column_end = rMatrix.index2_data().begin() + i_entry_end; + + // Look for the diagonal entry in the current row. + const auto it_column = std::lower_bound(it_column_begin, it_column_end, iRow); + + return (it_column == it_column_end || *it_column != iRow) + ? std::optional() + : rMatrix.value_data()[rMatrix.index1_data()[iRow] + std::distance(it_column_begin, it_column)]; +} + + +} // namespace detail + + +/// @internal +template +void Scaling::Cache(const typename TSparse::MatrixType& rMatrix) +{ + KRATOS_TRY + + using Value = typename TSparse::DataType; + using Reduction = CombinedReduction,SumReduction>; + + Value norm_coefficient = rMatrix.size1() + ? std::pow(static_cast(1) / static_cast(rMatrix.size1()), static_cast(2)) + : static_cast(1); + + const auto [abs_max, square_norm] = IndexPartition(rMatrix.size1()).template for_each([&rMatrix, norm_coefficient](auto i_row){ + const auto maybe_diagonal_entry = detail::FindDiagonal(rMatrix, i_row); + const Value diagonal_entry = maybe_diagonal_entry.has_value() ? *maybe_diagonal_entry : static_cast(0); + return std::make_tuple(std::abs(diagonal_entry), norm_coefficient * diagonal_entry * diagonal_entry); + }); + + this->Cache(abs_max, std::sqrt(square_norm)); + + KRATOS_CATCH("") +} + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.cpp new file mode 100644 index 000000000000..5ed7ce5014af --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.cpp @@ -0,0 +1,58 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp" // MultifreedomConstraint +#include "includes/variables.h" // GEOMETRIC_STIFFNESS_MATRIX, INTERNAL_FORCES_VECTOR + + +namespace Kratos { + + +LinearMultifreedomConstraint::LinearMultifreedomConstraint(const IndexType Id, + DofPointerVectorType&& rDofs, + const std::vector& rConstraintLabels, + const MatrixType& rRelationMatrix, + const VectorType& rConstraintGaps) + : MultifreedomConstraint(Id, std::move(rDofs), rConstraintLabels) +{ + // Sanity checks. + KRATOS_ERROR_IF_NOT(rRelationMatrix.size1() == rConstraintGaps.size()) + << "relation matrix (" << rRelationMatrix.size1() << "x" << rRelationMatrix.size2() << ") " + << "is incompatible with its constraint gap vector (" << rConstraintGaps.size() << ")"; + + KRATOS_ERROR_IF_NOT(this->GetDofs().size() == rRelationMatrix.size2()) + << "provided number of DoFs (" << this->GetDofs().size() << ") " + << "is incompatible with the relation matrix " + << "(" << rRelationMatrix.size1() << "x" << rRelationMatrix.size2() << ")"; + + KRATOS_ERROR_IF_NOT(rConstraintLabels.size() == rRelationMatrix.size1()) + << "constraint label vector (" << rConstraintLabels.size() << ") " + << "in incompatible with the relation matrix " + << "(" << rRelationMatrix.size1() << "x" << rRelationMatrix.size2() << ")"; + + // Store input data. + this->SetValue(GEOMETRIC_STIFFNESS_MATRIX, rRelationMatrix); + this->SetValue(INTERNAL_FORCES_VECTOR, rConstraintGaps); +} + + +void LinearMultifreedomConstraint::CalculateLocalSystem(MatrixType& rRelationMatrix, + VectorType& rConstraintGaps, + const ProcessInfo&) const +{ + rRelationMatrix = this->GetData().GetValue(GEOMETRIC_STIFFNESS_MATRIX); + rConstraintGaps = this->GetData().GetValue(INTERNAL_FORCES_VECTOR); +} + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp new file mode 100644 index 000000000000..3921aa74925f --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp @@ -0,0 +1,103 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.hpp" // MultifreedomConstraint +#include "includes/smart_pointers.h" // KRATOS_CLASS_POINTER_DEFINITION +#include "includes/code_location.h" // KRATOS_CODE_LOCATION + +// System includes +#include // std::numeric_limits + + +namespace Kratos { + + +/** @brief Class representing (part of) a linear constraint equation. + * @details The relation matrix is stored in the GEOMETRIC_STIFFNESS_MATRIX + * while the constraint gap in the INTERNAL_FORCES_VECTOR variable + * of the instance's @ref DataValueContainer. + */ +class KRATOS_API(KRATOS_CORE) LinearMultifreedomConstraint final + : public MultifreedomConstraint +{ +public: + KRATOS_CLASS_POINTER_DEFINITION(LinearMultifreedomConstraint); + + using MultifreedomConstraint::IndexType; + + using MultifreedomConstraint::DofPointerVectorType; + + using MultifreedomConstraint::MatrixType; + + using MultifreedomConstraint::VectorType; + + LinearMultifreedomConstraint() noexcept + : LinearMultifreedomConstraint(std::numeric_limits::max()) + {} + + /// @copydoc MultifreedomConstraint::MultifreedomConstraint(const IndexType) + LinearMultifreedomConstraint(const IndexType Id) noexcept + : LinearMultifreedomConstraint(Id, + DofPointerVectorType {}, + std::vector {}, + MatrixType {}, + VectorType {}) + {} + + /// @brief Construct a constraint instance will all necessary information. + /// @param Id Identifier of the constraint instance. Unrelated to the identifier of the constrain equation. + /// @param rDofs DoFs participating in the constrain equation(s). + /// @param rConstraintLabels Identifiers of the constraint equations. + /// @param rRelationMatrix Matrix storing the coefficients of the participating DoFs. Each row of the matrix + /// defines a constraint equation. The coefficients must be in the same order as the + /// DoFs in @p rDofs. The number of rows must match the size of @p rConstraintLabels. + /// @param rConstraintGaps Constants of the constraint equations. The constraint gap vector's size must match + /// the number of rows in @p rRelationMatrix. + LinearMultifreedomConstraint(const IndexType Id, + DofPointerVectorType&& rDofs, + const std::vector& rConstraintLabels, + const MatrixType& rRelationMatrix, + const VectorType& rConstraintGaps); + + void CalculateLocalSystem(MatrixType& rRelationMatrix, + VectorType& rConstraintGaps, + const ProcessInfo& rProcessInfo) const override; + + /// @name Unsupported Interface + /// @{ + + MasterSlaveConstraint::Pointer + Create(IndexType, + DofPointerVectorType&, + DofPointerVectorType&, + const MatrixType&, + const VectorType&) const override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported.";} + + MasterSlaveConstraint::Pointer + Create(IndexType, + NodeType&, + const VariableType&, + NodeType&, + const VariableType&, + const double, + const double) const override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported.";} + + /// @} +}; // class LinearMultifreedomConstraint + + +} // namespace Kratos \ No newline at end of file diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.cpp new file mode 100644 index 000000000000..826f3b9197d1 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.cpp @@ -0,0 +1,345 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.hpp" // MasterSlaveConstraintAssembler +#include "solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp" // GetDiagonalScaleFactor +#include "solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp" // MakeSparseTopology +#include "solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp" // ParseDiagonalScaling, GetDiagonalScaleFactor +#include "includes/define.h" // KRATOS_TRY, KRATOS_CATCH +#include "spaces/ublas_space.h" // TUblasSparseSpace +#include "utilities/sparse_matrix_multiplication_utility.h" // SparseMatrixMultiplicationUtility +#include "utilities/atomic_utilities.h" // AtomicAdd + + +namespace Kratos { + + +template +struct MasterSlaveConstraintAssembler::Impl +{ + std::vector mSlaveIds; + + std::vector mMasterIds; + + std::unordered_set mInactiveSlaveIds; + + std::unique_ptr mpDiagonalScaling; + + int mVerbosity; +}; + + +template +MasterSlaveConstraintAssembler::MasterSlaveConstraintAssembler() noexcept + : MasterSlaveConstraintAssembler(Parameters()) +{} + + +template +MasterSlaveConstraintAssembler::MasterSlaveConstraintAssembler(Parameters Settings) + : MasterSlaveConstraintAssembler(Settings, "unnamed") +{} + + +template +MasterSlaveConstraintAssembler::MasterSlaveConstraintAssembler(Parameters Settings, + std::string&& rInstanceName) + : Base(ConstraintImposition::MasterSlave, std::move(rInstanceName)), + mpImpl(new Impl) +{ + KRATOS_TRY + // Parse diagonal scaling and validate other settings. + Parameters default_parameters = this->GetDefaultParameters(); + Parameters default_diagonal_scaling = default_parameters["diagonal_scaling"].Clone(); + default_parameters.RemoveValue("diagonal_scaling"); + std::optional maybe_diagonal_scaling = Settings.Has("diagonal_scaling") ? Settings["diagonal_scaling"].Clone() : std::optional(); + if (maybe_diagonal_scaling.has_value()) Settings.RemoveValue("diagonal_scaling"); + Settings.ValidateAndAssignDefaults(default_parameters); + Settings.AddValue("diagonal_scaling", maybe_diagonal_scaling.has_value() ? maybe_diagonal_scaling.value() : default_diagonal_scaling); + mpImpl->mpDiagonalScaling = std::make_unique(Settings["diagonal_scaling"]); + mpImpl->mVerbosity = Settings["verbosity"].Get(); + KRATOS_CATCH("") +} + + +template +MasterSlaveConstraintAssembler::~MasterSlaveConstraintAssembler() = default; + + +template +void MasterSlaveConstraintAssembler::Allocate(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + typename Base::DofSet& rDofSet) +{ + KRATOS_TRY + std::vector> indices(rDofSet.size()); + std::vector mutexes(indices.size()); + const auto it_const_begin = rConstraints.begin(); + + #pragma omp parallel + { + Element::EquationIdVectorType slave_ids; + Element::EquationIdVectorType master_ids; + std::unordered_map> temp_indices; + + #pragma omp for nowait + for (int i_const = 0; i_const < static_cast(rConstraints.size()); ++i_const) { + auto it_const = it_const_begin + i_const; + it_const->EquationIdVector(slave_ids, master_ids, rProcessInfo); + + // Slave DoFs + for (auto &id_i : slave_ids) { + temp_indices[id_i].insert(master_ids.begin(), master_ids.end()); + } + } + + // Merging temporary indices into global rows. + for (auto& pair_temp_indices : temp_indices) { + std::scoped_lock lock(mutexes[pair_temp_indices.first]); + indices[pair_temp_indices.first].insert(pair_temp_indices.second.begin(), pair_temp_indices.second.end()); + } + } + + mpImpl->mSlaveIds.clear(); + mpImpl->mMasterIds.clear(); + for (int i = 0; i < static_cast(indices.size()); ++i) { + if (indices[i].size() == 0) // Master dof! + mpImpl->mMasterIds.push_back(i); + else // Slave dof + mpImpl->mSlaveIds.push_back(i); + indices[i].insert(i); // Ensure that the diagonal is there in T + } + + MakeSparseTopology(indices, + indices.size(), + this->GetRelationMatrix(), + /*EnsureDiagonal=*/false); + this->GetConstraintGapVector().resize(indices.size(), false); + + KRATOS_CATCH("") + + KRATOS_TRY + MergeMatrices(rLhs, this->GetRelationMatrix()); + KRATOS_CATCH("") +} + + +template +void MasterSlaveConstraintAssembler::Assemble(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename Base::DofSet& rDofSet, + const bool AssembleLhs, + const bool AssembleRhs) +{ + KRATOS_TRY + + /// @todo @p AssembleLhs and @p AssembleRhs are currently ignored, + /// and everything gets assembled. Think it through what needs + /// to be done and what can be omitted for each case. - matekelemen + + // Init. + mpImpl->mInactiveSlaveIds.clear(); + TSparse::SetToZero(this->GetRelationMatrix()); + TSparse::SetToZero(this->GetConstraintGapVector()); + std::vector mutexes(rDofSet.size()); + + // Declare local containers. + typename TDense::MatrixType local_relation_matrix; + typename TDense::VectorType local_constraint_gap_vector; + Element::EquationIdVectorType slave_equation_ids, master_equation_ids; + const int number_of_constraints = static_cast(rConstraints.size()); + + #pragma omp parallel firstprivate(local_relation_matrix, local_constraint_gap_vector, slave_equation_ids, master_equation_ids) + { + std::unordered_set tls_inactive_slave_dofs; + + #pragma omp for + for (int i_const = 0; i_const < number_of_constraints; ++i_const) { + auto it_const = rConstraints.begin() + i_const; + + // If the constraint is active + if (it_const->IsActive()) { + it_const->EquationIdVector(slave_equation_ids, + master_equation_ids, + rProcessInfo); + it_const->CalculateLocalSystem(local_relation_matrix, + local_constraint_gap_vector, + rProcessInfo); + + for (IndexType i = 0; i < slave_equation_ids.size(); ++i) { + const IndexType i_global = slave_equation_ids[i]; + + // Assemble matrix row. + { + std::scoped_lock lock(mutexes[i_global]); + MapRowContribution(this->GetRelationMatrix(), + local_relation_matrix, + i_global, + i, + master_equation_ids); + } + + // Assemble constant vector + KRATOS_ERROR_IF_NOT(i_global < this->GetConstraintGapVector().size()) + << "constraint gap vector size " << this->GetConstraintGapVector().size() + << " <= insertion index " << i_global << "\n"; + const double constant_value = local_constraint_gap_vector[i]; + AtomicAdd(this->GetConstraintGapVector()[i_global], typename TSparse::DataType(constant_value)); + } + } else { // Taking into account inactive constraints + tls_inactive_slave_dofs.insert(slave_equation_ids.begin(), slave_equation_ids.end()); + } + } + + // We merge all the sets in one thread + #pragma omp critical + { + mpImpl->mInactiveSlaveIds.insert(tls_inactive_slave_dofs.begin(), tls_inactive_slave_dofs.end()); + } + } + + // Setting the master dofs into the T and C system + for (auto eq_id : mpImpl->mMasterIds) { + std::scoped_lock lock(mutexes[eq_id]); + this->GetConstraintGapVector()[eq_id] = 0.0; + this->GetRelationMatrix()(eq_id, eq_id) = 1.0; + } + + // Setting inactive slave dofs in the T and C system + for (auto eq_id : mpImpl->mInactiveSlaveIds) { + std::scoped_lock lock(mutexes[eq_id]); + this->GetConstraintGapVector()[eq_id] = 0.0; + this->GetRelationMatrix()(eq_id, eq_id) = 1.0; + } + + KRATOS_CATCH("") +} + + +template +void MasterSlaveConstraintAssembler::Initialize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + [[maybe_unused]] typename Base::DofSet::iterator itDofBegin, + [[maybe_unused]] typename Base::DofSet::iterator itDofEnd) +{ + KRATOS_TRY + + // Output the relation matrix and constraint gaps if the verbosity is high enough. + if (4 <= mpImpl->mVerbosity) { + KRATOS_INFO(this->GetName() + ": ") + << "write relation matrix to " + << (std::filesystem::current_path() / "relation_matrix.mm") + << "\n"; + TSparse::WriteMatrixMarketMatrix("relation_matrix.mm", this->GetRelationMatrix(), false); + + KRATOS_INFO(this->GetName() + ": ") + << "write constraint gaps to " + << (std::filesystem::current_path() / "constraint_gaps.mm") + << "\n"; + TSparse::WriteMatrixMarketVector("constraint_gaps.mm", this->GetConstraintGapVector()); + } // if 4 <= verbosity + + // Compute the transposed matrix of the global relation matrix + { + // Storage for an intermediate matrix transpose(relation_matrix) * lhs_matrix + typename TSparse::MatrixType left_multiplied_lhs(this->GetRelationMatrix().size2(), rLhs.size2()); + { + // Transpose the relation matrix + typename TSparse::MatrixType transposed_relation_matrix(this->GetRelationMatrix().size2(), this->GetRelationMatrix().size1()); + SparseMatrixMultiplicationUtility::TransposeMatrix(transposed_relation_matrix, this->GetRelationMatrix(), 1.0); + + typename TSparse::VectorType b_modified(rRhs.size()); + BalancedProduct(transposed_relation_matrix, rRhs, b_modified); + TSparse::Copy(b_modified, rRhs); + + SparseMatrixMultiplicationUtility::MatrixMultiplication(transposed_relation_matrix, rLhs, left_multiplied_lhs); + } // deallocate transposed_relation_matrix + + SparseMatrixMultiplicationUtility::MatrixMultiplication(left_multiplied_lhs, this->GetRelationMatrix(), rLhs); + } // deallocate left_multiplied_lhs + + // Compute the scale factor for slave DoFs. + mpImpl->mpDiagonalScaling->template Cache(rLhs); + typename TSparse::DataType diagonal_scale_factor = mpImpl->mpDiagonalScaling->Evaluate(); + + // Apply diagonal values on slaves. + block_for_each(mpImpl->mSlaveIds, [this, &rLhs, &rRhs, diagonal_scale_factor](const auto iSlave){ + if (mpImpl->mInactiveSlaveIds.find(iSlave) == mpImpl->mInactiveSlaveIds.end()) { + rLhs(iSlave, iSlave) = diagonal_scale_factor; + rRhs[iSlave] = 0.0; + } + }); + + KRATOS_CATCH("") +} + + +template +bool +MasterSlaveConstraintAssembler::FinalizeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + PMGStatusStream::Report& rReport, + PMGStatusStream& rStream) +{ + rReport.maybe_constraint_residual = 0; + rReport.constraints_converged = true; + rStream.Submit(rReport.Tag(2), mpImpl->mVerbosity); + return true; +} + + +template +void MasterSlaveConstraintAssembler::Apply(typename TSparse::VectorType& rSolution) const +{ + KRATOS_TRY + typename TSparse::VectorType constrained_solution(rSolution.size()); + TSparse::SetToZero(constrained_solution); + BalancedProduct(this->GetRelationMatrix(), + rSolution, + constrained_solution); + rSolution.swap(constrained_solution); + KRATOS_CATCH("") +} + + +template +void MasterSlaveConstraintAssembler::Clear() +{ + Base::Clear(); + mpImpl->mSlaveIds = decltype(mpImpl->mSlaveIds)(); + mpImpl->mMasterIds = decltype(mpImpl->mMasterIds)(); + mpImpl->mInactiveSlaveIds = decltype(mpImpl->mInactiveSlaveIds)(); +} + + +template +Parameters MasterSlaveConstraintAssembler::GetDefaultParameters() +{ + return Parameters(R"({ + "method" : "master_slave", + "diagonal_scaling" : "norm", + "verbosity" : 1 + })"); +} + + +template class MasterSlaveConstraintAssembler,TUblasDenseSpace>; + +template class MasterSlaveConstraintAssembler,TUblasDenseSpace>; + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.hpp new file mode 100644 index 000000000000..fa206144cd06 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/master_slave_constraint_assembler.hpp @@ -0,0 +1,80 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project Includes +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp" // ConstraintAssembler +#include "includes/kratos_parameters.h" // Parameters + + +namespace Kratos { + + +template +class MasterSlaveConstraintAssembler final : public ConstraintAssembler +{ +public: + using Base = ConstraintAssembler; + + MasterSlaveConstraintAssembler() noexcept; + + MasterSlaveConstraintAssembler(Parameters Settings); + + MasterSlaveConstraintAssembler(Parameters Settings, + std::string&& rInstanceName); + + ~MasterSlaveConstraintAssembler(); + + /// @copydoc Base::Allocate + void Allocate(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + typename Base::DofSet& rDofSet) override; + + /// @copydoc Base::Assemble + void Assemble(const typename Base::ConstraintArray& rConstraints, + const ProcessInfo& rProcessInfo, + typename Base::DofSet& rDofSet, + const bool AssembleLhs, + const bool AssembleRhs) override; + + /// @copydoc Base::Initialize + void Initialize(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + typename Base::DofSet::iterator itDofBegin, + typename Base::DofSet::iterator itDofEnd) override; + + /// @copydoc Base::FinalizeSolutionStep + bool FinalizeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + PMGStatusStream::Report& rReport, + PMGStatusStream& rStream) override; + + /// @copydoc Base::Apply + void Apply(typename TSparse::VectorType& rSolution) const override; + + /// @copydoc Base::Clear + void Clear() override; + + static Parameters GetDefaultParameters(); + +private: + struct Impl; + std::unique_ptr mpImpl; +}; // class MasterSlaveConstraintAssembler + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.cpp new file mode 100644 index 000000000000..9c7a12c28ad1 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.cpp @@ -0,0 +1,90 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.hpp" // MultifreedomConstraint +#include "includes/variables.h" // CONSTRAINT_LABELS, GEOMETRIC_STIFFNESS_MATRIX + +// System includes +#include // std::copy, std::transform + + +namespace Kratos { + + +MultifreedomConstraint::MultifreedomConstraint(const IndexType Id, + DofArray&& rDofs, + const std::vector& rConstraintLabels) + : MasterSlaveConstraint(Id), + mDofs(std::move(rDofs)) +{ + Vector constraint_labels(rConstraintLabels.size()); + std::copy(rConstraintLabels.begin(), rConstraintLabels.end(), constraint_labels.begin()); + this->SetValue(CONSTRAINT_LABELS, constraint_labels); + this->SetValue(GEOMETRIC_STIFFNESS_MATRIX, ZeroMatrix(0)); +} + + +void MultifreedomConstraint::GetDofList(DofPointerVectorType& rSlaveDofs, + DofPointerVectorType& rMasterDofs, + const ProcessInfo&) const +{ + rSlaveDofs.clear(); + rMasterDofs.resize(mDofs.size()); + std::copy(mDofs.begin(), + mDofs.end(), + rMasterDofs.begin()); +} + + +void MultifreedomConstraint::SetDofList(const DofPointerVectorType& rSlaveDofs, + const DofPointerVectorType& rMasterDofs, + const ProcessInfo&) +{ + mDofs.resize(rSlaveDofs.size() + rMasterDofs.size()); + std::copy(rSlaveDofs.begin(), + rSlaveDofs.end(), + mDofs.begin()); + std::copy(rMasterDofs.begin(), + rMasterDofs.end(), + mDofs.begin() + rSlaveDofs.size()); +} + + +void MultifreedomConstraint::EquationIdVector(EquationIdVectorType& rSlaveDofs, + EquationIdVectorType& rMasterDofs, + const ProcessInfo& rProcessInfo) const +{ + rSlaveDofs.clear(); + rMasterDofs.resize(mDofs.size()); + std::transform(mDofs.begin(), + mDofs.end(), + rMasterDofs.begin(), + [](const Dof* p_dof) {return p_dof->EquationId();}); +} + + +void MultifreedomConstraint::save(Serializer& rSerializer) const +{ + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, MasterSlaveConstraint); + rSerializer.save("dofs", mDofs); +} + + +void MultifreedomConstraint::load(Serializer& rSerializer) +{ + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, MasterSlaveConstraint); + rSerializer.load("dofs", mDofs); +} + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.hpp new file mode 100644 index 000000000000..a50576935f89 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/multifreedom_constraint.hpp @@ -0,0 +1,174 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "includes/master_slave_constraint.h" // MasterSlaveConstraint +#include "includes/code_location.h" // KRATOS_CODE_LOCATION +#include "includes/variables.h" // GEOMETRIC_STIFFNESS_MATRIX + +// System includes +#include // std::vector +#include // std::numeric_limits + + +namespace Kratos { + + +/** @brief Class representing (part of) a linearizable constraint equation on @ref Dof "DoFs". + * @details MultifreedomConstraint offers a slightly different interface than @ref MasterSlaveConstraint, + * which it inherits from. The key difference is that @p MultifreedomConstraint does not impose + * a partitioning of its DoFs into slaves and masters. + * + * Another important difference is how constraint equations are identified. When a constraint + * equations is broken up into several @p MasterSlaveConstraints, instances belonging to the + * same constraint equation are the ones that have identical slave DoF IDs. + * Since @p MultifreedomConstraint provides no slave DoFs, each instance stores an array of + * constraint equation identifiers. This vector is stored in the @ref CONSTRAINT_LABELS variable + * of the @ref DataValueContainer belonging to the constraint instance. + * @code + * MultifreedomConstraint& r_constraint = ...; + * auto& r_constraint_ids = r_constraint[CONSTRAINT_LABELS]; + * @endcode + * Instances sharing constraint labels belong to the same constraint equation. + * + * @note Since this class does not partition its DoFs into slaves and masters, any member function returning + * data related to slaves/masters will leave slaves empty and only populate masters. + * - @ref MultifreedomConstraint::GetDofList + * - @ref MultifreedomConstraint::SetDofList + * - @ref MultifreedomConstraint::EquationIdVector + */ +class KRATOS_API(KRATOS_CORE) MultifreedomConstraint + : public MasterSlaveConstraint +{ +public: + KRATOS_CLASS_POINTER_DEFINITION(MultifreedomConstraint); + + using MasterSlaveConstraint::IndexType; + + using MasterSlaveConstraint::VariableType; + + using MasterSlaveConstraint::DofPointerVectorType; + + using MasterSlaveConstraint::MatrixType; + + using MasterSlaveConstraint::VectorType; + + using DofArray = std::vector*>; + + MultifreedomConstraint() noexcept + : MultifreedomConstraint(std::numeric_limits::max()) + {} + + explicit MultifreedomConstraint(const IndexType Id) noexcept + : MultifreedomConstraint(Id, DofPointerVectorType {}, std::vector {}) + {} + + /** @brief Construct with topological data. + * @param Id Identifier of the constraint instance. Note that, in general, this ID + * is not related to the ID of the constraint this instance is part of. + * @param rDofs Pointers to the @ref Dof "DoFs" this constraint is defined on. + * @param rConstraintLabels Identifiers of the constraint equations this instance is part of. + * Instances sharing constraint labels are assumed to be part of the + * same constraint equations, and are summed up during assembly. + */ + MultifreedomConstraint(const IndexType Id, + DofArray&& rDofs, + const std::vector& rConstraintLabels); + + /// @copydoc MasterSlaveConstraint::GetDofList + /// @note Slave DoFs are left empty and only master DoFs are populated. + void GetDofList(DofPointerVectorType& rSlaveDofs, + DofPointerVectorType& rMasterDofs, + const ProcessInfo&) const override; + + /// @copydoc MasterSlaveConstraint::SetDofList + /// @note Slave DoFs are appended to the list of master DoFs. + void SetDofList(const DofPointerVectorType& rSlaveDofs, + const DofPointerVectorType& rMastertDofs, + const ProcessInfo&) override; + + /// @copydoc MasterSlaveConstraint::EquationIdVector + /// @note Only master DoFs are populated. + void EquationIdVector(EquationIdVectorType& rSlaveDoFs, + EquationIdVectorType& rMasterDoFs, + const ProcessInfo& rProcessInfo) const override; + + /// @brief copydoc MasterSlaveConstraint::GetSlaveDofsVector + const DofPointerVectorType& GetSlaveDofsVector() const override + {return mDummy;} + + const DofPointerVectorType& GetMasterDofsVector() const override + {return mDofs;} + + void SetMasterDofsVector(const DofPointerVectorType& rDofs) override + {mDofs = rDofs;} + + + /// This function is NoOp. + void ResetSlaveDofs(const ProcessInfo&) override {} + + /// @copydoc MasterSlaveConstraint::Apply + void Apply(const ProcessInfo&) override {} + + /// @name Unsupported Interface + /// @{ + + MasterSlaveConstraint::Pointer + Create(IndexType, + DofPointerVectorType&, + DofPointerVectorType&, + const MatrixType&, + const VectorType&) const override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported.";} + + MasterSlaveConstraint::Pointer + Create(IndexType, + NodeType&, + const VariableType&, + NodeType&, + const VariableType&, + const double, + const double) const override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not supported.";} + + /// @} + +protected: + const DofArray& GetDofs() const noexcept + { + return mDofs; + } + + const MatrixType& GetHessian() const + { + return this->GetData().GetValue(GEOMETRIC_STIFFNESS_MATRIX); + } + + MatrixType& GetHessian() + { + return this->Data().GetValue(GEOMETRIC_STIFFNESS_MATRIX); + } + +private: + friend class Serializer; + + void save(Serializer& rSerializer) const override; + + void load(Serializer& rDeserializer) override; + + DofArray mDofs, mDummy; +}; // class MultifreedomConstraint + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/noop_constraint_assembler.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/noop_constraint_assembler.hpp new file mode 100644 index 000000000000..aa122f96dd04 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/noop_constraint_assembler.hpp @@ -0,0 +1,62 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project Includes +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp" // ConstraintAssembler +#include "includes/kratos_parameters.h" // Parameters + + +namespace Kratos { + + +/// @brief Dummy constraint assembler class that does not impose constraints. +template +class NoOpConstraintAssembler final : public ConstraintAssembler +{ +public: + using Base = ConstraintAssembler; + + NoOpConstraintAssembler() + : NoOpConstraintAssembler(Parameters()) + {} + + NoOpConstraintAssembler(Parameters Settings) + : NoOpConstraintAssembler(Settings, "unnamed") + { + } + + NoOpConstraintAssembler(Parameters Settings, std::string&& rInstanceName) noexcept + : Base(ConstraintImposition::None, std::move(rInstanceName)), + mVerbosity(Settings.Has("verbosity") ? Settings["verbosity"].Get() : 1) + {} + + bool FinalizeSolutionStep(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rSolution, + typename TSparse::VectorType& rRhs, + PMGStatusStream::Report& rReport, + PMGStatusStream& rStream) override + { + rReport.maybe_constraint_residual = 0; + rReport.constraints_converged = true; + rStream.Submit(rReport.Tag(2), mVerbosity); + return true; + } + +private: + int mVerbosity; +}; // class NoOpConstraintAssembler + + +} // namespace Kratos + diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.cpp new file mode 100644 index 000000000000..280d9dbb5b37 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.cpp @@ -0,0 +1,626 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/p_grid.hpp" // PGrid +#include "solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_utilities.hpp" // MakePRestrictionOperator +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.hpp" // ConstraintAssemblerFactory +#include "solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp" // ApplyDirichletConditions +#include "solving_strategies/builder_and_solvers/p_multigrid/status_stream.hpp" +#include "spaces/ublas_space.h" // TUblasSparseSpace, TUblasDenseSpace +#include "factories/linear_solver_factory.h" // LinearSolverFactory +#include "includes/kratos_components.h" // KratosComponents +#include "utilities/profiler.h" // KRATOS_PROFILE_SCOPE +#include "utilities/sparse_matrix_multiplication_utility.h" // SparseMatrixMultiplicationUtility +#include "utilities/builtin_timer.h" // BuiltinTimer + +// System includes +#include // std::numeric_limits + + +namespace Kratos { + + +template +PGrid::PGrid(Parameters Settings, + const unsigned CurrentDepth, + Parameters SmootherSettings, + Parameters LeafSolverSettings, + Parameters DiagonalScalingSettings) + : mRestrictionOperator(), + mProlongationOperator(), + mLhs(), + mSolution(), + mRhs(), + mDofSet(), + mIndirectDofSet(), + mDofMap(), + mpConstraintAssembler(), + mpDiagonalScaling(), + mpSolver(), + mMaybeChild(), + mVerbosity(), + mDepth(CurrentDepth) +{ + // Sanity checks. + KRATOS_ERROR_IF_NOT(mDepth) << "PGrid can only have positive depths (the original system is at depth=0)"; + + KRATOS_TRY + Settings.ValidateAndAssignDefaults(this->GetDefaultParameters()); + + // Check float precision. + using value_type = typename TSparse::DataType; + const std::string precision = Settings["precision"].Get(); + if constexpr (std::is_same_v) { + KRATOS_ERROR_IF_NOT(precision == "double") + << "attempting to construct a PGrid with inconsistent sparse space type " + << "(requested precision: \"" << precision << "\" " + << " build precision: \"double\")"; + } else if constexpr (std::is_same_v) { + KRATOS_ERROR_IF_NOT(precision == "single") + << "attempting to construct a PGrid with inconsistent sparse space type " + << "(requested precision: \"" << precision << "\" " + << " build precision: \"single\")"; + } else { + static_assert(!std::is_same_v, "unhandled sparse space type"); + } + + mpConstraintAssembler = ConstraintAssemblerFactory(Settings["constraint_imposition_settings"], + "Grid " + std::to_string(mDepth) + " constraints"); + mpDiagonalScaling = std::make_unique(DiagonalScalingSettings); + mVerbosity = Settings["verbosity"].Get(); + + const int max_depth = Settings["max_depth"].Get(); + KRATOS_ERROR_IF_NOT(0 <= max_depth) << Settings << "\n\"max_depth\" must be non-negative"; + if (mDepth < static_cast(max_depth)) { + KRATOS_TRY + KRATOS_ERROR_IF_NOT(SmootherSettings.Has("solver_type")); + const std::string solver_name = SmootherSettings["solver_type"].Get(); + using SolverFactoryRegistry = KratosComponents>; + KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name)) + << "\"" << solver_name << "\" is not a valid linear solver name in the registry. " + << "Make sure you imported the application it is defined in and that the spelling is correct."; + const auto& r_factory = SolverFactoryRegistry::Get(solver_name); + mpSolver = r_factory.Create(SmootherSettings); + mMaybeChild = std::unique_ptr(new PGrid(Settings, + mDepth + 1u, + SmootherSettings, + LeafSolverSettings, + DiagonalScalingSettings)); + KRATOS_CATCH(std::to_string(mDepth + 1u)); + } else { + KRATOS_ERROR_IF_NOT(LeafSolverSettings.Has("solver_type")); + const std::string solver_name = LeafSolverSettings["solver_type"].Get(); + using SolverFactoryRegistry = KratosComponents>; + + if (!SolverFactoryRegistry::Has(solver_name)) { + std::stringstream message; + message << "PMultigridBuilderAndSolver: " + << "\"" << solver_name << "\" is not a valid linear solver name in the registry. " + << "Make sure you imported the application it is defined in and that the spelling is correct. " + << "Registered options are:\n"; + for ([[maybe_unused]] const auto& [r_name, r_entry] : SolverFactoryRegistry::GetComponents()) { + message << "\t" << r_name << "\n"; + } + KRATOS_ERROR << message.str(); + } // if not SolverFactoryRegistry::Has(solver_name) + + const auto& r_factory = SolverFactoryRegistry::Get(solver_name); + mpSolver = r_factory.Create(LeafSolverSettings); + } + KRATOS_CATCH("") +} + + +template +PGrid::PGrid(Parameters Settings, + Parameters SmootherSettings, + Parameters LeafSolverSettings, + Parameters DiagonalScalingSettings) + : PGrid(Settings, + 1u, + SmootherSettings, + LeafSolverSettings, + DiagonalScalingSettings) +{ +} + + +template +PGrid::PGrid() + : PGrid(Parameters(), + Parameters(R"({"solver_type" : "gauss_seidel"})"), + Parameters(R"({"solver_type" : "amgcl"})"), + Parameters(R"("norm")")) +{ +} + + +template +template +void PGrid::MakeLhsTopology(ModelPart& rModelPart, + const typename TParentSparse::MatrixType& rParentLhs, + [[maybe_unused]] const ConstraintAssembler& rParentConstraintAssembler, + const IndirectDofSet& rParentDofSet) +{ + // ConstraintAssembler::Allocate is deliberately not invoked here because the relation matrix + // and constraint gap vector are passed from the fine level in PGrid::Assemble, and mapped + // to this level by the restriction operator. +} + + +template +template +void PGrid::Assemble(ModelPart& rModelPart, + const typename TParentSparse::MatrixType* pParentLhs, + const typename TParentSparse::VectorType* pParentRhs, + const ConstraintAssembler& rParentConstraintAssembler, + IndirectDofSet& rParentDofSet) +{ + KRATOS_TRY + using SparseUtils = SparseMatrixMultiplicationUtility; + + // Assemble LHS matrix. + if constexpr (AssembleLHS) { + BuiltinTimer timer; + KRATOS_ERROR_IF_NOT(pParentLhs); + + // The restriction operator immediately constructs the linear equivalent, + // because arbitrary-level coarsening strategies are not supported yet. The problem + // is that when going from a generic polynomial level q to some other lower polynomial + // level p!=1, new nodes must be introduced that do not exist in the original fine mesh. + // This wouldn't be a huge problem by itself, but deciding on DoF indices on the coarse + // level would be quite painful and require keeping track of the coarse grid's topology + // in some manner. + // + // Note: + // One might argue that the construction of the restriction operator could happen + // during the allocation stage, but the problem is that some constraint imposition + // methods (such as lagrange or its augmented version) may mutate the DoFs during their + // assembly. As a result, the restriction operator implicitly depends on the constraint + // imposition method of the fine grid, meaning it must be constructed AFTER constraint + // assembly. + MakePRestrictionOperator::max(),typename TSparse::DataType>( + rModelPart, + pParentLhs->size1(), + rParentDofSet, + mRestrictionOperator, + mpVariableList, + mDofSet, + mIndirectDofSet, + mDofMap); + + // Compute the coarse LHS matrix. + typename TSparse::MatrixType left_multiplied_lhs; + SparseUtils::MatrixMultiplication(mRestrictionOperator, *pParentLhs, left_multiplied_lhs); + SparseUtils::TransposeMatrix(mProlongationOperator, mRestrictionOperator, 1.0); + SparseUtils::MatrixMultiplication(left_multiplied_lhs, mProlongationOperator, mLhs); + + KRATOS_INFO_IF("Grid " + std::to_string(mDepth), 2 <= this->mVerbosity) + << ": grid construction took " << timer << "\n"; + } // if AssembleLHS + + // Assemble RHS vector. + if constexpr (AssembleRHS) { + // Sanity checks + KRATOS_ERROR_IF_NOT(pParentRhs); + KRATOS_ERROR_IF_NOT(mRestrictionOperator.size2() == pParentRhs->size()) + << "expecting an RHS vector of size " << mRestrictionOperator.size2() + << " but got " << pParentRhs->size(); + + // Allocate the coarse vectors. + mSolution.resize(mRestrictionOperator.size1(), false); + mRhs.resize(mRestrictionOperator.size1(), false); + } // if AssembleRHS + + // Compute coarse constraints. + // How this is actually implemented unfortunately depends on what the imposition + // method is. The original master-slave imposition stores a different version of + // the relation matrix, while the noop imposition stores nothing. Since there are + // two grid levels here that may use different imposition methods, all combinations + // must be handled separately. + switch (rParentConstraintAssembler.GetImposition()) { + // The parent does not impose constraints, so neither must the current level. + case ConstraintImposition::None: + KRATOS_ERROR_IF_NOT(mpConstraintAssembler->GetImposition() == ConstraintImposition::None) + << "PMultigridBuilderAndSolver: grid " << mDepth + << " imposes constraints (" << mpConstraintAssembler->GetValue(mpConstraintAssembler->GetImpositionVariable()) << ")" + << " but its parent does not"; + break; + + // The parent imposes constraints via master-slave elimination. The original + // implementation stores a transformation matrix instead of the relation matrix, + // so coarsening it is complicated and expensive => no support for it for now. + /// @todo implement constraint imposition when the parent on a p-grid uses + // master-slave imposition but the child grid uses augmented lagrange. + case ConstraintImposition::MasterSlave: + KRATOS_ERROR << "PMultigridBuilderAndSolver: constraint imposition using master-slave elimination " + << "is only supported if there is no coarse hierarchy (depth=0). Consider imposition " + << "using augmented lagrange multipliers."; + break; + + // The parent imposes constraints via augmented lagrange multipliers. + // The only supported imposition on the current level in this case is also + // augmented lagrange (for similar reasons why master-slave doesn't support anything). + // There's also a theoretical reason behind the choice of dropping support for + // master-slave elimination on this level though. Picking augmented lagrange on the + // parent grid probably means that the user expects that constraints may become + // linearly dependent at some point. If that's the case, master-slave elimination will + // definitely break the linear solver while augmented lagrange might provide acceptable + // results if configured for penalty. + case ConstraintImposition::AugmentedLagrange: + if (AssembleLHS) { + SparseUtils::MatrixMultiplication(rParentConstraintAssembler.GetRelationMatrix(), + mProlongationOperator, + mpConstraintAssembler->GetRelationMatrix()); + + typename TSparse::MatrixType tmp; + SparseUtils::MatrixMultiplication(mRestrictionOperator, + rParentConstraintAssembler.GetHessian(), + tmp); + SparseUtils::MatrixMultiplication(tmp, + mProlongationOperator, + mpConstraintAssembler->GetHessian()); + } + + if (AssembleRHS) { + mpConstraintAssembler->GetConstraintGapVector() = rParentConstraintAssembler.GetConstraintGapVector(); + } + break; + + // No other impositions are supported for now. + default: + KRATOS_ERROR << "PMultigridBuilderAndSolver: unsupported constraint imposition at depth " << mDepth + << " (parent: " << rParentConstraintAssembler.GetValue(rParentConstraintAssembler.GetImpositionVariable()) + << " child: " << mpConstraintAssembler->GetValue(mpConstraintAssembler->GetImpositionVariable()) << ")"; + } // switch rParentConstraintAssembler.GetImposition() + + KRATOS_CATCH("") +} + + +template +void PGrid::ApplyDirichletConditions(typename IndirectDofSet::const_iterator itParentDofBegin, + [[maybe_unused]] typename IndirectDofSet::const_iterator itParentDofEnd) +{ + + if (mIndirectDofSet.empty()) return; + + // Apply dirichlet conditions on the restriction operator. + KRATOS_TRY + block_for_each(mIndirectDofSet.begin(), + mIndirectDofSet.end(), + [this, itParentDofBegin](const Dof& r_dof){ + const std::size_t i_dof = r_dof.EquationId(); + const typename TSparse::IndexType i_entry_begin = mRestrictionOperator.index1_data()[i_dof]; + const typename TSparse::IndexType i_entry_end = mRestrictionOperator.index1_data()[i_dof + 1]; + + if (r_dof.IsFixed()) { + // Zero out the whole row, except the entry related to the dof on the fine grid. + const auto i_fine_dof = mDofMap[i_dof]; + for (typename TSparse::IndexType i_entry=i_entry_begin; i_entry(1); + } else { + mRestrictionOperator.value_data()[i_entry] = static_cast(0); + } + } // for i_entry in range(i_entry_begin, i_entry_end) + } /*if r_dof.IsFixed()*/ else { + // Zero out the column which is associated with the zero'ed row. + for (typename TSparse::IndexType i_entry=i_entry_begin; i_entryIsFixed()) { + mRestrictionOperator.value_data()[i_entry] = 0.0; + } + } // for i_entry in range(i_entry_begin, i_entry_end) + } /*not r_dof.IsFixed()*/ + }); + KRATOS_CATCH("") + + // Apply dirichlet conditions on the prolongation operator. + // @todo make this more efficient. + KRATOS_TRY + mProlongationOperator = decltype(mProlongationOperator)(); + SparseMatrixMultiplicationUtility::TransposeMatrix(mProlongationOperator, mRestrictionOperator, 1.0); + KRATOS_CATCH("") + + // Apply dirichlet conditions on the LHS. + KRATOS_TRY + mpDiagonalScaling->Cache(mLhs); + const auto diagonal_scale = mpDiagonalScaling->Evaluate(); + Kratos::ApplyDirichletConditions( + mLhs, + mRhs, + mIndirectDofSet.begin(), + mIndirectDofSet.end(), + diagonal_scale); + KRATOS_CATCH("") +} + + +template +void PGrid::ApplyConstraints() +{ + KRATOS_TRY + mpConstraintAssembler->Initialize(mLhs, + mRhs, + mIndirectDofSet.begin(), + mIndirectDofSet.end()); + KRATOS_CATCH("") +} + + +template +template +void PGrid::Initialize(ModelPart& rModelPart, + const typename TParentSparse::MatrixType&, + const typename TParentSparse::VectorType&, + const typename TParentSparse::VectorType&) +{ + KRATOS_TRY + if (mpSolver->AdditionalPhysicalDataIsNeeded()) + mpSolver->ProvideAdditionalData(mLhs, + mSolution, + mRhs, + mIndirectDofSet, + rModelPart); + KRATOS_CATCH("") +} + + +template +void PGrid::ExecuteMultigridLoop(PMGStatusStream& rStream, + PMGStatusStream::Report& rReport) +{ + KRATOS_TRY + + rReport.multigrid_converged = false; + rReport.multigrid_iteration = 0ul; + rReport.maybe_multigrid_residual.reset(); + + // The multigrid hierarchy depth is currently capped at 1, + // so the linear solver is used here instead of invoking + // lower grids. + mpSolver->InitializeSolutionStep(mLhs, mSolution, mRhs); + rReport.multigrid_converged = mpSolver->Solve(mLhs, mSolution, mRhs); + mpSolver->FinalizeSolutionStep(mLhs, mSolution, mRhs); + + KRATOS_CATCH("") +} + + +template +void PGrid::ExecuteConstraintLoop(PMGStatusStream& rStream, + PMGStatusStream::Report& rReport) +{ + bool constraints_finished = false; + auto initial_residual = TSparse::TwoNorm(mRhs); + initial_residual = initial_residual ? initial_residual : 1; + + // Impose constraints and solve the coarse system. + KRATOS_TRY + do { + rReport.constraints_converged = false; + rReport.maybe_constraint_residual.reset(); + + // Initialize the constraint assembler. + mpConstraintAssembler->InitializeSolutionStep(mLhs, mSolution, mRhs); + + // Get an update on the solution with respect to the current right hand side. + this->ExecuteMultigridLoop(rStream, rReport); + + constraints_finished = mpConstraintAssembler->FinalizeSolutionStep(mLhs, mSolution, mRhs, rReport, rStream); + + // Update state log. + if (!constraints_finished) { + rStream.Submit(rReport.Tag(3), mVerbosity); + ++rReport.constraint_iteration; + } + } while (!constraints_finished); + + // Update the residual. + BalancedProduct(mLhs, mSolution, mRhs, static_cast(-1)); + + // Update state log. + rReport.maybe_multigrid_residual = TSparse::TwoNorm(mRhs) / initial_residual; + rStream.Submit(rReport.Tag(2), mVerbosity); + + mpConstraintAssembler->Finalize(mLhs, mSolution, mRhs, mIndirectDofSet); + KRATOS_CATCH("") +} + + +template +template +bool PGrid::ApplyCoarseCorrection(typename TParentSparse::VectorType& rParentSolution, + const typename TParentSparse::VectorType& rParentRhs, + const ConstraintAssembler& rParentConstraintAssembler, + PMGStatusStream& rStream) +{ + #ifndef NDEBUG + KRATOS_TRY + CheckMatrix(mRestrictionOperator); + KRATOS_CATCH("") + KRATOS_TRY + CheckMatrix(mProlongationOperator); + KRATOS_CATCH("") + #endif + + PMGStatusStream::Report status_report { + /*grid_level=*/ static_cast(mDepth), + /*multigrid_converged=*/ false, + /*multigrid_iteration=*/ 0ul, + /*multigrid_residual=*/ {}, + /*constraints_converged=*/ false, + /*constraint_iteration=*/ 0ul, + /*maybe_constraint_residual=*/ {} + }; + + TSparse::SetToZero(mSolution); + this->Restrict(rParentRhs, mRhs); + this->ExecuteConstraintLoop(rStream, status_report); + this->Prolong(mSolution, rParentSolution); + + return status_report.multigrid_converged && status_report.constraints_converged; +} + + +template +template +void PGrid::Finalize(ModelPart& rModelPart, + const typename TParentSparse::MatrixType&, + const typename TParentSparse::VectorType&, + const typename TParentSparse::VectorType&) +{ + mpConstraintAssembler->Finalize(mLhs, mSolution, mRhs, mIndirectDofSet); +} + + +template +template +void PGrid::Restrict(const typename TParentSparse::VectorType& rParentVector, + typename TSparse::VectorType& rOutput) const +{ + KRATOS_TRY + rOutput.resize(mRestrictionOperator.size1(), false), + TSparse::SetToZero(rOutput); + BalancedProduct(mRestrictionOperator, + rParentVector, + rOutput); + KRATOS_CATCH("") +} + + +template +template +void PGrid::Prolong(const typename TSparse::VectorType& rInput, + typename TParentSparse::VectorType& rParentVector) const +{ + KRATOS_TRY + rParentVector.resize(mProlongationOperator.size1(), false); + TParentSparse::SetToZero(rParentVector); + BalancedProduct(mProlongationOperator, + rInput, + rParentVector); + KRATOS_CATCH("") +} + + +template +void PGrid::Clear() +{ + mRestrictionOperator = decltype(mRestrictionOperator)(); + mProlongationOperator = decltype(mProlongationOperator)(); + mLhs = decltype(mLhs)(); + mSolution = decltype(mSolution)(); + mRhs = decltype(mRhs)(); + mIndirectDofSet = decltype(mIndirectDofSet)(); + mDofSet = decltype(mDofSet)(); + mpConstraintAssembler.reset(); + mpSolver.reset(); + + if (mMaybeChild.has_value()) { + mMaybeChild.value().reset(); + } +} + + +template +Parameters PGrid::GetDefaultParameters() +{ + Parameters output = Parameters(R"({ + "max_depth" : 0, + "verbosity" : 1, + "precision" : "", + "constraint_imposition_settings" : { + "method" : "augmented_lagrange", + "max_iterations" : 1 + } + })"); + + using value_type = typename TSparse::DataType; + if constexpr (std::is_same_v) { + output["precision"].SetString("double"); + } else if constexpr (std::is_same_v) { + output["precision"].SetString("single"); + } else { + static_assert(!std::is_same_v, "unhandled sparse space type"); + } + + return output; +} + + +#define KRATOS_INSTANTIATE_PGRID_MEMBERS(TSparse, TDense, TParentSparse) \ + template void PGrid::MakeLhsTopology(ModelPart&, \ + const typename TParentSparse::MatrixType&, \ + const ConstraintAssembler&, \ + const PointerVectorSet>&); \ + template void PGrid::Assemble(ModelPart&, \ + const typename TParentSparse::MatrixType*, \ + const typename TParentSparse::VectorType*, \ + const ConstraintAssembler&, \ + PointerVectorSet>&); \ + template void PGrid::Assemble(ModelPart&, \ + const typename TParentSparse::MatrixType*, \ + const typename TParentSparse::VectorType*, \ + const ConstraintAssembler&, \ + PointerVectorSet>&); \ + template void PGrid::Assemble(ModelPart&, \ + const typename TParentSparse::MatrixType*, \ + const typename TParentSparse::VectorType*, \ + const ConstraintAssembler&, \ + PointerVectorSet>&); \ + template void PGrid::Assemble(ModelPart&, \ + const typename TParentSparse::MatrixType*, \ + const typename TParentSparse::VectorType*, \ + const ConstraintAssembler&, \ + PointerVectorSet>&); \ + template void PGrid::Initialize(ModelPart&, \ + const TParentSparse::MatrixType&, \ + const TParentSparse::VectorType&, \ + const TParentSparse::VectorType&); \ + template bool PGrid::ApplyCoarseCorrection(TParentSparse::VectorType&, \ + const TParentSparse::VectorType&, \ + const ConstraintAssembler&, \ + PMGStatusStream&); \ + template void PGrid::Finalize(ModelPart&, \ + const TParentSparse::MatrixType&, \ + const TParentSparse::VectorType&, \ + const TParentSparse::VectorType&); \ + template void PGrid::Restrict(const typename TParentSparse::VectorType&, \ + typename TSparse::VectorType&) const; \ + template void PGrid::Prolong(const typename TSparse::VectorType&, \ + typename TParentSparse::VectorType&) const + +#define KRATOS_INSTANTIATE_PGRID(TSparse, TDense) \ + template class PGrid; \ + KRATOS_INSTANTIATE_PGRID_MEMBERS(TSparse, TDense, TUblasSparseSpace); \ + KRATOS_INSTANTIATE_PGRID_MEMBERS(TSparse, TDense, TUblasSparseSpace) + +KRATOS_INSTANTIATE_PGRID(TUblasSparseSpace, TUblasDenseSpace); + +KRATOS_INSTANTIATE_PGRID(TUblasSparseSpace, TUblasDenseSpace); + +#undef KRATOS_INSTANTIATE_GRID +#undef KRATOS_INSTANTIATE_GRID_MEMBERS + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.hpp new file mode 100644 index 000000000000..f604aef658b9 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_grid.hpp @@ -0,0 +1,194 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler.hpp" // ConstraintAssembler +#include "solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp" // DiagonalScaling +#include "solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_utilities.hpp" // detail::DofData +#include "linear_solvers/linear_solver.h" // LinearSolver, Reorderer +#include "includes/dof.h" // Dof + +// System includes +#include // std::shared_ptr +#include // std::optional +#include // std::vector + + +namespace Kratos { + + +/** @brief Coarse grid level of @PMultigridBuilderAndSolver. + * @details Settings: + * @code + * { + * "max_depth" : 0, + * "precision" : "double", + * "verbosity" : 1, + * "constraint_imposition_settings" : { + * "method" : "augmented_lagrange", + * "max_iterations" : 1 + * } + * } + * @endcode + * - @p "max_depth" Maximum number of coarse grid levels. If set to 0 (default), there are no + * coarse grids and no multigrid preconditioning is done. This also means that + * that the top-level grid uses the linear solver instead of the smoother. + * - @p "precision" Floating point precision of system matrices and vectors. Either @p "double" (default) + * or @p "single". The potential benefit of using single precision systems is reducing + * memory consumtion on accelerator hardware (e.g.: VRAM on GPUs). This is viable because + * coarse corrections in a p-multigrid preconditioner need not have high accuracy. + * - @p "verbosity" Level of verbosity. Every level includes lower verbosity levels and + * adds new events to report: + * - @p 0 No messages, not even in case of failure. + * - @p 1 Warnings and failure messages. + * - @p 2 Aggregated status reports. + * - @p 3 Per-iteration status reports. + * - @p 4 Output system matrices and vectors. + * - @p "constraint_imposition_settings" : Settings of the top-level constraint assembler. + * Refer to @ref ConstraintAssemblerFactory for more + * information. + */ +template +class PGrid +{ +public: + using LinearSolverType = LinearSolver>; + + using IndirectDofSet = PointerVectorSet>; + + PGrid(); + + PGrid(Parameters Settings, + Parameters SmootherSettings, + Parameters LeafSolverSettings, + Parameters DiagonalScalingSettings); + + PGrid(PGrid&&) noexcept = default; + + PGrid& operator=(PGrid&&) noexcept = default; + + template + void MakeLhsTopology(ModelPart& rModelPart, + const typename TParentSparse::MatrixType& rParentLhs, + const ConstraintAssembler& rParentConstraintAssembler, + const IndirectDofSet& rParentDofSet); + + template + void Assemble(ModelPart& rModelPart, + const typename TParentSparse::MatrixType* pParentLhs, + const typename TParentSparse::VectorType* pParentRhs, + const ConstraintAssembler& rParentConstraintAssembler, + IndirectDofSet& rParentDofSet); + + void ApplyDirichletConditions(typename IndirectDofSet::const_iterator itParentDofBegin, + typename IndirectDofSet::const_iterator itParentDofEnd); + + void ApplyConstraints(); + + template + void Initialize(ModelPart& rModelPart, + const typename TParentSparse::MatrixType& rParentLhs, + const typename TParentSparse::VectorType& rParentSolution, + const typename TParentSparse::VectorType& rParentRhs); + + template + bool ApplyCoarseCorrection(typename TParentSparse::VectorType& rParentSolution, + const typename TParentSparse::VectorType& rParentResidual, + const ConstraintAssembler& rParentConstraintAssembler, + PMGStatusStream& rStream); + + template + void Finalize(ModelPart& rModelPart, + const typename TParentSparse::MatrixType& rParentLhs, + const typename TParentSparse::VectorType& rParentSolution, + const typename TParentSparse::VectorType& rParentRhs); + + template + void Restrict(const typename TParentSparse::VectorType& rParentVector, + typename TSparse::VectorType& rOutput) const; + + template + void Prolong(const typename TSparse::VectorType& rInput, + typename TParentSparse::VectorType& rParentVector) const; + + void Clear(); + + static Parameters GetDefaultParameters(); + + std::optional GetChild() const + { + return mMaybeChild.has_value() ? mMaybeChild.value().get() : std::optional(); + } + + const typename TSparse::VectorType& GetSolution() const + { + return mSolution; + } + +private: + PGrid(Parameters Settings, + const unsigned CurrentDepth, + Parameters SmootherSettings, + Parameters LeafSolverSettings, + Parameters DiagonalScalingSettings); + + PGrid(const PGrid&) = delete; + + PGrid& operator=(const PGrid&) = delete; + + void ExecuteMultigridLoop(PMGStatusStream& rStream, + PMGStatusStream::Report& rReport); + + void ExecuteConstraintLoop(PMGStatusStream& rStream, + PMGStatusStream::Report& rReport); + + typename TSparse::MatrixType mRestrictionOperator; + + typename TSparse::MatrixType mProlongationOperator; + + typename TSparse::MatrixType mLhs; + + typename TSparse::VectorType mSolution; + + typename TSparse::VectorType mRhs; + + VariablesList::Pointer mpVariableList; + + /// @details Array of @ref Dof "DoFs" unique to the current grid level. + /// DoFs need a pointer to a @ref NodalData object, which is why + /// pairs of @ref NodalData and @ref Dof are stored instead of just + /// DoFs. + std::vector mDofSet; + + IndirectDofSet mIndirectDofSet; + + std::vector mDofMap; + + std::shared_ptr> mpConstraintAssembler; + + std::unique_ptr mpDiagonalScaling; + + typename LinearSolverType::Pointer mpSolver; + + std::optional> mMaybeChild; + + int mVerbosity; + + unsigned mDepth; +}; // class PGrid + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp new file mode 100644 index 000000000000..7d98e48cb71b --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.cpp @@ -0,0 +1,1201 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp" // PMultigridBuilderAndSolver +#include "solving_strategies/builder_and_solvers/p_multigrid/constraint_assembler_factory.hpp" // ConstraintAssemblerFactory +#include "solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp" // MakeSparseTopology +#include "solving_strategies/builder_and_solvers/p_multigrid/diagonal_scaling.hpp" // DiagonalScaling, ParseDiagonalScaling, GetDiagonalScalingFactor +#include "solving_strategies/builder_and_solvers/p_multigrid/p_grid.hpp" // PGrid +#include "solving_strategies/builder_and_solvers/p_multigrid/status_stream.hpp" // PMGStatusStream +#include "includes/model_part.h" // ModelPart +#include "spaces/ublas_space.h" // TUblasSparseSpace, TUblasDenseSpace +#include "linear_solvers/linear_solver.h" // LinearSolver +#include "factories/linear_solver_factory.h" // LinearSolverFactory +#include "includes/kratos_components.h" // KratosComponents +#include "utilities/proxies.h" // MakeProxy +#include "utilities/profiler.h" // KRATOS_PROFILE_SCOPE, KRATOS_PROFILE_SCOPE_MILLI +#include "utilities/builtin_timer.h" // BuiltinTimer + +// System includes +#include // std::optional +#include // std::unordered_set +#include // std::variant + + +namespace Kratos { + + +// --------------------------------------------------------- // +// PIMPL +// --------------------------------------------------------- // + + +template +struct PMultigridBuilderAndSolver::Impl +{ + using Interface = PMultigridBuilderAndSolver; + + // --------------------------------------------------------- // + // Member Variables + // --------------------------------------------------------- // + + Interface* mpInterface; + + std::optional mpMaybeModelPart; + + std::shared_ptr> mpConstraintAssembler; + + std::optional,TUblasDenseSpace>, + PGrid,TUblasDenseSpace> + >> mMaybeHierarchy; + + std::unique_ptr mpDiagonalScaling; + + int mMaxIterations; + + typename TSparse::DataType mTolerance; + + int mMaxDepth; + + int mVerbosity; + + // --------------------------------------------------------- // + // Special Member Functions + // --------------------------------------------------------- // + + Impl(Interface* pInterface) + : mpInterface(pInterface), + mpMaybeModelPart(), + mpConstraintAssembler(), + mMaybeHierarchy(), + mpDiagonalScaling(), + mMaxIterations(0), + mTolerance(std::numeric_limits::max()), + mMaxDepth(-1), + mVerbosity(0) + {} + + // --------------------------------------------------------- // + // Solution + // --------------------------------------------------------- // + + + void ExecuteMultigridLoop(typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + const typename Interface::TSystemVectorType& rRhs, + typename Interface::TSystemVectorType& rSolutionUpdate, + typename Interface::TSystemVectorType& rResidual, + const typename TSparse::DataType InitialResidualNorm, + PMGStatusStream& rStream, + PMGStatusStream::Report& rReport) + { + KRATOS_TRY + // Inner loop for multigrid. + rReport.multigrid_converged = false; + rReport.multigrid_iteration = 0ul; + + while ( rReport.multigrid_iteration < static_cast(mMaxIterations) + && !rReport.multigrid_converged) { + rReport.maybe_multigrid_residual.reset(); + + // Solve the coarse grid and apply its correction. + if (mMaybeHierarchy.has_value()) { + std::visit([&rSolutionUpdate, &rResidual, &rStream, this](auto& r_hierarchy){ + return r_hierarchy.template ApplyCoarseCorrection( + rSolutionUpdate, + rResidual, + *mpConstraintAssembler, + rStream);}, + mMaybeHierarchy.value()); + + // Update the fine solution. + TSparse::UnaliasedAdd(rSolution, 1.0, rSolutionUpdate); + + // Update the fine residual. + BalancedProduct(rLhs, rSolutionUpdate, rResidual, static_cast(-1)); + } // if mMaybeHierarchy + + // Perform smoothing on the fine grid. + TSparse::SetToZero(rSolutionUpdate); //< do I need this? + mpInterface->GetLinearSystemSolver()->InitializeSolutionStep(rLhs, rSolutionUpdate, rResidual); + mpInterface->GetLinearSystemSolver()->Solve(rLhs, rSolutionUpdate, rResidual); + mpInterface->GetLinearSystemSolver()->FinalizeSolutionStep(rLhs, rSolutionUpdate, rResidual); + + // Update the fine solution. + TSparse::UnaliasedAdd(rSolution, 1.0, rSolutionUpdate); + + // Update the fine residual. + BalancedProduct(rLhs, rSolutionUpdate, rResidual, static_cast(-1)); + + // Emit status and check for convergence. + rReport.maybe_multigrid_residual = TSparse::TwoNorm(rResidual) / InitialResidualNorm; + rReport.multigrid_converged = rReport.maybe_multigrid_residual.value() < mTolerance; + if ( rReport.multigrid_iteration + 1 < static_cast(mMaxIterations) + && !rReport.multigrid_converged) + rStream.Submit(rReport.Tag(3), mVerbosity); + + ++rReport.multigrid_iteration; + } // while i_multigrid_iteration <= mMaxIterations + + if (rReport.multigrid_iteration) rReport.multigrid_iteration -= 1; + KRATOS_CATCH("") + } + + + [[nodiscard]] PMGStatusStream::Report + ExecuteConstraintLoop(typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs, + PMGStatusStream& rStream) + { + KRATOS_TRY + bool constraint_status = false; + PMGStatusStream::Report status_report { + /*grid_level=*/ 0ul, + /*multigrid_converged=*/ false, + /*multigrid_iteration=*/ 0ul, + /*maybe_multigrid_residual=*/ {}, + /*constraints_converged=*/ false, + /*constraint_iteration=*/ 0ul, + /*maybe_constraint_residual=*/ {} + }; + + typename TSparse::VectorType residual(rRhs.size()), + residual_update(rRhs.size()), + solution_update(rSolution.size()); + + // Compute the initial residual norm if it will be used. + typename TSparse::DataType initial_residual_norm = 1; + if (0 < mTolerance || 3 <= mVerbosity) { + initial_residual_norm = TSparse::TwoNorm(rRhs); + initial_residual_norm = initial_residual_norm ? initial_residual_norm : 1; + } + + // Outer loop for constraints. + do { + status_report.maybe_multigrid_residual.reset(); + + // Initialize the constraint assembler and update residuals. + mpConstraintAssembler->InitializeSolutionStep(rLhs, rSolution, rRhs); + TSparse::Copy(rRhs, residual); + BalancedProduct(rLhs, rSolution, residual, static_cast(-1)); + + // Get an update on the solution with respect to the current residual. + this->ExecuteMultigridLoop(rLhs, + rSolution, + rRhs, + solution_update, + residual, + initial_residual_norm, + rStream, + status_report); + + // Check for constraint convergence. + constraint_status = mpConstraintAssembler->FinalizeSolutionStep(rLhs, + rSolution, + rRhs, + status_report, + rStream); + if (constraint_status) { + rStream.Submit(status_report.Tag(2), mVerbosity); + status_report.maybe_constraint_residual.reset(); + break; + } else { + rStream.Submit(status_report.Tag(3), mVerbosity); + status_report.maybe_constraint_residual.reset(); + } + + ++status_report.constraint_iteration; + } while (true); + + return status_report; + KRATOS_CATCH("") + } + + + /// @brief Initialize the linear solver and solve the provided system. + bool Solve(typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs, + ModelPart& rModelPart) + { + // Prepare and initialize members. + KRATOS_TRY + if (mpInterface->GetLinearSolver().AdditionalPhysicalDataIsNeeded()) { + mpInterface->GetLinearSolver().ProvideAdditionalData(rLhs, + rSolution, + rRhs, + mpInterface->GetDofSet(), + rModelPart); + } + + if (mMaybeHierarchy.has_value()) { + std::visit([&rModelPart, &rLhs, &rSolution, &rRhs](auto& r_grid){ + r_grid.template Initialize( + rModelPart, + rLhs, + rSolution, + rRhs); + }, + mMaybeHierarchy.value()); + } // if mMaybeHierarchy + + std::optional status_stream = PMGStatusStream( + /*rStream=*/ std::cout, + /*rBuilderAndSolver=*/ *mpInterface, + /*rModelPart=*/ rModelPart, + /*rRootLhs=*/ rLhs, + /*rRootSolution=*/ rSolution, + /*rRootRhs=*/ rRhs, + /*UseAnsiColors=*/ true); + const auto status_report = this->ExecuteConstraintLoop( + rLhs, + rSolution, + rRhs, + status_stream.value()); + status_stream.reset(); + + // Clean up. + if (mMaybeHierarchy.has_value()) { + std::visit([&rModelPart, &rLhs, &rSolution, &rRhs](auto& r_grid){ + r_grid.template Finalize( + rModelPart, + rLhs, + rSolution, + rRhs); + }, + mMaybeHierarchy.value()); + } // if mMaybeHierarchy + + mpConstraintAssembler->Finalize(rLhs, + rSolution, + rRhs, + mpInterface->GetDofSet()); + + return status_report.multigrid_converged && status_report.constraints_converged; + KRATOS_CATCH("") + } + + + // --------------------------------------------------------- // + // Topology + // --------------------------------------------------------- // + + + template + static void CollectDoFs(const TProxy& rEntities, + const ProcessInfo& rProcessInfo, + typename Interface::TSchemeType& rScheme, + LockObject* pLockBegin, + [[maybe_unused]] LockObject* pLockEnd, + TIndexSet* pRowSetBegin, + [[maybe_unused]] TIndexSet* pRowSetEnd) + { + KRATOS_TRY + KRATOS_PROFILE_SCOPE(KRATOS_CODE_LOCATION); + using TLS = Element::EquationIdVectorType; + block_for_each(rEntities, + TLS(), + [&rScheme, &rProcessInfo, pLockBegin, pRowSetBegin](const auto& r_entity, TLS& r_tls){ + if (r_entity.GetEntity().IsActive()) { + rScheme.EquationId(r_entity.GetEntity(), r_tls, rProcessInfo); + for (const auto equation_id : r_tls) { + [[maybe_unused]] std::scoped_lock lock(pLockBegin[equation_id]); + auto& r_row_indices = pRowSetBegin[equation_id]; + r_row_indices.insert(r_tls.begin(), r_tls.end()); + } // for equation_id in r_tls + } // if r_entity.IsActive + }); + KRATOS_CATCH("") + } + + + void MakeLhsTopology(const typename Interface::TSchemeType::Pointer& rpScheme, + typename Interface::TSystemMatrixType& rLhs, + ModelPart& rModelPart) + { + KRATOS_PROFILE_SCOPE_MILLI(KRATOS_CODE_LOCATION); + KRATOS_TRY + + using IndexSet = std::unordered_set; + std::vector indices(mpInterface->GetEquationSystemSize()); + + { + std::vector mutexes(mpInterface->GetEquationSystemSize()); + + // Collect DoFs from elements. + Impl::CollectDoFs(MakeProxy(rModelPart), + rModelPart.GetProcessInfo(), + *rpScheme, + mutexes.data(), + mutexes.data() + mutexes.size(), + indices.data(), + indices.data() + indices.size()); + + // Collect DoFs from conditions. + Impl::CollectDoFs(MakeProxy(rModelPart), + rModelPart.GetProcessInfo(), + *rpScheme, + mutexes.data(), + mutexes.data() + mutexes.size(), + indices.data(), + indices.data() + indices.size()); + } + + // Compute and allocate LHS topology. + MakeSparseTopology(indices, + indices.size(), + rLhs, + /*EnsureDiagonal=*/true); + + // Construct the coarse hierarchy's topology. + if (mMaybeHierarchy.has_value()) { + std::visit([&rModelPart, &rLhs, this](auto& r_grid){ + r_grid.template MakeLhsTopology(rModelPart, + rLhs, + *mpConstraintAssembler, + mpInterface->GetDofSet()); + }, + mMaybeHierarchy.value()); + } // if mMaybeHierarchy + + KRATOS_CATCH("") + } + + + // --------------------------------------------------------- // + // Assembly + // --------------------------------------------------------- // + + + /// @brief Compute and assemble local contributions from elements and conditions into the global system. + /// @details This function body mainly handles looping over elements/conditions and its parallelization. + /// The actual local system calculation, as well as assembly, is deferred to @ref MapEntityContribution. + template + void Assemble(ModelPart& rModelPart, + typename Interface::TSchemeType& rScheme, + std::optional pMaybeLhs, + std::optional pMaybeRhs) + { + KRATOS_TRY + + // Sanity checks. + if constexpr (AssembleLHS) { + KRATOS_ERROR_IF_NOT(pMaybeLhs.has_value() && pMaybeLhs.value() != nullptr) + << "Requested an assembly of the left hand side, but no matrix is provided to assemble into."; + } + + if constexpr (AssembleRHS) { + KRATOS_ERROR_IF_NOT(pMaybeRhs.has_value() && pMaybeRhs.value() != nullptr) + << "Requested an assembly of the right hand side, but no vector is provided to assemble into."; + } + + // Assemble constraints. + // Constraint assembly MUST happen before assembling the unconstrained system, + // because constraints may have to propagate dirichlet conditions, and dirichlet + // conditions are partly imposed during element assembly. + mpConstraintAssembler->Assemble( + rModelPart.MasterSlaveConstraints(), + rModelPart.GetProcessInfo(), + mpInterface->GetDofSet(), + AssembleLHS, + AssembleRHS); + + // Function-wide variables. + const ProcessInfo& r_process_info = rModelPart.GetProcessInfo(); + typename Interface::TSystemMatrixType* pLhs = pMaybeLhs.has_value() ? pMaybeLhs.value() : nullptr; + typename Interface::TSystemVectorType* pRhs = pMaybeRhs.has_value() ? pMaybeRhs.value() : nullptr; + auto p_locks = std::make_unique>(AssembleLHS ? pLhs->size1() : 0ul); + + const int element_count = rModelPart.Elements().size(); + const int condition_count = rModelPart.Conditions().size(); + const auto it_element_begin = rModelPart.Elements().begin(); + const auto it_condition_begin = rModelPart.Conditions().begin(); + + // Collect contributions from constitutive entities. + #pragma omp parallel + { + // Thread-local variables. + Element::EquationIdVectorType equation_indices; + typename Interface::LocalSystemMatrixType lhs_contribution; + typename Interface::LocalSystemVectorType rhs_contribution; + + // Collect contributions from elements. + #pragma omp for schedule(guided, 512) nowait + for (int i_entity=0; i_entity( + *(it_element_begin + i_entity), + rScheme, + r_process_info, + equation_indices, + &lhs_contribution, + &rhs_contribution, + pLhs, + pRhs, + p_locks->data()); + } // pragma omp for + + // Collect contributions from conditions. + #pragma omp for schedule(guided, 512) + for (int i_entity=0; i_entity( + *(it_condition_begin + i_entity), + rScheme, + r_process_info, + equation_indices, + &lhs_contribution, + &rhs_contribution, + pLhs, + pRhs, + p_locks->data()); + } // pragma omp for + } // pragma omp parallel + + // Assemble coarse hierarchy. + if (this->mMaybeHierarchy.has_value()) { + std::visit([&rModelPart, &pMaybeLhs, &pMaybeRhs, this](auto& r_grid){ + r_grid.template Assemble( + rModelPart, + pMaybeLhs.has_value() ? pMaybeLhs.value() : nullptr, + pMaybeRhs.has_value() ? pMaybeRhs.value() : nullptr, + *this->mpConstraintAssembler, + mpInterface->GetDofSet()); + }, + this->mMaybeHierarchy.value()); + } // if mMaybeHierarchy + + KRATOS_CATCH("") + } + +}; // class PMultigridBuilderAndSolver::Impl + + +// --------------------------------------------------------- // +// Constructors +// --------------------------------------------------------- // + + +template +PMultigridBuilderAndSolver::~PMultigridBuilderAndSolver() +{ + delete mpImpl; +} + + +template +PMultigridBuilderAndSolver::PMultigridBuilderAndSolver() + : Interface(), + mpImpl(new Impl(this)) +{ +} + + +template +PMultigridBuilderAndSolver::PMultigridBuilderAndSolver(const typename LinearSolverType::Pointer& pSolver, + Parameters Settings) + : Interface(pSolver), + mpImpl(new Impl(this)) +{ + this->AssignSettings(Settings); +} + + +template +typename BuilderAndSolver>::Pointer +PMultigridBuilderAndSolver::Create(typename LinearSolverType::Pointer pSolver, + Parameters Settings) const +{ + KRATOS_TRY + return typename Interface::Pointer(new PMultigridBuilderAndSolver(pSolver, Settings)); + KRATOS_CATCH("") +} + + +// --------------------------------------------------------- // +// Allocation and Initialization +// --------------------------------------------------------- // + + +std::optional FindNodeIndex(std::size_t NodeId, + Kratos::ModelPart::NodesContainerType::const_iterator itNodeBegin, + Kratos::ModelPart::NodesContainerType::const_iterator itNodeEnd) noexcept +{ + const auto it = std::lower_bound(itNodeBegin, + itNodeEnd, + NodeId, + [](const Kratos::Node& r_node, std::size_t target_id){ + return r_node.Id() < target_id; + }); + if (it == itNodeEnd || it->Id() != NodeId) { + return {}; + } else { + return std::distance(itNodeBegin, it); + } +} + + +template +void PMultigridBuilderAndSolver::SetUpDofSet(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart) +{ + KRATOS_TRY; + + using DofsVectorType = ModelPart::DofsVectorType; + this->GetDofSet().clear(); + const auto& r_process_info = rModelPart.GetProcessInfo(); + + std::vector> hanging_nodes(rModelPart.Nodes().size()); + std::fill(hanging_nodes.begin(), + hanging_nodes.end(), + static_cast(1)); + + // Fill the global DOF set array + #pragma omp parallel + { + DofsVectorType tls_dofs, tls_constraint_dofs; + + // We create the temporal set in current thread and we reserve some space on it + std::unordered_set dofs_tmp_set; + dofs_tmp_set.reserve(20000); + + // Add the DOFs from the model part elements + const auto& r_elements_array = rModelPart.Elements(); + const int number_of_elements = static_cast(r_elements_array.size()); + + #pragma omp for schedule(guided, 512) nowait + for (int i_entity=0; i_entity(rModelPart.MasterSlaveConstraints().size()); + + #pragma omp for schedule(guided, 512) nowait + for (int i_entity=0; i_entityGetDofList(tls_dofs, tls_constraint_dofs, r_process_info); + + // Different behavior depending on whether MasterSlaveConstraint or MultifreedomConstraint. + if (tls_dofs.empty() && !tls_constraint_dofs.empty()) { + // Multifreedom constraint. + dofs_tmp_set.insert(tls_constraint_dofs.begin(), tls_constraint_dofs.end()); + } else { + // Master-slave constraint. + dofs_tmp_set.insert(tls_dofs.begin(), tls_dofs.end()); + dofs_tmp_set.insert(tls_constraint_dofs.begin(), tls_constraint_dofs.end()); + } + } // for i_entity in range(number_of_constraints) + + // Merge all the sets in one thread + #pragma omp critical + { + this->GetDofSet().insert(dofs_tmp_set.begin(), dofs_tmp_set.end()); + } + } // pragma omp parallel + + // Make sure that conditions act exclusively on collected DoFs. + #ifdef KRATOS_DEBUG + block_for_each(rModelPart.Conditions().begin(), + rModelPart.Conditions().end(), + DofsVectorType(), + [&r_process_info, this](const Condition& r_condition, auto& r_tls_dofs){ + if (r_condition.IsActive()) { + r_condition.GetDofList(r_tls_dofs, r_process_info); + for (Dof* p_dof : r_tls_dofs) { + KRATOS_ERROR_IF_NOT(this->GetDofSet().count(*p_dof)) + << "condition " << r_condition.Id() << " " + << "acts on " << p_dof->GetVariable().Name() << " " + << "of node " << p_dof->Id() << ", " + << "which is not part of any element or constraint."; + } // for p_dof in r_tls_dofs + } // if r_condition.IsActive() + }); + #endif + + #ifdef KRATOS_DEBUG + // If reactions are to be calculated, we check if all the dofs have reactions defined + // This is to be done only in debug mode + for (const auto& r_dof : this->GetDofSet()) { + KRATOS_ERROR_IF_NOT(r_dof.HasReaction()) + << "Reaction variable not set for " + << "DoF " << r_dof << " " + << "of node "<< r_dof.Id(); + } + #endif + + KRATOS_CATCH(""); +} + + +template +void PMultigridBuilderAndSolver::SetUpSystem(ModelPart& rModelPart) +{ + KRATOS_TRY + this->mEquationSystemSize = this->GetDofSet().size(); + + // Set equation indices of DoFs. + IndexPartition(this->GetDofSet().size()).for_each([&, this](std::size_t Index){ + (this->GetDofSet().begin() + Index)->SetEquationId(Index); + }); + KRATOS_CATCH("") +} + + +template +void PMultigridBuilderAndSolver::ResizeAndInitializeVectors(typename Interface::TSchemeType::Pointer pScheme, + typename Interface::TSystemMatrixPointerType& rpLhs, + typename Interface::TSystemVectorPointerType& rpSolution, + typename Interface::TSystemVectorPointerType& rpRhs, + ModelPart& rModelPart) +{ + KRATOS_TRY + + // Construct empty containers if necessary. + if (!rpLhs) + rpLhs.reset(new typename Interface::TSystemMatrixType); + TSparse::SetToZero(*rpLhs); + + if (!rpSolution) + rpSolution.reset(new typename Interface::TSystemVectorType); + if (rpSolution->size() != this->mEquationSystemSize) + rpSolution->resize(this->mEquationSystemSize, false); + TSparse::SetToZero(*rpSolution); + + if (!rpRhs) + rpRhs.reset(new typename Interface::TSystemVectorType); + if (rpRhs->size() != this->mEquationSystemSize) + rpRhs->resize(this->mEquationSystemSize, false); + TSparse::SetToZero(*rpRhs); + + // Construct LHS topology if necessary or requested. + if (rpLhs->size1() == 0 || this->GetReshapeMatrixFlag() == true) { + mpImpl->MakeLhsTopology(pScheme, *rpLhs, rModelPart); + + // Make constraint topology. + mpImpl->mpConstraintAssembler->Allocate(rModelPart.MasterSlaveConstraints(), + rModelPart.GetProcessInfo(), + *rpLhs, + *rpSolution, + *rpRhs, + this->GetDofSet()); + } else { + if (rpLhs->size1() != this->mEquationSystemSize || rpLhs->size2() != this->mEquationSystemSize) { + KRATOS_ERROR <<"The equation system size has changed during the simulation. This is not permitted."<SetDofSetIsInitializedFlag(true); + + KRATOS_CATCH("") +} + + +// --------------------------------------------------------- // +// Hooks +// --------------------------------------------------------- // + + +template +void PMultigridBuilderAndSolver::InitializeSolutionStep(ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) +{ + mpImpl->mpMaybeModelPart = &rModelPart; +} + + +template +void PMultigridBuilderAndSolver::FinalizeSolutionStep(ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) +{ + mpImpl->mpMaybeModelPart.reset(); +} + + +// --------------------------------------------------------- // +// Assembly +// --------------------------------------------------------- // + + +template +void PMultigridBuilderAndSolver::Build(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rRhs) +{ + KRATOS_PROFILE_SCOPE_MILLI(KRATOS_CODE_LOCATION); + KRATOS_ERROR_IF(!pScheme) << "missing scheme" << std::endl; + KRATOS_TRY + mpImpl->mpMaybeModelPart = &rModelPart; + mpImpl->template Assemble(rModelPart, + *pScheme, + &rLhs, + &rRhs); + KRATOS_CATCH("") +} + + +template +void PMultigridBuilderAndSolver::BuildLHS(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs) +{ + KRATOS_PROFILE_SCOPE_MILLI(KRATOS_CODE_LOCATION); + KRATOS_ERROR_IF(!pScheme) << "missing scheme"; + KRATOS_TRY + mpImpl->template Assemble(rModelPart, + *pScheme, + &rLhs, + nullptr); + KRATOS_CATCH("") +} + + +template +void PMultigridBuilderAndSolver::BuildRHS(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemVectorType& rRhs) +{ + KRATOS_PROFILE_SCOPE_MILLI(KRATOS_CODE_LOCATION); + KRATOS_ERROR_IF(!pScheme) << "missing scheme"; + KRATOS_TRY + mpImpl->template Assemble(rModelPart, + *pScheme, + nullptr, + &rRhs); + block_for_each(this->GetDofSet(), [&rRhs](Dof& rDof){ + if (rDof.IsFixed()) + rRhs[rDof.EquationId()] = 0.0; + }); + KRATOS_CATCH("") +} + + +// --------------------------------------------------------- // +// Constraint Imposition +// --------------------------------------------------------- // + + +template +void PMultigridBuilderAndSolver::ApplyDirichletConditions(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + [[maybe_unused]] typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) +{ + KRATOS_TRY + KRATOS_PROFILE_SCOPE(KRATOS_CODE_LOCATION); + + // Early exit on an empty system. + if (this->GetDofSet().empty()) return; + + const auto it_dof_set_begin = this->GetDofSet().begin(); + const auto it_dof_set_end = this->GetDofSet().begin() + this->GetDofSet().size(); + + mpImpl->mpDiagonalScaling->template Cache(rLhs); + const auto diagonal_scale = mpImpl->mpDiagonalScaling->Evaluate(); + Kratos::ApplyDirichletConditions(rLhs, + rRhs, + it_dof_set_begin, + it_dof_set_end, + diagonal_scale); + if (mpImpl->mMaybeHierarchy.has_value()) { + std::visit([this](auto& r_grid){ + r_grid.ApplyDirichletConditions(this->GetDofSet().begin(), + this->GetDofSet().end()); + }, + mpImpl->mMaybeHierarchy.value()); + } // if mMaybeHierarchy + + KRATOS_CATCH("") +} + + +template +void PMultigridBuilderAndSolver::ApplyConstraints(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rRhs) +{ + KRATOS_TRY + KRATOS_PROFILE_SCOPE_MILLI(KRATOS_CODE_LOCATION); + mpImpl->mpConstraintAssembler->Initialize(rLhs, + rRhs, + this->GetDofSet().begin(), + this->GetDofSet().end()); + if (mpImpl->mMaybeHierarchy.has_value()) { + std::visit([](auto& r_hierarchy){r_hierarchy.ApplyConstraints();}, + mpImpl->mMaybeHierarchy.value()); + } + KRATOS_CATCH("") +} + + +// --------------------------------------------------------- // +// Compound Assembly and Solution +// --------------------------------------------------------- // + + +template +void PMultigridBuilderAndSolver::BuildAndSolve(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) +{ + KRATOS_TRY + + // Assemble unconstrained system. + Build(pScheme, rModelPart, rLhs, rRhs); + + // Apply multifreedom constraints. + ApplyConstraints(pScheme, rModelPart, rLhs, rRhs); + + // Apply Dirichlet conditions. + ApplyDirichletConditions(pScheme, rModelPart, rLhs, rSolution, rRhs); + + this->SystemSolve(rLhs, rSolution, rRhs); + KRATOS_CATCH("") +} + + +template +void PMultigridBuilderAndSolver::SystemSolve(typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) +{ + KRATOS_TRY + // Solve constrained assembled system. + BuiltinTimer timer; + if (!mpImpl->Solve(rLhs, rSolution, rRhs, *mpImpl->mpMaybeModelPart.value()) && 1 <= mpImpl->mVerbosity) { + std::cerr << this->Info() << "Grid 0: failed to solve the assembled system\n"; + } + KRATOS_INFO_IF(this->Info(), 2 <= mpImpl->mVerbosity) + << ": Linear solution took " << timer << "\n"; + KRATOS_CATCH("") +} + + +// --------------------------------------------------------- // +// Postprocessing +// --------------------------------------------------------- // + + + +template +void PMultigridBuilderAndSolver::CalculateReactions(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) +{ + KRATOS_PROFILE_SCOPE_MILLI(KRATOS_CODE_LOCATION); + + KRATOS_TRY + BuiltinTimer timer; + + TSparse::SetToZero(rRhs); + mpImpl->template Assemble(rModelPart, + *pScheme, + nullptr, + &rRhs); + block_for_each(this->GetDofSet(), [&rRhs](Dof& rDof){ + rDof.GetSolutionStepReactionValue() = -rRhs[rDof.EquationId()]; + }); + + KRATOS_INFO_IF(this->Info(), 2 <= mpImpl->mVerbosity) + << "Computing reactions took " << timer << "\n"; + KRATOS_CATCH("") +} + + +// --------------------------------------------------------- // +// Misc +// --------------------------------------------------------- // + + +template +void PMultigridBuilderAndSolver::AssignSettings(const Parameters Settings) +{ + // Mutable parameters. + Parameters settings = Settings; + + KRATOS_TRY + // Parse diagonal scaling and validate other settings. + Parameters default_parameters = this->GetDefaultParameters(); + Parameters default_diagonal_scaling = default_parameters["diagonal_scaling"].Clone(); + default_parameters.RemoveValue("diagonal_scaling"); + std::optional maybe_diagonal_scaling = settings.Has("diagonal_scaling") ? settings["diagonal_scaling"].Clone() : std::optional(); + if (maybe_diagonal_scaling.has_value()) settings.RemoveValue("diagonal_scaling"); + settings.ValidateAndAssignDefaults(default_parameters); + settings.AddValue("diagonal_scaling", maybe_diagonal_scaling.has_value() ? maybe_diagonal_scaling.value() : default_diagonal_scaling); + mpImpl->mpDiagonalScaling = std::make_unique(settings["diagonal_scaling"]); + KRATOS_CATCH("") + + KRATOS_TRY + Interface::AssignSettings(settings); + KRATOS_CATCH("") + + // Set the relative residual tolerance and max W cycles. + mpImpl->mMaxIterations = settings["max_iterations"].Get(); + mpImpl->mTolerance = settings["tolerance"].Get(); + + // Set multifreedom constraint imposition strategy. + mpImpl->mpConstraintAssembler = ConstraintAssemblerFactory(settings["constraint_imposition_settings"], + "Grid 0 constraints"); + + // Construct smoother, solver and coarse hierarchy. + Parameters smoother_settings = settings["smoother_settings"]; + Parameters coarse_hierarchy_settings = settings["coarse_hierarchy_settings"]; + Parameters leaf_solver_settings = settings["linear_solver_settings"]; + + if (coarse_hierarchy_settings["max_depth"].Get()) { + KRATOS_TRY + KRATOS_ERROR_IF_NOT(smoother_settings.Has("solver_type")); + const std::string solver_name = smoother_settings["solver_type"].Get(); + using SolverFactoryRegistry = KratosComponents>; + KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name)) + << "\"" << solver_name << "\" is not a valid linear solver name in the registry. " + << "Make sure you imported the application it is defined in and that the spelling is correct."; + const auto& r_factory = SolverFactoryRegistry::Get(solver_name); + this->mpLinearSystemSolver = r_factory.Create(smoother_settings); + + // Construct the coarse hierarchy. + const std::string coarse_build_precision = coarse_hierarchy_settings["precision"].Get(); + if (coarse_build_precision == "double") { + using CoarseSparseSpace = TUblasSparseSpace; + using CoarseDenseSpace = TUblasDenseSpace; + using GridType = PGrid; + + coarse_hierarchy_settings.ValidateAndAssignDefaults(GridType().GetDefaultParameters()); + if (coarse_hierarchy_settings["max_depth"].Get()) { + mpImpl->mMaybeHierarchy = GridType(coarse_hierarchy_settings, + smoother_settings, + leaf_solver_settings, + settings["diagonal_scaling"]); + } + } /* if coarse_build_precision == "double" */ else if (coarse_build_precision == "single") { + using CoarseSparseSpace = TUblasSparseSpace; + using CoarseDenseSpace = TUblasDenseSpace; + using GridType = PGrid; + + coarse_hierarchy_settings.ValidateAndAssignDefaults(GridType().GetDefaultParameters()); + mpImpl->mMaybeHierarchy = GridType(coarse_hierarchy_settings, + smoother_settings, + leaf_solver_settings, + settings["diagonal_scaling"]); + } /* elif coarse_build_precision == "single" */ else { + KRATOS_ERROR << "unsupported coarse precision: \"" << coarse_build_precision << "\". " + << "Options are:\n" + << "\t\"single\"" + << "\t\"double\""; + } + KRATOS_CATCH("") + } /* if coarse_hierarchy_settings["max_depth"].Get() */ else { + KRATOS_TRY + KRATOS_ERROR_IF_NOT(leaf_solver_settings.Has("solver_type")); + const std::string solver_name = leaf_solver_settings["solver_type"].Get(); + using SolverFactoryRegistry = KratosComponents>; + KRATOS_ERROR_IF_NOT(SolverFactoryRegistry::Has(solver_name)) + << "\"" << solver_name << "\" is not a valid linear solver name in the registry. " + << "Make sure you imported the application it is defined in and that the spelling is correct."; + const auto& r_factory = SolverFactoryRegistry::Get(solver_name); + this->mpLinearSystemSolver = r_factory.Create(leaf_solver_settings); + KRATOS_CATCH("") + } + + // Other settings. + KRATOS_TRY + mpImpl->mVerbosity = settings["verbosity"].Get(); + KRATOS_CATCH("") +} + + +template +void PMultigridBuilderAndSolver::Clear() +{ + Interface::Clear(); + mpImpl->mpConstraintAssembler->Clear(); + + if (mpImpl->mMaybeHierarchy.has_value()) { + std::visit([](auto& r_grid) {r_grid.Clear();}, + mpImpl->mMaybeHierarchy.value()); + } // if mMaybeHierarchy + + this->SetDofSetIsInitializedFlag(false); +} + + +template +Parameters PMultigridBuilderAndSolver::GetDefaultParameters() const +{ + Parameters parameters = Parameters(R"({ + "name" : "p_multigrid", + "diagonal_scaling" : "max", + "max_iterations" : 1e2, + "tolerance" : 1e-8, + "verbosity" : 1, + "smoother_settings" : { + "solver_type" : "" + }, + "linear_solver_settings" : { + "solver_type" : "amgcl" + }, + "constraint_imposition_settings" : { + "method" : "master_slave_elimination" + }, + "coarse_hierarchy_settings" : {} + })"); + parameters.SetValue("coarse_hierarchy_settings", PGrid,TUblasDenseSpace>::GetDefaultParameters()); + parameters.RecursivelyAddMissingParameters(Interface::GetDefaultParameters()); + return parameters; +} + + +template +std::string PMultigridBuilderAndSolver::Info() const +{ + return "PMultigridBuilderAndSolver"; +} + + +template +std::size_t PMultigridBuilderAndSolver::GetEquationSystemSize() const noexcept +{ + return Interface::mEquationSystemSize; +} + +template +LinearSolver& PMultigridBuilderAndSolver::GetLinearSolver() noexcept +{ + return *Interface::mpLinearSystemSolver; +} + + +template +void PMultigridBuilderAndSolver::ProjectGrid(int GridLevel, + const typename TSparse::MatrixType& rRootLhs, + const typename TSparse::VectorType& rRootSolution, + const typename TSparse::VectorType& rRootRhs) +{ + KRATOS_TRY + + typename TSparse::VectorType projected; + + // Project residuals. + if (0 < GridLevel) { + KRATOS_ERROR_IF_NOT(mpImpl->mMaybeHierarchy.has_value()) + << "grid level " << GridLevel << " does not exist"; + const auto i_grid_level = GridLevel - 1; + + // Construct a flat vector of coarse grids. + std::variant< + std::vector,TUblasDenseSpace>*>, + std::vector,TUblasDenseSpace>*> + > coarse_grids; + + std::visit([&coarse_grids](const auto& r_coarse_grid){ + using GridType = std::remove_cv_t>; + coarse_grids = std::vector(); + }, mpImpl->mMaybeHierarchy.value()); + + const auto vector_fill_visitor = [&coarse_grids](const auto& r_coarse_grid) { + std::visit([&r_coarse_grid](auto& r_vector){ + using GridType = std::remove_cv_t>; + using ValueType = typename std::remove_cv_t>::value_type; + if constexpr (std::is_same_v) { + const GridType* p_grid = &r_coarse_grid; + do { + r_vector.push_back(p_grid); + const auto p_maybe_child = p_grid->GetChild(); + p_grid = p_maybe_child.has_value() ? p_maybe_child.value() : nullptr; + } while (p_grid); + } + }, coarse_grids); + }; + + std::visit(vector_fill_visitor, mpImpl->mMaybeHierarchy.value()); + + // Initialize projected solution vector. + std::visit([&projected, i_grid_level](const auto& r_vector){ + KRATOS_ERROR_IF_NOT(static_cast(i_grid_level) < r_vector.size()) + << "grid level " << i_grid_level + 1 << " does not exist"; + const auto& r_solution = r_vector[i_grid_level]->GetSolution(); + projected.resize(r_solution.size(), false); + std::fill(projected.begin(), projected.end(), static_cast(0)); + }, coarse_grids); + + // Project solution to the root grid. + { + typename TSparse::VectorType tmp; + for (int i_grid=i_grid_level; 0<=i_grid; --i_grid) { + std::visit([&projected, &tmp, i_grid](const auto& r_vector) { + const auto& r_solution = r_vector[i_grid]->GetSolution(); + IndexPartition(r_solution.size()).for_each([&r_solution, &projected](std::size_t i_component){ + projected[i_component] += r_solution[i_component]; + }); + r_vector[i_grid]->template Prolong(projected, tmp); + }, coarse_grids); + tmp.swap(projected); + } + } + + } /*if 0 < GridLevel*/ else { + projected.resize(rRootRhs.size(), false); + TSparse::SetToZero(projected); + } + + // Add diff of the root grid. + TSparse::UnaliasedAdd(projected, static_cast(1), rRootSolution); + + // Compute residual. + auto residual = rRootRhs; + BalancedProduct(rRootLhs, + projected, + residual, + static_cast(-1)); + + // Update DoFs. + IndexPartition(projected.size()).for_each([this, &projected, &residual](std::size_t i_dof) { + auto& r_dof = *(this->GetDofSet().begin() + i_dof); + (this->GetDofSet().begin() + i_dof)->GetSolutionStepReactionValue() = residual[i_dof]; + r_dof.GetSolutionStepValue() += projected[i_dof]; + }); + + KRATOS_CATCH("") +} + + +// --------------------------------------------------------- // +// Template Instantiations +// --------------------------------------------------------- // + +template class KRATOS_API(KRATOS_CORE) PMultigridBuilderAndSolver, + TUblasDenseSpace>; + +template class KRATOS_API(KRATOS_CORE) PMultigridBuilderAndSolver, + TUblasDenseSpace>; + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp new file mode 100644 index 000000000000..55b3251ea75d --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp @@ -0,0 +1,446 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "solving_strategies/builder_and_solvers/builder_and_solver.h" // BuilderAndSolver +#include "includes/smart_pointers.h" // KRATOS_CLASS_POINTER_DEFINITION +#include "includes/model_part.h" // ModelPart +#include "includes/code_location.h" // KRATOS_CODE_LOCATION +#include "linear_solvers/linear_solver.h" // LinearSolver + +// STL Includes +#include // unique_ptr + + +namespace Kratos { + + + +/** @page P-Multigrid + * + * P-Multigrid refers to a multigrid method whose coarsening strategy is based on shape functions of high order elements (quadratic elements already qualify in this case). This implementation also supports models with multifreedom constraints. + * + * @section Overview + * + * The system's entry point is the @ref PMultigridBuilderAndSolver. As the name suggests, the main responsibilities of this class are + * - allocating the left hand side matrix, as well as the right hand side vector and the solution vector + * - assembling the left hand side matrix, right hand side vector + * - solving the resulting linear system of equations. + * + * Due to the current design of @ref Scheme and @ref BuilderAndSolver, there is a number of other tasks as well that either tie in the primary purpose or help in pre- and postprocessing: + * - allocation, assembly and imposiotion of @ref MasterSlaveConstraint "multifreedom constraints" + * - applying Dirichlet conditions + * - partial reassembly of the linear system + * - computing reactions. + * + * The other exposed family of classes from this system is @ref MultifreedomConstraint. Its main difference compared to + * @ref MasterSlaveConstraint (which it derives from) is that it only models the constraint equations and does not impose + * a partition of slave- and master @ref Dof "DoFs". + * @see LinearMultifreedomConstraint. + * @see LinkConstraint + * + * @note The rest of the classes in this system are meant for internal use and must not be exposed to the rest of Kratos. + * + * Although the primary purpose of this system is to exploit the structure arising from a model using high order elements, + * the multigrid feature can be completely disabled to use it as a standard @ref BuilderAndSolver. One reason for + * doing this would be taking advantage of different constraint imposition methods, such as + * @ref AugmentedLagrangeConstraintAssembler "augmented Lagrange". + * + * @section Coarse Hierarchy + * + * @note The current implementation only supports a two-grid method, since the selection of high order elements + * in Kratos is limited, and does not yet justify an arbitrary depth. That said, adding support for it + * should be possible without major interface changes, but would involve minor changes in + * @ref PMultigridBuilderAndSolver and @ref PGrid, as well as major changes in @ref MakePRestrictionOperator. + * + * The root grid (i.e.: the finest level) is stored in and represented by @ref PMultigridBuilderAndSolver, while + * coarse grids are represented by @ref PGrid in a linked list. The reason for this difference is the additional + * set of responsibilities of the root grid, namely the allocation and assembly of the finest level. Coarse grids + * do not perform assembly, but construct restriction operators that they then apply on the parent grid to compute + * their own system. + * + * Another important distinction is that the coarse grids can have floating point types different than the root grid. + * This can be useful when the user has access to accelerator hardware (i.e.: GPUs). Coarse grids need not solve their + * own problems with high precision so they might as well use single precision floating point numbers to save VRAM, + * for example. + * + * @section Constraints + * + * Constraint assembly and imposition is extracted through @ref ConstraintAssembler "a dedicated interface" that + * currently supports @ref MasterSlaveConstraintAssembler "master-slave elimination", + * @ref AugmentedLagrangeConstraintAssembler "augmented Lagrange", and @ref NoOpConstraintAssembler "a dummy" + * for debugging. + * + * @note If the multigrid feature is enabled, the current implementation only supports augmented Lagrange imposition. + * + * @section Linear Solvers + * + * Unlike other @ref BuilderAndSolver "BuilderAndSolvers", @ref PMultigridBuilderAndSolver does not use the linear + * solver provided to it from the python layer. The reason is that it must construct a solver for each grid level + * separately. These fall into two categories + * - smoothers for the finer grids (including the root grid) + * - linear solver for the coarsest grid (usually an AMG solver). + * + * Instead of passing linear solver instances, the user must provide two sets of parameters for the two different + * solver categories, after which the grids take care of constructing their own instances. + */ + + +///@name Kratos Classes +///@{ + +/** @brief Class used as a standard @ref BuilderAndSolver with a built-in optional preconditioner for high-order elements. See @ref P-Multigrid. + * + * @details + * @section Structure + * PMultigridBuilderAndSolver consists of 3 main components: + * - The top-level grid, which is responsible for assembling the main linear system + * of equations from @ref Element "elements" and @ref Condition "conditions". + * - The top-level @ref ConstraintAssembler "constraint assembler" responsible for assembling + * and imposing constraint equations defined by @ref MasterSlaveConstraint "multifreedom constraints". + * - A hierarchy of coarser grid levels, each with its own @ref LinearSolver "linear solver" and + * @ref ConstraintAssembler "constraint assembler". + * + * @copydoc PMultigridBuilderAndSolver::GetDefaultParameters + * + * @subsection Top-Level-Grid Top-Level Grid + * Like any other grid level, the fine grid has its own linear solver and constraint assembler, but + * they can be configured separately from coarse grid levels. + * Settings: + * @code + * { + * "name" : "p_multigrid", + * "diagonal_scaling" : "max", + * "max_iterations" : 1e2, + * "tolerance" : 1e-8, + * "verbosity" : 1, + * ... + * } + * @endcode + * - @p "name" Name of this builder-and-solver to refer to in @ref PythonSolver solvers. + * - @p "diagonal_scaling" Expression to scale diagonal entries of the LHS that are + * related to @ref Dof "DoFs" constrained by Dirichlet conditions. + * It can either be a numeric literal, or a string that defines + * the scaling factor as a function of: + * - @p "max" Infinity norm (maximum absolute value) of the LHS' main diagonal. + * - @p "norm" 2-norm of the LHS' main diagonal. + * - @p "max_iterations" Max number of multigrid iterations to carry out while attempting + * to reduce the residual to the requested value. + * - @p "tolerance" Target relative residual norm to achieve. + * - @p "verbosity" Level of verbosity. Every level includes lower verbosity levels and + * adds new events to report: + * - @p 0 No messages, not even in case of failure. + * - @p 1 Warnings and failure messages. + * - @p 2 Aggregated status reports. + * - @p 3 Per-iteration status reports. + * - @p 4 Output system matrices and vectors. + * - @p 5 Write a VTU output with the solution and residuals at each iteration. + * + * @subsection Top-Level-Constraints Top-Level Constraint Assembler + * The fine-level constraint assembler is responsible for building and imposing + * constraint equations defined by @ref MasterSlaveConstraint "multifreedom constraints" of + * the input @ref ModelPart "model part". + * Settings: + * @code + * { + * ... + * "constraint_imposition_settings" { + * "method" "master_slave_elimination" + * } + * ... + * } + * @endcode + * - @p "constraint_imposition_settings" Settings of the top-level constraint assembler. + * Refer to @ref ConstraintAssemblerFactory for more + * information. + * + * @subsection Coarse-Hierarchy Coarse Hierarchy + * Coarse grid levels have their own linear solvers as well as constraint assemblers. However, + * all of them share the same configuration (except the coarses grid's linear solver, which is + * defined by @p "linear_solver_settings" ). Coarse grids and their constraint assemblers do + * not actually construct their systems by assembly, but by applying their restriction operators + * to their parent grid's matrices/vectors. Refer to @ref PGrid for more information. + * Default settings: + * @code + * { + * ... + * "coarse_hierarchy_settings" {} + * ... + * } + * @endcode + * - @p "coarse_hierarchy_settings" Coarse grid hierarchy configuration. + * + * @copydoc PGrid + * + * @subsection Solvers-Smoother Solvers and Smoothers + * PMultigridBuilderAndSolver uses two linear solvers in general. A smoother defined by + * @p "smoother_settings" for all grid levels except the coarses, and a proper linear solver + * on the coarsest grid, defined by @p "linear_solver_settings". + * Note that if @p "coarse_hierarchy_settings.max_depth" is @p 0, the coarsest grid level is + * the top-level one, in which case the smoother is not used. Both solvers are constructed by + * @ref LinearSolverFactory, so check there for more information. + * Settings: + * @code + * { + * ... + * "smoother_settings" { + * "solver_type" : "" + * }, + * "linear_solver_settings" : { + * "solver_type" : "amgcl" + * } + * } + * @endcode + */ +template +class KRATOS_API(KRATOS_CORE) PMultigridBuilderAndSolver + : public BuilderAndSolver> +{ +private: + using LinearSolverType = LinearSolver; + + using Interface = BuilderAndSolver; + +public: + /// @internal + KRATOS_CLASS_POINTER_DEFINITION(PMultigridBuilderAndSolver); + + /// @details Required by PIMPL. + /// @internal + ~PMultigridBuilderAndSolver() override; + + PMultigridBuilderAndSolver(); + + /// @brief Construct from a linear solver and parameters. + /// @note The provided linear solver is not used. Check the input settings + /// on how to configure linear solvers for this class. + /// @see PMultigridBuilderAndSolver::GetDefaultParameters + PMultigridBuilderAndSolver(const typename LinearSolverType::Pointer& pSolver, + Parameters Settings); + + PMultigridBuilderAndSolver(PMultigridBuilderAndSolver&& rOther) noexcept = default; + + PMultigridBuilderAndSolver& operator=(PMultigridBuilderAndSolver&& rRhs) noexcept = default; + + /// @brief Construct from a linear solver and parameters. + /// @see PMultigridBuilderAndSolver::GetDefaultParameters + typename BuilderAndSolver::Pointer + Create(typename LinearSolverType::Pointer pNewLinearSystemSolver, + Parameters ThisParameters) const override; + + /// @name Allocation and Initialization + /// @{ + + /// @copydoc BuilderAndSolver::SetUpDofSet + void SetUpDofSet(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart) override; + + /// @brief Assign equation indices to @ref Dof "degrees of freedom". + void SetUpSystem(ModelPart& rModelPart) override; + + /// @copydoc BuilderAndSolver::ResizeAndInitializeVectors + void ResizeAndInitializeVectors(typename Interface::TSchemeType::Pointer pScheme, + typename Interface::TSystemMatrixPointerType& rpLhs, + typename Interface::TSystemVectorPointerType& rpSolution, + typename Interface::TSystemVectorPointerType& rpRhs, + ModelPart& rModelPart) override; + + /// @} + /// @name Hooks + /// @{ + + void InitializeSolutionStep(ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) override; + + void FinalizeSolutionStep(ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) override; + + /// @} + /// @name Assembly + /// @{ + + /// @copydoc BuilderAndSolver::BuildLHS + void BuildLHS(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs) override; + + /// @copydoc BuilderAndSolver::BuildRHS + void BuildRHS(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemVectorType& rRhs) override; + + /// @copydoc BuilderAndSolver::Build + void Build(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rRhs) override; + + /// @} + /// @name Constraint Imposition + /// @{ + + /// @copydoc BuilderAndSolver::ApplyDirichletConditions + void ApplyDirichletConditions(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) override; + + /// @copydoc BuilderAndSolver::ApplyConstraints + void ApplyConstraints(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rRhs) override; + + /// @} + /// @name Solution + /// @{ + + /// @copydoc BuilderAndSolver::BuildAndSolve + void BuildAndSolve(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& A, + typename Interface::TSystemVectorType& Dx, + typename Interface::TSystemVectorType& b) override; + + /// @copydoc BuilderAndSolver::SystemSolve(typename Interface::TSystemMatrixType&, typename Interface::TSystemVectorType&, typename Interface::TSystemVectorType&) + void SystemSolve(typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) override; + + /// @} + /// @name Postprocessing + /// @{ + + /// @copydoc BuilderAndSolver::CalculateReactions + void CalculateReactions(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) override; + + /// @} + /// @name Misc + /// @{ + + /// @copydoc BuilderAndSolver::Clear + void Clear() override; + + /// @brief Default configuration. + /// @details @code + /// { + /// "name" : "p_multigrid", + /// "diagonal_scaling" : "max", + /// "max_iterations" : 1e2, + /// "tolerance" : 1e-8, + /// "verbosity" : 1, + /// "constraint_imposition_settings" : { + /// "method" : "master_slave_elimination" + /// }, + /// "coarse_hierarchy_settings" : {} + /// "smoother_settings" : { + /// "solver_type" : "" + /// }, + /// "linear_solver_settings" : { + /// "solver_type" : "amgcl" + /// }, + /// } + /// @endcode + Parameters GetDefaultParameters() const override; + + /// @copydoc BuilderAndSolver::Info + std::string Info() const override; + + /// @} + /// @name Debug and Analysis + /// @{ + + /// @brief Project state variables and residuals from a coarse grid to the root grid. + /// @param GridLevel Grid to project. + /// @param rRootLhs LHS matrix of the root grid. + /// @param rRootSolution Current state vector of the root grid. + /// @param rRootRhs RHS vector of the root grid. + /// @details State variables are stored in the corresponding @ref Dof "DoFs"' values, + /// while residuals are written to reactions. + /// @warning This function changes DoFs' values and reactions, so they should be copied + /// before invoking this function and restored afterwards. + void ProjectGrid(int GridLevel, + const typename TSparseSpace::MatrixType& rRootLhs, + const typename TSparseSpace::VectorType& rRootSolution, + const typename TSparseSpace::VectorType& rRootRhs); + + /// @} + /// @name Unsupported Interface + /// @{ + + /// @warning Not implemented. + void BuildRHSAndSolve(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not implemented\n";} + + /// @warning Not implemented. + void ApplyRHSConstraints(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemVectorType& rRhs) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not implemented\n";} + + /// @warning Not implemented. + void BuildLHS_CompleteOnFreeRows(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& A) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not implemented\n";} + + /// @warning Not implemented. + void BuildAndSolveLinearizedOnPreviousIteration(typename Interface::TSchemeType::Pointer pScheme, + ModelPart& rModelPart, + typename Interface::TSystemMatrixType& rLhs, + typename Interface::TSystemVectorType& rSolution, + typename Interface::TSystemVectorType& rRhs, + const bool MoveMesh) override + {KRATOS_ERROR << KRATOS_CODE_LOCATION.CleanFunctionName() << " is not implemented\n";} + + /// @} + +protected: + void AssignSettings(const Parameters Settings) override; + +private: + std::size_t GetEquationSystemSize() const noexcept; + + LinearSolverType& GetLinearSolver() noexcept; + + PMultigridBuilderAndSolver(const PMultigridBuilderAndSolver& rOther) = delete; + + PMultigridBuilderAndSolver& operator=(const PMultigridBuilderAndSolver& rRhs) = delete; + + /// @details MSVC desperately wants to know about the destructor of @p Impl when trying to + /// use PIMPL with an @p std::unique_ptr, up to the point that it ICE-s if I try + /// dodging it with a lambda. Hopefully a raw pointer will get the message + /// across its thick skull that it does not need to know what the destructor + /// does if it only ever sees a bloody pointer. + /// @todo Change @p Impl* to @p std::unique_ptr when (if) MSVC sorts its shit out. + struct Impl; + Impl* mpImpl; +}; // class PMultigridBuilderAndSolver + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_utilities.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_utilities.hpp new file mode 100644 index 000000000000..decaf2db41dd --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_utilities.hpp @@ -0,0 +1,631 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "geometries/geometry_data.h" // GeometryData +#include "geometries/geometry.h" // Geometry +#include "spaces/ublas_space.h" // TUblasSparseSpace +#include "includes/lock_object.h" // LockObject + +// System includes +#include // std::tuple +#include // std::vector +#include // std::unordered_map +#include // std::numeric_limits +#include // std::uint8_t +#include // std::remove_if + + +namespace Kratos { + + +namespace detail { +template +struct UnionReduction +{ + using value_type = TTable; + + using return_type = TTable; + + return_type mValue; + + return_type&& GetValue() + { + return std::move(mValue); + } + + void LocalReduce(const value_type& rValue) + { + mValue.insert(rValue.begin(), rValue.end()); + } + + void ThreadSafeReduce(UnionReduction& rOther) + { + KRATOS_CRITICAL_SECTION + { + for (auto& r_item : rOther) { + mValue.insert(std::move(r_item)); + } + } + } +}; // class UnionReduction +} // namespace detail + + +/** @brief Compute the p-restriction operator of a single geometry. + * @tparam OrderReduction Defines how many polynomial order to reduce the incoming + * geometry, at the maximum. If 0, the function returns a correctly sized + * identity matrix. Reduction to linear geometries will always happen if + * @p OrderReduction is greater or equal than the geometry's order. + * @tparam TValue Value type of entries. + * @tparam TIndex + * @tparam TNode Node type (usually @ref Node). + * @tparam TOutputIterator Output iterator of triplets. Nonzero entries in the restriction operator + * are provided as @p std::tuple triplets of {coarse node index, fine node index, value}. + * @param rGeometry Geometry to construct a restriction operator for. + * @param itOutput Output iterator where nonzero entries in the restriction operator are written to. Entries are + * represented as {coarse node index, fine node index, value} triplets (COO matrix). + * @details This function computes the (local) restriction operator @f$ R @f$ for a geometry that + * transforms the input (higher order @f$ q @f$) geometry's shape functions to an equivalent + * lower order (@f$ p @f$). Due to the linear relationship between shape functions and + * degrees-of-freedom, this mapping holds for DoFs as well: + * @f[ + * u^p_i = R_{ij} u^q_j + * @f] + * The reduction of the input geometry's order is defined by `OrderReduction` and is limited + * to reducing to a linear geometry. Any `OrderReduction` greater than that will only reduce + * to linear geometries as well. + */ +template +void MakePRestrictionOperator(const Geometry& rGeometry, + TOutputIterator itOutput) +{ + using G = GeometryData::KratosGeometryType; + const G geometry_type = rGeometry.GetGeometryType(); + + using Triplet = std::tuple; + + switch(geometry_type) { + /// - @ref GeometryData::KratosGeometryType::Kratos_Point2D "Kratos_Point2D" and @ref GeometryData::KratosGeometryType::Kratos_Point3D "Kratos_Point3D" + /// can't be higher order, so they just map back to themselves. + /// \f[\begin{bmatrix} + /// 1 + /// \end{bmatrix}\f] + case G::Kratos_Point2D: + case G::Kratos_Point3D: { + *itOutput = Triplet(0, 0, 1); + break; + } // Kratos_Point2D || Kratos_Point3D + + /** - @ref GeometryData::KratosGeometryType::Kratos_Line2D2 "Kratos_Line2D2" and @ref GeometryData::KratosGeometryType::Kratos_Line3D2 "Kratos_Line3D2" + * linear segment maps to itself. + * \f[\begin{bmatrix} + * 1 & \\ + * & 1 + * \end{bmatrix}\f] + */ + case G::Kratos_Line3D2: + case G::Kratos_Line2D2: { + *itOutput++ = Triplet(0, 0, 1); + *itOutput = Triplet(1, 1, 1); + break; + } // Kratos_Line2D2 || Kratos_Line3D2 + + /** - @ref GeometryData::KratosGeometryType::Kratos_Line2D3 "Kratos_Line2D3" and @ref GeometryData::KratosGeometryType::Kratos_Line3D3 "Kratos_Line3D3" + * quadratic line segments reduces to linear line segments + * \ref GeometryData::KratosGeometryType::Kratos_Line2D2 "Kratos_Line2D2" and @ref GeometryData::KratosGeometryType::Kratos_Line3D2 "Kratos_Line3D2". + * \f[\begin{bmatrix} + * 1 & & \frac{1}{2} \\ + * & 1 & \frac{1}{2} + * \end{bmatrix}\f] + */ + case G::Kratos_Line3D3: + case G::Kratos_Line2D3: { + if constexpr (OrderReduction == 0) { + *itOutput++ = Triplet(0, 0, 1.0); + *itOutput++ = Triplet(1, 1, 1.0); + *itOutput = Triplet(2, 2, 1.0); + } /*if OrderReduction == 0*/ else { + *itOutput++ = Triplet(0, 0, 1.0); + *itOutput++ = Triplet(0, 2, 0.5); + *itOutput++ = Triplet(1, 1, 1.0); + *itOutput = Triplet(1, 2, 0.5); + } // if OrderReduction != 0 + break; + } // Kratos_Line2D3 || Kratos_Line3D3 + + /** - @ref GeometryData::KratosGeometryType::Kratos_Triangle2D3 "Kratos_Triangle2D3" and @ref GeometryData::KratosGeometryType::Kratos_Triangle3D3 "Kratos_Triangle3D3" + * linear triangles map to themselves. + * \f[\begin{bmatrix} + * 1 & & \\ + * & 1 & \\ + * & & 1 + * \end{bmatrix}\f] + */ + case G::Kratos_Triangle2D3: + case G::Kratos_Triangle3D3: { + *itOutput++ = Triplet(0, 0, 1); + *itOutput++ = Triplet(1, 1, 1); + *itOutput = Triplet(2, 2, 1); + break; + } // Kratos_Triangle2D3 || Kratos_Triangle3D3 + + /** - @ref GeometryData::KratosGeometryType::Kratos_Triangle2D6 "Kratos_Triangle2D6" and @ref GeometryData::KratosGeometryType::Kratos_Triangle3D6 "Kratos_Triangle3D6" + * quadratic triangles map to linear triangles + * @ref GeometryData::KratosGeometryType::Kratos_Triangle2D3 "Kratos_Triangle2D3" and @ref GeometryData::KratosGeometryType::Kratos_Triangle3D3 "Kratos_Triangle3D3". + * \f[\begin{bmatrix} + * 1 & & & \frac{1}{2} & & \frac{1}{2} \\ + * & 1 & & \frac{1}{2} & \frac{1}{2} & \\ + * & & 1 & & \frac{1}{2} & \frac{1}{2} + * \end{bmatrix}\f] + */ + case G::Kratos_Triangle2D6: + case G::Kratos_Triangle3D6: { + if constexpr (OrderReduction == 0) { + *itOutput++ = Triplet(0, 0, 1.0); + *itOutput++ = Triplet(1, 1, 1.0); + *itOutput++ = Triplet(2, 2, 1.0); + *itOutput++ = Triplet(3, 3, 1.0); + *itOutput++ = Triplet(4, 4, 1.0); + *itOutput = Triplet(5, 5, 1.0); + } /*if OrderReduction == 0*/ else { + *itOutput++ = Triplet(0, 0, 1.0); + *itOutput++ = Triplet(0, 3, 0.5); + *itOutput++ = Triplet(0, 5, 0.5); + *itOutput++ = Triplet(1, 1, 1.0); + *itOutput++ = Triplet(1, 3, 0.5); + *itOutput++ = Triplet(1, 4, 0.5); + *itOutput++ = Triplet(2, 2, 1.0); + *itOutput++ = Triplet(2, 4, 0.5); + *itOutput = Triplet(2, 5, 0.5); + } // if OrderReduction != 0 + break; + } // Kratos_Triangle2D6 || Kratos_Triangle3D6 + + /** - @ref GeometryData::KratosGeometryType::Kratos_Tetrahedra3D4 "Kratos_Tetrahedra3D4" + * linear tetrahedron maps to itself. + * \f[\begin{bmatrix} + * 1 & & & \\ + * & 1 & & \\ + * & & 1 & \\ + * & & & 1 + * \end{bmatrix}\f] + */ + case G::Kratos_Tetrahedra3D4: { + *itOutput++ = Triplet(0, 0, 1); + *itOutput++ = Triplet(1, 1, 1); + *itOutput++ = Triplet(2, 2, 1); + *itOutput = Triplet(3, 3, 1); + break; + } // Kratos_Tetrahedra3D4 + + /** - @ref GeometryData::KratosGeometryType::Kratos_Tetrahedra3D10 "Kratos_Tetrahedra3D10" + * quadratic tetrahedron maps to linear tetrahedron @ref GeometryData::KratosGeometryType::Kratos_Tetrahedra3D4 "Kratos_Tetrahedra3D4". + * \f[\begin{bmatrix} + * 1 & & & & \frac{1}{2} & & \frac{1}{2} & \frac{1}{2} & & \\ + * & 1 & & & \frac{1}{2} & \frac{1}{2} & & & \frac{1}{2} & \\ + * & & 1 & & & \frac{1}{2} & \frac{1}{2} & & & \frac{1}{2} \\ + * & & & 1 & & & & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} + * \end{bmatrix}\f] + */ + case G::Kratos_Tetrahedra3D10: { + if constexpr (OrderReduction == 0) { + *itOutput++ = Triplet(0, 0, 1.0); + *itOutput++ = Triplet(1, 1, 1.0); + *itOutput++ = Triplet(2, 2, 1.0); + *itOutput++ = Triplet(3, 3, 1.0); + *itOutput++ = Triplet(4, 4, 1.0); + *itOutput++ = Triplet(5, 5, 1.0); + *itOutput++ = Triplet(6, 6, 1.0); + *itOutput++ = Triplet(7, 7, 1.0); + *itOutput++ = Triplet(8, 8, 1.0); + *itOutput = Triplet(9, 9, 1.0); + } /*if OrderReduction == 0*/ else { + *itOutput++ = Triplet(0, 0, 1.0); + *itOutput++ = Triplet(0, 4, 0.5); + *itOutput++ = Triplet(0, 6, 0.5); + *itOutput++ = Triplet(0, 7, 0.5); + *itOutput++ = Triplet(1, 1, 1.0); + *itOutput++ = Triplet(1, 4, 0.5); + *itOutput++ = Triplet(1, 5, 0.5); + *itOutput++ = Triplet(1, 8, 0.5); + *itOutput++ = Triplet(2, 2, 1.0); + *itOutput++ = Triplet(2, 5, 0.5); + *itOutput++ = Triplet(2, 6, 0.5); + *itOutput++ = Triplet(2, 9, 0.5); + *itOutput++ = Triplet(3, 3, 1.0); + *itOutput++ = Triplet(3, 7, 0.5); + *itOutput++ = Triplet(3, 8, 0.5); + *itOutput = Triplet(3, 9, 0.5); + } + break; + } // Kratos_Tetrahedra3D10 + + default: KRATOS_ERROR << "unsupported geometry: " << rGeometry.Name() << "\n"; + } // switch geometry_type +} + + +namespace detail { + + +inline auto FindNodeIndex(ModelPart& rModelPart, + Node& rNode) noexcept +{ + return std::distance( + rModelPart.Nodes().begin(), + std::lower_bound(rModelPart.Nodes().begin(), + rModelPart.Nodes().end(), + rNode, + [](const Node& r_left, const Node& r_right){ + return r_left.Id() < r_right.Id(); + }) + ); +} + + +struct DofData +{ + NodalData first; + Dof second; +}; // struct DofData + + +} // namespace detail + + +/// @brief Compute the p-multigrid restriction operator for the provided mesh. +/// @tparam TValue Number type of stored values in the sparse matrix. +/// @tparam OrderReduction +/// @param rModelPart @ref ModelPart representing the system to be solved. +/// @param FineSystemSize Number of rows in the fine system's left hand side matrix. +/// @param rRestrictionOperator Output matrix to which the restriction operator will be written. +/// @param rDofSet @ref Output vector containing the @ref Dof "DoFs" of the coarse system. +/// @param rIndirectDofSet An array of pointers pointing to @ref Dof "dofs" in @p rDofSet. +/// @param rDofMap An index map relating fine DoFs to their coarse counterparts. +/// @warning This function assumes that elements use every @ref Dof of their nodes. This may +/// not always be true (e.g.: coupled analyses on overlapping domains and shared +/// @ref Node "nodes" but different Dofs). +template +void MakePRestrictionOperator(ModelPart& rModelPart, + const std::size_t FineSystemSize, + const PointerVectorSet>& rParentIndirectDofSet, + typename TUblasSparseSpace::MatrixType& rRestrictionOperator, + const VariablesList::Pointer& rpVariableList, + std::vector& rDofSet, + PointerVectorSet>& rIndirectDofSet, + std::vector& rDofMap) +{ + static_assert(OrderReduction == std::numeric_limits::max(), + "The current implementation requires geometries to be always reduced to their linear equivalents in a single step."); + + KRATOS_TRY + + using GlobalIndex = typename TUblasSparseSpace::MatrixType::size_type; + using LocalIndex = std::uint8_t; + + // Construct an extended representation of the restriction operator. + // Each entry in the array represents a row in the matrix, with the + // stored map defining its nonzero entries. Empty rows are ignored and + // won't be in the final operator. + // Since the final coarse DoF index assignment can only happen once we + // know which rows have entries and which don't, we can't assign global + // column indices while constructing the restriction operator. Instead, + // fine row indices are stored, from which coarse row indices can later + // be computed. + std::vector, + Dof* // <== pointer to the related DoF on the fine grid + >> rows(FineSystemSize); + + // Construct a COO representation of the restriction operator. + { + std::vector locks(FineSystemSize); + using Triplet = std::tuple; + + // A vector of bools indicating whether the node at the same position + // is a hanging node (i.e.: not part of any element/condition) or not. + // In case you're wondering, I'm not using an std::vector> + // because I don't want someone to change it into a vector of naked + // bools in the future, which would lead to catastrophe. Bools and + // std::uint8_t s are memory and performance-wise identical anyway. + std::vector> hanging_nodes(rModelPart.Nodes().size()); + block_for_each(hanging_nodes, [](auto& r_flag) {r_flag = 1;}); + + struct TLS + { + std::vector local_restriction_operator; + std::vector included_dofs; + }; // struct TLS + + // Collect terms from elements. + KRATOS_TRY + block_for_each(rModelPart.Elements(), + TLS(), + [&rows, &locks, &hanging_nodes, &rParentIndirectDofSet, &rModelPart](const Element& r_element, TLS& r_tls) { + if (r_element.IsActive()) { + const auto& r_geometry = r_element.GetGeometry(); + + // Fetch the local restriction operator of the element. + r_tls.local_restriction_operator.clear(); + MakePRestrictionOperator(r_geometry, std::back_inserter(r_tls.local_restriction_operator)); + + for (const auto& [i_row, i_column, value] : r_tls.local_restriction_operator) { + auto& r_row_dofs = r_geometry[i_row].GetDofs(); + const auto& r_column_dofs = r_geometry[i_column].GetDofs(); + KRATOS_DEBUG_ERROR_IF_NOT(r_row_dofs.size() == r_column_dofs.size()); + + if (!r_row_dofs.empty()) { + // Find the first DoF in the parent's set. + auto it_first_dof = std::lower_bound(rParentIndirectDofSet.begin(), + rParentIndirectDofSet.end(), + **r_row_dofs.begin(), + [](const Dof& r_parent_dof, const Dof& r_dof){ + return r_parent_dof.Id() < r_dof.Id(); + }); + if (it_first_dof == rParentIndirectDofSet.end()) continue; + + // Find which DoFs are included from the node. + r_tls.included_dofs.resize(r_row_dofs.size()); + std::fill(r_tls.included_dofs.begin(), r_tls.included_dofs.end(), false); + + { + const auto i_row_node = r_row_dofs.front()->Id(); + auto it_parent_dof = it_first_dof; + for (std::size_t i_component=0ul; i_component& r_row_dof = *r_row_dofs[i_component]; + if (it_parent_dof->Id() == i_row_node) { + if (it_parent_dof->GetVariable().Key() == r_row_dof.GetVariable().Key()) { + r_tls.included_dofs[i_component] = true; + ++it_parent_dof; + } + } else break; + } // for i_component in range(r_row_dofs.size()) + } + + const unsigned component_count = r_row_dofs.size(); + for (unsigned i_component=0u; i_component* p_fine_row_dof = r_row_dofs[i_component].get(); + [[maybe_unused]] std::scoped_lock row_lock(locks[p_fine_row_dof->EquationId()]); + + const std::size_t i_row_dof = p_fine_row_dof->EquationId(); + rows[i_row_dof].second = p_fine_row_dof; + + const std::size_t i_column_dof = r_column_dofs[i_component]->EquationId(); + [[maybe_unused]] const auto [it_emplace, inserted] = rows[i_row_dof].first.emplace(i_column_dof, value); + + // Check whether the insertion was successful, and if not, + // make sure that the existing restriction component is + // identical to what we tried to insert. + KRATOS_DEBUG_ERROR_IF_NOT(inserted || value == it_emplace->second); + } // for i_component in range(component_count) + } // if !r_row_dofs.empty() + } // for triplet in r_tls + + // Mark nodes as not hanging. + for (Node& r_node : r_element.GetGeometry()) { + const std::size_t i_node = detail::FindNodeIndex(rModelPart, r_node); + KRATOS_DEBUG_ERROR_IF_NOT(i_node < hanging_nodes.size()); + hanging_nodes[i_node] = 0u; + } // for r_node in r_element.GetGeometry() + } // if r_element.IsActive() + }); // for r_element in rModelPart.Elements() + KRATOS_CATCH("") + + // Here comes the tricky part. + // Some models have nodes that are not connected to any element or condition, + // but have multifreedom constraints to other nodes that are. These connected + // external DoFs must also be included in the coarse set of DoFs. + // To find these DoFs, loop through the input DoF set and check whether their + // nodes are hanging. If they are, they fall into the category mentioned above. + KRATOS_TRY + block_for_each(rParentIndirectDofSet, [&hanging_nodes, &locks, &rows, &rModelPart](Dof& r_dof){ + const std::size_t node_id = r_dof.Id(); + const std::size_t i_node = std::distance( + rModelPart.Nodes().begin(), + std::lower_bound(rModelPart.Nodes().begin(), + rModelPart.Nodes().end(), + node_id, + [](const Node& r_left, const std::size_t i_node){ + return r_left.Id() < i_node; + })); + + // Carry the DoFs of such nodes to the coarse system. + if (hanging_nodes[i_node]) { + const std::size_t i_dof = r_dof.EquationId(); + [[maybe_unused]] std::scoped_lock lock(locks[i_dof]); + rows[i_dof].second = &r_dof; + rows[i_dof].first.emplace(i_dof, 1.0); + } // if hanging_nodes[i_node] + }); + KRATOS_CATCH("") + } // destroy locks + + // Construct a fine DoF index => coarse DoF index map. + rDofMap.clear(); + rDofMap.resize(FineSystemSize, std::numeric_limits::max()); + + KRATOS_TRY + std::size_t i_coarse_dof = 0ul; + std::size_t entry_count = 0ul; + for (std::size_t i_fine_dof=0ul; i_fine_dof erase them. + rows.erase(std::remove_if(rows.begin(), + rows.end(), + [](const auto& r_pair) {return r_pair.first.empty();}), + rows.end()); + KRATOS_ERROR_IF_NOT(rows.size() == i_coarse_dof); + + // Resize the restriction operator. + rRestrictionOperator = typename TUblasSparseSpace::MatrixType(i_coarse_dof, FineSystemSize, entry_count); + KRATOS_CATCH("") + + // todo + // Normally, the output DofSet would be filled here but since constructing + // Dofs is a pain in the ass and the current implementation of constraint assemblers + // don't need it, I'll skip this task for now. I'm leaving the DofSet in the arguments though + // because it will eventually become necessary when imposition with lagrange multipliers + // (not augmented ones) is implemented. + // Filling the Dof pointers is relatively easy with the current implementation though. + // This function only allows a direct restriction to a linear mesh, and linear geometries' + // nodes (i.e.: corner nodes) are always a strict subset of their high order counterparts. + // This means that fine Dofs can be reused for the coarse system. + KRATOS_TRY + // Collect solution step variables and store them in a hash map. + std::unordered_map::KeyType,const Variable*> solution_step_variable_map; + + for (const auto& r_variable_data : rModelPart.GetNodalSolutionStepVariablesList()) { + const std::string& r_name = r_variable_data.Name(); + + if (r_variable_data.IsNotComponent()) { + for (auto suffix : {"_X", "_Y", "_Z", "_XX", "_YY", "_ZZ", "_XY", "_YZ", "_XZ"}) { + const std::string component_name = r_name + suffix; + if (KratosComponents>::Has(component_name)) { + const auto key = r_variable_data.Key(); + const Variable& r_variable = KratosComponents>::Get(component_name); + solution_step_variable_map.emplace(key, &r_variable); + } + } + } else { + const auto key = r_variable_data.Key(); + KRATOS_ERROR_IF_NOT(KratosComponents>::Has(r_name)) + << r_name << " is not a registered variable"; + const Variable& r_variable = KratosComponents>::Get(r_name); + solution_step_variable_map.emplace(key, &r_variable); + } + } // for r_variable_data in rModelPart.GetNodalSolutionStepVariablesList + + rDofSet.clear(); + rIndirectDofSet.clear(); + + rDofSet.reserve(rows.size()); + rIndirectDofSet.reserve(rows.size()); + + Dof::IndexType i_dof = 0; + for ([[maybe_unused]] auto& [r_row, rp_dof] : rows) { + { + // The CI's rocky linux is using an ancient version of GCC (8.5) that + // has some bug with non-copyable classes in std::pair, so the following + // line has to be rewritten in a manner that this fossil understands: + //rDofSet.emplace_back(NodalData(rp_dof->Id(), rpVariableList), Dof()); + detail::DofData entry {1ul, Dof()}; + entry.first = NodalData(rp_dof->Id(), rpVariableList); + detail::DofData wtf = std::move(entry); + rDofSet.push_back(std::move(wtf)); + } + rDofSet.back().first.SetSolutionStepData(*rp_dof->GetSolutionStepsData()); + + const auto& r_variable_data = rp_dof->GetVariable(); + auto it_variable = solution_step_variable_map.find(r_variable_data.Key()); + + if (it_variable == solution_step_variable_map.end()) { + const auto key = r_variable_data.Key(); + const std::string& r_name = r_variable_data.Name(); + KRATOS_ERROR_IF_NOT(KratosComponents>::Has(r_name)); + const Variable& r_variable = KratosComponents>::Get(r_name); + it_variable = solution_step_variable_map.emplace(key, &r_variable).first; + } + + const auto& r_reaction_data = rp_dof->GetReaction(); + auto it_reaction = solution_step_variable_map.find(r_reaction_data.Key()); + + if (it_reaction == solution_step_variable_map.end()) { + const auto key = r_reaction_data.Key(); + const std::string& r_name = r_reaction_data.Name(); + KRATOS_ERROR_IF_NOT(KratosComponents>::Has(r_name)); + const Variable& r_variable = KratosComponents>::Get(r_name); + it_reaction = solution_step_variable_map.emplace(key, &r_variable).first; + } + + rDofSet.back().second = Dof(&rDofSet.back().first, + *it_variable->second, + *it_reaction->second); + rDofSet.back().second.SetEquationId(i_dof++); + + if (rp_dof->IsFixed()) { + rDofSet.back().second.FixDof(); + } + + rIndirectDofSet.insert(rIndirectDofSet.end(), &rDofSet.back().second); + } + KRATOS_CATCH("") + + // Fill restriction operator row extents. + // Note: it's a cumulative sum: can't do it in parallel and the standard + // doesn't have an algorithm flexible enough for this purpose. + KRATOS_TRY + rRestrictionOperator.index1_data()[0] = 0; + for (std::size_t i_coarse_dof=0ul; i_coarse_dof>; + IndexPartition(rows.size()).for_each(TLS(), + [&rows, &rRestrictionOperator](const std::size_t i_coarse_dof, + TLS& r_tls) { + const std::size_t i_entry_begin = rRestrictionOperator.index1_data()[i_coarse_dof]; + [[maybe_unused]] const std::size_t i_entry_end = rRestrictionOperator.index1_data()[i_coarse_dof + 1]; + + KRATOS_DEBUG_ERROR_IF(i_entry_end - i_entry_begin != rows[i_coarse_dof].first.size()); + + // Copy and sort the current row. + const auto& r_row = rows[i_coarse_dof].first; + r_tls.resize(r_row.size()); + std::copy(r_row.begin(), r_row.end(), r_tls.begin()); + std::sort(r_tls.begin(), + r_tls.end(), + [](const auto& r_left, const auto& r_right){ + return r_left.first < r_right.first; + }); + + // Copy row into the CSR matrix. + std::size_t i_entry = i_entry_begin; + for (const auto& [i_fine_column, entry] : r_tls) { + KRATOS_DEBUG_ERROR_IF_NOT(i_entry < rRestrictionOperator.index2_data().size()) << i_entry; + KRATOS_DEBUG_ERROR_IF_NOT(i_entry < rRestrictionOperator.value_data().size()) << i_entry; + rRestrictionOperator.index2_data()[i_entry] = i_fine_column; + rRestrictionOperator.value_data()[i_entry] = entry; + ++i_entry; + } // for i_fine_column, entry in r_row + }); // for i_coarse_dof in range(rows.size()) + KRATOS_CATCH("") + + KRATOS_TRY + rRestrictionOperator.set_filled(rRestrictionOperator.index1_data().size(), + rRestrictionOperator.index2_data().size()); + KRATOS_CATCH("") + + KRATOS_CATCH("") +} + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp new file mode 100644 index 000000000000..061142145c07 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/sparse_utilities.hpp @@ -0,0 +1,672 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// Project includes +#include "solving_strategies/schemes/scheme.h" // Scheme +#include "spaces/ublas_space.h" // TUblasSparseSpace +#include "utilities/profiler.h" // KRATOS_PROFILE_SCOPE + +// System includes +#include // std::uint8_t +#include // CHAR_BIT + + +namespace Kratos { + + +// "Hey compiler trust me, my arrays are aligned so try vectorizing my loops over them." +#if defined(__GNUC__) || defined(__clang__) +#define KRATOS_GET_ALIGNED_INDEX_ARRAY(POINTER_TYPE, POINTER_NAME, POINTER_VALUE) \ + POINTER_TYPE __restrict POINTER_NAME = (POINTER_TYPE) __builtin_assume_aligned(POINTER_VALUE, CHAR_BIT * sizeof(POINTER_TYPE)); +#else +// Someone can figure out how to encourage AVX on MSVC, I won't. +#define KRATOS_GET_ALIGNED_INDEX_ARRAY(POINTER_TYPE, POINTER_NAME, POINTER_VALUE) \ + POINTER_TYPE POINTER_NAME = POINTER_VALUE; +#endif + + +/** @brief Computes \f$ A += c * B \f$. + * @tparam TValue Value type of the operand matrices. + * @param rLeft Left hand side matrix to add the scaled right hand matrix to. + * @param rRight Right hand matrix to scale and add to the left hand side. + * @param Coefficient Coefficient to scale the right hand side matrix by. + * @warning The sparsity pattern of @p rLeft must be a superset of @p rRight. + */ +template +static void +InPlaceMatrixAdd(typename TUblasSparseSpace::MatrixType& rLeft, + const typename TUblasSparseSpace::MatrixType& rRight, + const TValue Coefficient = 1.0) +{ + // Sanity checks. + KRATOS_ERROR_IF_NOT(rLeft.size1() == rRight.size1() && rLeft.size2() == rRight.size2()) + << "cannot add matrices of different shapes " + << "(" << rLeft.size1() << "x" << rLeft.size2() << ") + " + << "(" << rRight.size1() << "x" << rRight.size2() << ")"; + + // Early exit on empty update matrix. + if (!rRight.size1() || !rRight.size2() || !rRight.nnz()) return; + + KRATOS_ERROR_IF_NOT(rRight.nnz() <= rLeft.nnz()); + + #define KRATOS_MATRIX_SUM(sum_operator) \ + IndexPartition(rLeft.size1()).for_each([&rLeft, &rRight, Coefficient](const IndexType i_row) { \ + (void)Coefficient; /*<== suppress unused capture warnings.*/ \ + const auto i_right_entry_begin = rRight.index1_data()[i_row]; \ + const auto i_right_entry_end = rRight.index1_data()[i_row + 1]; \ + \ + const auto i_left_entry_begin = rLeft.index1_data()[i_row]; \ + const auto i_left_entry_end = rLeft.index1_data()[i_row + 1]; \ + const auto it_left_column_begin = rLeft.index2_data().begin() + i_left_entry_begin; \ + const auto it_left_column_end = rLeft.index2_data().begin() + i_left_entry_end; \ + \ + auto it_left_column = it_left_column_begin; \ + for (auto i_right_entry=i_right_entry_begin; i_right_entry!=i_right_entry_end; ++i_right_entry) { \ + const auto i_column = rRight.index2_data()[i_right_entry]; \ + \ + /* Find the entry in the left matrix corresponding to the entry in the right one. */ \ + it_left_column = std::lower_bound(it_left_column, it_left_column_end, i_column); \ + KRATOS_ERROR_IF(it_left_column == it_left_column_end || *it_left_column != i_column) \ + << "left hand matrix has no entry in row " << i_row << " at column " << i_column; \ + \ + const auto i_left_entry = std::distance(rLeft.index2_data().begin(), it_left_column); \ + rLeft.value_data()[i_left_entry] sum_operator rRight.value_data()[i_right_entry]; \ + } /* for i_right_entry in range(i_right_begin, i_right_end) */ \ + }) /*for i_row in range(rLeft.size1())*/ + + if (Coefficient == 1.0) { + KRATOS_TRY + KRATOS_MATRIX_SUM(+=); + KRATOS_CATCH("") + } /*if Coefficient == 1.0*/ else if (Coefficient == -1.0) { + KRATOS_TRY + KRATOS_MATRIX_SUM(-=); + KRATOS_CATCH("") + } /*if Coefficient == -1.0*/ else { + KRATOS_TRY + KRATOS_MATRIX_SUM(+= Coefficient *); + KRATOS_CATCH("") + } // if Coefficient != 1.0 + #undef KRATOS_MATRIX_SUM +} + + +/// @brief Compute the union of two sparse matrix' sparsity patterns. +template +static void +MergeMatrices(typename TUblasSparseSpace::MatrixType& rLeft, + const typename TUblasSparseSpace::MatrixType& rRight) +{ + using SpaceType = TUblasSparseSpace; + using MatrixType = typename SpaceType::MatrixType; + using IndexType = typename MatrixType::index_array_type::value_type; + MatrixType output; + + KRATOS_TRY + + // Sanity checks. + KRATOS_ERROR_IF_NOT(rLeft.size1() == rRight.size1() && rLeft.size2() == rRight.size2()) + << "cannot merge incompatible matrices " + << "(" << rLeft.size1() << "x" << rLeft.size2() << ")" + << " and " + << "(" << rRight.size1() << "x" << rRight.size2() << ")"; + + if (!rRight.size1() || !rRight.size2() || !rRight.nnz()) return; + + // Declare new containers for the merged matrix. + typename MatrixType::index_array_type row_extents(rLeft.index1_data().size()); + typename MatrixType::index_array_type column_indices; + block_for_each(row_extents, [](auto& r_item){r_item = 0;}); + + // Merge rows into separate containers. + { + std::vector> rows(rLeft.size1()); + IndexPartition(rLeft.size1()).for_each([&rows, &rLeft, &rRight](const IndexType i_row){ + const IndexType i_left_begin = rLeft.index1_data()[i_row]; + const IndexType i_left_end = rLeft.index1_data()[i_row + 1]; + const IndexType i_right_begin = rRight.index1_data()[i_row]; + const IndexType i_right_end = rRight.index1_data()[i_row + 1]; + rows[i_row].reserve(i_left_end - i_left_begin + i_right_end - i_right_begin); + + rows[i_row].insert(rows[i_row].end(), + rLeft.index2_data().begin() + i_left_begin, + rLeft.index2_data().begin() + i_left_end); + rows[i_row].insert(rows[i_row].end(), + rRight.index2_data().begin() + i_right_begin, + rRight.index2_data().begin() + i_right_end); + std::sort(rows[i_row].begin(), + rows[i_row].end()); + rows[i_row].erase(std::unique(rows[i_row].begin(), + rows[i_row].end()), + rows[i_row].end()); + rows[i_row].shrink_to_fit(); + }); + + // Compute new row extents. + for (IndexType i_row=0; i_row(rLeft.size1()).for_each([&rows, &row_extents, &column_indices](const IndexType i_row){ + const IndexType i_begin = row_extents[i_row]; + std::copy(rows[i_row].begin(), + rows[i_row].end(), + column_indices.begin() + i_begin); + }); + } + + // Construct the new matrix. + rLeft = MatrixType(rLeft.size1(), rLeft.size2(), column_indices.size()); + rLeft.index1_data().swap(row_extents); + rLeft.index2_data().swap(column_indices); + block_for_each(rLeft.value_data(), [](auto& r_item){r_item = 0;}); + rLeft.set_filled(rLeft.size1() + 1, column_indices.size()); + + KRATOS_CATCH("") +} + + +template +void MakeSparseTopology(TRowMapContainer& rRows, + const IndexType ColumnCount, + typename TUblasSparseSpace::MatrixType& rMatrix, + bool EnsureDiagonal) +{ + KRATOS_TRY + + if (EnsureDiagonal) { + IndexPartition(rRows.size()).for_each([&rRows](IndexType i_row){ + rRows[i_row].emplace(i_row); + }); + } + + const IndexType row_count = rRows.size(); + IndexType entry_count = 0ul; + + { + auto row_extents = rMatrix.index1_data(); + + // Resize row extents. + row_extents.resize(row_count + 1, false); + + // Fill row extents and collect the total number of entries to store. + row_extents[0] = 0; + for (int i = 0; i < static_cast(row_count); i++) { + row_extents[i + 1] = row_extents[i] + rRows[i].size(); + entry_count += rRows[i].size(); + } // for i in range(row_count) + + // Resize the output matrix and all its containers. + rMatrix = typename TUblasSparseSpace::MatrixType(rRows.size(), ColumnCount, entry_count); + rMatrix.index1_data().swap(row_extents); + } + + auto& r_row_extents = rMatrix.index1_data(); + auto& r_column_indices = rMatrix.index2_data(); + + // Copy column indices. + IndexPartition(row_count).for_each([&r_row_extents, &r_column_indices, &rRows](const IndexType i_row){ + const IndexType i_entry_begin = r_row_extents[i_row]; + const IndexType i_entry_end = r_row_extents[i_row + 1]; + KRATOS_DEBUG_ERROR_IF(r_column_indices.size() < i_entry_end); + std::copy(rRows[i_row].begin(), + rRows[i_row].end(), + r_column_indices.begin() + i_entry_begin); + + if constexpr (!SortedRows) { + std::sort(r_column_indices.begin() + i_entry_begin, + r_column_indices.begin() + i_entry_end); + } + }); + + KRATOS_TRY + block_for_each(rMatrix.value_data(), [](TValue& r_entry){r_entry=0;}); + rMatrix.set_filled(row_count + 1, entry_count); + KRATOS_CATCH("") + + KRATOS_CATCH("") +} + + +/// @brief Map LHS contributions from a local dense matrix to a sparse global matrix. +template +void MapRowContribution(typename TSparse::MatrixType& rLhs, + const typename TDense::MatrixType& rLocalLhs, + const IndexType iRow, + const IndexType iLocalRow, + const Element::EquationIdVectorType& rEquationIds) noexcept +{ + auto& r_entries = rLhs.value_data(); + const auto& r_row_extents = rLhs.index1_data(); + const auto& r_column_indices = rLhs.index2_data(); + + const auto i_row_begin = r_row_extents[iRow]; + const auto i_row_end = r_row_extents[iRow + 1]; + const typename TDense::IndexType row_size = rEquationIds.size(); + + for (typename TDense::IndexType i_local_column=0; i_local_column +void MapContribution(typename TSparse::VectorType& rRhs, + const typename TDense::VectorType& rContribution, + const Element::EquationIdVectorType& rEquationIds) noexcept +{ + const typename TDense::IndexType local_size = rContribution.size(); + + for (typename TDense::IndexType i_local = 0; i_local < local_size; i_local++) { + const typename TSparse::IndexType i_global = rEquationIds[i_local]; + AtomicAdd(rRhs[i_global], static_cast(rContribution[i_local])); + } +} + + +/// @brief Map LHS contributions from local to global space. +template +void MapContribution(typename TSparse::MatrixType& rLhs, + const typename TDense::MatrixType& rContribution, + const Element::EquationIdVectorType& rEquationIds, + LockObject* pLockBegin) noexcept +{ + const typename TDense::IndexType local_size = rContribution.size1(); + for (typename TDense::IndexType i_local = 0; i_local < local_size; i_local++) { + const IndexType i_global = rEquationIds[i_local]; + std::scoped_lock lock(pLockBegin[i_global]); + MapRowContribution(rLhs, + rContribution, + i_global, + i_local, + rEquationIds); + } +} + + +/// @brief Map LHS and RHS contributions from local to global space. +template +void MapContribution(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + const typename TDense::MatrixType& rLhsContribution, + const typename TDense::VectorType& rRhsContribution, + const Element::EquationIdVectorType& rEquationIds, + LockObject* pLockBegin) noexcept +{ + const typename TDense::IndexType local_size = rLhsContribution.size1(); + for (typename TDense::IndexType i_local = 0; i_local < local_size; i_local++) { + const typename TSparse::IndexType i_global = rEquationIds[i_local]; + std::scoped_lock lock(pLockBegin[i_global]); + rRhs[i_global] += rRhsContribution[i_local]; + MapRowContribution(rLhs, + rLhsContribution, + i_global, + i_local, + rEquationIds); + } +} + + +/// @brief Compute the local contributions of an @ref Element or @ref Condition and assemble them into the global system. +/// @details Whether from elements or conditions, whether assembling into the LHS, RHS, or both, +/// all cases are handled in this one function to reduce code duplication and have logically +/// coherent parts of the code in one place. +template +void MapEntityContribution(TEntity& rEntity, + Scheme& rScheme, + const ProcessInfo& rProcessInfo, + Element::EquationIdVectorType& rEquationIndices, + typename TDense::MatrixType* pLhsContribution, + typename TDense::VectorType* pRhsContribution, + typename TSparse::MatrixType* pLhs, + typename TSparse::VectorType* pRhs, + LockObject* pLockBegin) +{ + if constexpr (AssembleLHS || AssembleRHS) { + if (rEntity.IsActive()) { + if constexpr (AssembleLHS) { + if constexpr (AssembleRHS) { + // Assemble LHS and RHS. + rScheme.CalculateSystemContributions(rEntity, + *pLhsContribution, + *pRhsContribution, + rEquationIndices, + rProcessInfo); + MapContribution(*pLhs, + *pRhs, + *pLhsContribution, + *pRhsContribution, + rEquationIndices, + pLockBegin); + } /*if AssembleRHS*/ else { + // Assemble LHS only. + rScheme.CalculateLHSContribution(rEntity, + *pLhsContribution, + rEquationIndices, + rProcessInfo); + MapContribution(*pLhs, + *pLhsContribution, + rEquationIndices, + pLockBegin); + } // if !AssembleRHS + } /*if AssembleLHS*/ else { + // Assemble RHS only. + rScheme.CalculateRHSContribution(rEntity, + *pRhsContribution, + rEquationIndices, + rProcessInfo); + MapContribution(*pRhs, + *pRhsContribution, + rEquationIndices); + } // if !AssembleLHS + } // if rEntity.IsActive + } // if AssembleLHS or AssembleRHS +} + + +template +void ApplyDirichletConditions(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + const TItRowDof itRowDofBegin, + const TItRowDof itRowDofEnd, + const TItColumnDof itColumnDofBegin, + [[maybe_unused]] const TItColumnDof itColumnDofEnd, + const typename TSparse::DataType DiagonalScaleFactor) +{ + // Type checks. + static_assert(std::is_same_v< + typename std::iterator_traits::value_type, + Dof>); + + KRATOS_TRY + KRATOS_PROFILE_SCOPE(KRATOS_CODE_LOCATION); + + block_for_each(itRowDofBegin, + itRowDofEnd, + [&rLhs, &rRhs, itColumnDofBegin, DiagonalScaleFactor](const Dof& r_dof){ + const IndexType i_dof = r_dof.EquationId(); + const typename TSparse::IndexType i_entry_begin = rLhs.index1_data()[i_dof]; + const typename TSparse::IndexType i_entry_end = rLhs.index1_data()[i_dof + 1]; + bool found_diagonal = false; + + if (r_dof.IsFixed()) { + // Zero out the whole row, except the diagonal. + for (typename TSparse::IndexType i_entry=i_entry_begin; i_entry(0); + } + } // for i_entry in range(i_entry_begin, i_entry_end) + + rRhs[i_dof] = static_cast(0); + } /*if r_dof.IsFixed()*/ else { + // Zero out the column which is associated with the zero'ed row. + for (typename TSparse::IndexType i_entry=i_entry_begin; i_entryIsFixed()) { + rLhs.value_data()[i_entry] = static_cast(0); + } + + } // for i_entry in range(i_entry_begin, i_entry_end) + } /*not r_dof.IsFixed()*/ + + KRATOS_ERROR_IF_NOT(found_diagonal) + << "missing diagonal in row " << i_dof << " " + << "related to dof " << r_dof.GetVariable().Name() << " " + << "of node " << r_dof.Id(); + }); + + KRATOS_CATCH("") +} + + +template +void ApplyDirichletConditions(typename TSparse::MatrixType& rLhs, + typename TSparse::VectorType& rRhs, + const TItDof itDofBegin, + const TItDof itDofEnd, + const typename TSparse::DataType DiagonalScaleFactor) +{ + return ApplyDirichletConditions( + rLhs, + rRhs, + itDofBegin, + itDofEnd, + itDofBegin, + itDofEnd, + DiagonalScaleFactor); +} + + +/// @internal +struct MatrixChecks +{ + static constexpr std::uint8_t None = 0; + static constexpr std::uint8_t DiagonalExists = 1; + static constexpr std::uint8_t DiagonalIsNonNegative = 1 << 1; + static constexpr std::uint8_t DiagonalIsPositive = 1 << 2; + static constexpr std::uint8_t IsDiagonallyDominant = 1 << 3; + static constexpr std::uint8_t RowsAreSorted = 1 << 4; + static constexpr std::uint8_t ColumnsAreSorted = 1 << 5; + static constexpr std::uint8_t All = 1 + + (1 << 1) + + (1 << 2) + + (1 << 3) + + (1 << 4) + + (1 << 5); +}; // struct MatrixCheck + + +/// @internal +template +void CheckMatrix(const typename TUblasSparseSpace::MatrixType& rMatrix) +{ + if constexpr (Checks == MatrixChecks::None) return; + + KRATOS_ERROR_IF_NOT(rMatrix.size1() + 1 == rMatrix.index1_data().size()) + << "input matrix has inconsistent row extents"; + KRATOS_ERROR_IF_NOT(rMatrix.index1_data()[rMatrix.size1()] == rMatrix.nnz()) + << "row extents of the input matrix do not consistently cover its contents"; + + if constexpr (Checks & MatrixChecks::ColumnsAreSorted) { + KRATOS_ERROR_IF_NOT(std::is_sorted(rMatrix.index1_data().begin(), rMatrix.index1_data().end())) + << "the input matrix' columns are not sorted"; + } + + for (IndexType i_row=0ul; i_row maybe_diagonal; + + for (auto i_entry=i_entry_begin; i_entry(0) <= diagonal) + << "diagonal in row " << i_row << " is negative (" << diagonal << ")"; + + if constexpr (Checks & MatrixChecks::DiagonalIsPositive) + KRATOS_ERROR_IF_NOT(static_cast(0) < diagonal) + << "diagonal in row " << i_row << " is not positive (" << diagonal << ")"; + + if constexpr (Checks & MatrixChecks::IsDiagonallyDominant) { + const auto abs_diagonal = std::abs(diagonal); + for (auto i_entry=i_entry_begin; i_entry +void BalancedProduct(const typename TLHSSparse::MatrixType& rLhs, + const typename TRHSSparse::VectorType& rRhs, + typename TOutputSparse::VectorType& rOutput, + const typename TOutputSparse::DataType Coefficient = static_cast(1)) +{ + KRATOS_TRY + + // Create partition for entries in the matrix. + const auto thread_count = ParallelUtilities::GetNumThreads(); + std::vector partition(thread_count + 1); + + partition.front() = 0ul; + const auto chunk_size = rLhs.nnz() / partition.size(); + for (typename TLHSSparse::IndexType i_end=1ul; i_end(thread_count).for_each( \ + [&rLhs, &rRhs, &rOutput, &partition, Coefficient](const typename TLHSSparse::IndexType i_chunk){ \ + (void)Coefficient; /*<== suppress unused capture warnings*/ \ + /* Define the entry range to compute on. */ \ + const auto i_chunk_begin = partition[i_chunk]; \ + const auto i_chunk_end = partition[i_chunk + 1]; \ + \ + /* Find the initial row's index. */ \ + const auto it_initial_row = std::lower_bound(rLhs.index1_data().begin(), \ + rLhs.index1_data().end(), \ + static_cast(i_chunk_begin)); \ + typename TLHSSparse::IndexType i_row = std::distance(rLhs.index1_data().begin(), it_initial_row); \ + if (rLhs.index1_data()[i_row] != i_chunk_begin) --i_row; \ + \ + do { \ + const auto i_row_begin = rLhs.index1_data()[i_row]; \ + const auto i_row_end = rLhs.index1_data()[i_row + 1]; \ + \ + const auto i_begin = std::max(i_row_begin, i_chunk_begin); \ + const auto i_end = std::min(i_row_end, i_chunk_end); \ + const typename TLHSSparse::IndexType chunk_size = i_end - i_begin; \ + \ + auto contribution = static_cast(0); \ + \ + KRATOS_GET_ALIGNED_INDEX_ARRAY(const typename TLHSSparse::IndexType*, it_column, &*(rLhs.index2_data().begin() + i_begin)); \ + KRATOS_GET_ALIGNED_INDEX_ARRAY(const typename TLHSSparse::DataType*, it_entry, &*(rLhs.value_data().begin() + i_begin)); \ + KRATOS_GET_ALIGNED_INDEX_ARRAY(const typename TRHSSparse::DataType*, it_rhs, &*rRhs.begin()); \ + \ + for (typename TLHSSparse::IndexType i=0; i>) { \ + AtomicAdd(rOutput[i_row], contribution); \ + } else if constexpr (std::is_same_v>) { \ + AtomicSub(rOutput[i_row], contribution); \ + } else if constexpr (std::is_same_v>) { \ + AtomicAdd(rOutput[i_row], Coefficient * contribution); \ + } else { \ + static_assert(std::is_same_v, "unsupported operator"); \ + } \ + \ + ++i_row; \ + if (i_end == i_chunk_end) \ + break; \ + } while (true); \ + }) /* for i_chunk in range(thread_count) */ + + if (Coefficient == static_cast(1)) { + KRATOS_BALANCED_MATRIX_VECTOR_PRODUCT(std::plus); + } else if (Coefficient == static_cast(-1)) { + KRATOS_BALANCED_MATRIX_VECTOR_PRODUCT(std::minus); + } else { + KRATOS_BALANCED_MATRIX_VECTOR_PRODUCT(std::multiplies); + } + + #undef KRATOS_BALANCED_MATRIX_VECTOR_PRODUCT + KRATOS_CATCH("") +} + + +#undef KRATOS_GET_ALIGNED_INDEX_ARRAY + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/status_stream.cpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/status_stream.cpp new file mode 100644 index 000000000000..53af788cb046 --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/status_stream.cpp @@ -0,0 +1,351 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +// Project includes +#include "solving_strategies/builder_and_solvers/p_multigrid/status_stream.hpp" // PMGStatusStream +#include "includes/kratos_components.h" +#include "solving_strategies/builder_and_solvers/p_multigrid/p_multigrid_builder_and_solver.hpp" // PMultigridBuilderAndSolver +#include "spaces/ublas_space.h" // TUblasSparseSpace +#include "input_output/vtu_output.h" // VtuOutput +#include "includes/model_part.h" // ModelPart + +// System includes +#include // std::cout +#include // std::stringstream +#include // std::setw, std::setprecision, std::scientific +#include // std::variant +#include // std::tuple +#include // std::filesystem::path, std::filesystem::make_directory, std::filesystem::is_directory +#include // std::unordered_set + + +namespace Kratos { + + +using DoubleBnS = PMultigridBuilderAndSolver, + TUblasDenseSpace>; + + +using SingleBnS = PMultigridBuilderAndSolver, + TUblasDenseSpace>; + + +std::unique_ptr MakeVtuOutput(ModelPart& rModelPart, + const PointerVectorSet>& rDofSet) +{ + auto p_output = std::make_unique(rModelPart); + + KRATOS_TRY + + std::unordered_set variable_names; + for (const Dof& r_dof : rDofSet) { + variable_names.emplace(r_dof.GetVariable().Name()); + variable_names.emplace(r_dof.GetReaction().Name()); + } + + std::vector*> variables(variable_names.size()); + for (const std::string& r_variable_name : variable_names){ + KRATOS_ERROR_IF_NOT(KratosComponents>::Has(r_variable_name)) + << r_variable_name << " is not a registered variable name"; + p_output->AddHistoricalVariable(KratosComponents>::Get(r_variable_name)); + } + + KRATOS_CATCH("") + + return p_output; +} + + +struct PMGStatusStream::Impl { + void PrintHeader() { + this->PrintHorizontalLine(); + (*mpStream) << "| Grid | Const It | Const Res | It | Res |\n"; + this->PrintHorizontalLine(); + } + + void PrintHorizontalLine() { + // "| Grid | Const It | Const Res | It | Res |\n"; + (*mpStream) << "+ ---- + -------- + ----------- + ---- + ----------- +\n"; + } + + void PrintReport(const PMGStatusStream::Report& rReport) { + std::ostream& r_stream = (*mpStream); + std::stringstream tmp_stream; + + constexpr auto unconverged_color = "\033[0;31m"; //< red + constexpr auto converged_color = "\033[0;32m"; //< green + + // Grid level. + tmp_stream << std::setw(4) << rReport.grid_level; + r_stream << "| " << tmp_stream.str() << " "; + + // Constraint iteration. + tmp_stream = std::stringstream(); + tmp_stream << std::setw(8) << rReport.constraint_iteration; + r_stream << "| " << tmp_stream.str() << " "; + + // Constraint residual. + if (rReport.maybe_constraint_residual.has_value()) { + tmp_stream = std::stringstream(); + if (mUseAnsiColors) { + if (rReport.constraints_converged) { + tmp_stream << converged_color; + } else { + tmp_stream << unconverged_color; + } + } + tmp_stream << std::setw(11) << std::setprecision(4) << std::scientific << rReport.maybe_constraint_residual.value(); + if (mUseAnsiColors) tmp_stream << "\033[0m"; + r_stream << "| " << tmp_stream.str() << " "; + } else { + r_stream << "| "; + } + + // Multigrid iteration. + tmp_stream = std::stringstream(); + tmp_stream << std::setw(4) << rReport.multigrid_iteration; + r_stream << "| " << tmp_stream.str() << " "; + + // Multigrid residual. + tmp_stream = std::stringstream(); + if (mUseAnsiColors) { + if (rReport.multigrid_converged) { + tmp_stream << converged_color; + } else { + tmp_stream << unconverged_color; + } + } + + if (rReport.maybe_multigrid_residual.has_value()) + tmp_stream << std::setw(11) << std::setprecision(4) << std::scientific << rReport.maybe_multigrid_residual.value(); + else + tmp_stream << " "; + + if (mUseAnsiColors) tmp_stream << "\033[0m"; + r_stream << "| " << tmp_stream.str() << " "; + + r_stream << "|\n"; + } + + ~Impl() + { + if (mMaybeLastReport.has_value()) + this->PrintHorizontalLine(); + + if (!mVtuPaths.empty()) { + Parameters contents(R"({"file-series-version" : "1.0", "files" : []})"); + Parameters files = contents["files"]; + for (std::size_t i_path=0ul; i_pathName() + ".vtu.series")); + series_journal << contents.PrettyPrintJsonString(); + } + } + + /// @brief Returns true if the input report is identical to the last one that got issued. + bool DuplicateFilter(const PMGStatusStream::Report& rReport) + { + if (!mMaybeLastReport.has_value()) { + this->PrintHeader(); + mMaybeLastReport = rReport; + return false; + } + + PMGStatusStream::Report& r_last = mMaybeLastReport.value(); + bool is_duplicate = rReport.constraint_iteration == r_last.constraint_iteration + && rReport.constraints_converged == r_last.constraints_converged + && rReport.grid_level == r_last.grid_level + && rReport.multigrid_converged == r_last.multigrid_converged + && rReport.multigrid_iteration == r_last.multigrid_iteration; + + if (rReport.maybe_multigrid_residual.has_value()) { + if (r_last.maybe_multigrid_residual.has_value()) { + is_duplicate = is_duplicate && rReport.maybe_multigrid_residual.value() == r_last.maybe_multigrid_residual.value(); + } else is_duplicate = false; + } else if (r_last.maybe_multigrid_residual.has_value()) is_duplicate = false; + + if (rReport.maybe_constraint_residual.has_value()) { + if (r_last.maybe_constraint_residual.has_value()) { + is_duplicate = is_duplicate && rReport.maybe_constraint_residual.value() == r_last.maybe_constraint_residual.value(); + } else is_duplicate = false; + } else if (r_last.maybe_constraint_residual.has_value()) is_duplicate = false; + + if (is_duplicate) return true; + else { + mMaybeLastReport = rReport; + return false; + } + } + + /// @brief Write state field and residuals while solving the system. + /// @details Any time a report is submitted (and the verbosity is high enough), + /// project the current solution from the submitter's grid onto the + /// root grid, and write the state field and residuals to a VTU file. + void WriteIntermediateState(const PMGStatusStream::Report& rReport) + { + if (!mpMaybeVtuOutput.has_value()) { + std::visit([this](const auto* p_builder_and_solver) { + mpMaybeVtuOutput = MakeVtuOutput(*mpModelPart, p_builder_and_solver->GetDofSet()); + }, mpBuilderAndSolver); + } + + // PMultigridBuilderAndSolver::ProjectGrid stores the state and residual values + // in the Dofs' values and reactions that otherwise should not be overwritten. + // => store the original values and restore them after the VTU output is done. + std::vector state, reaction; + std::visit([&state, &reaction](const auto* p_builder_and_solver){ + const auto& r_dof_set = p_builder_and_solver->GetDofSet(); + state.resize(r_dof_set.size()); + reaction.resize(r_dof_set.size()); + IndexPartition(r_dof_set.size()).for_each([&r_dof_set, &state, &reaction](std::size_t i_dof){ + const Dof& r_dof = *(r_dof_set.begin() + i_dof); + state[i_dof] = r_dof.GetSolutionStepValue(); + reaction[i_dof] = r_dof.GetSolutionStepReactionValue(); + }); // for i_dof in range(len(r_dof_set)) + }, mpBuilderAndSolver); + + KRATOS_TRY + std::visit([this, &rReport](auto* p_builder_and_solver) { + std::visit([&rReport, p_builder_and_solver] (const auto tuple) { + using BSMatrixType = typename std::pointer_traits::element_type::TSystemMatrixType; + using SystemMatrixPointer = std::tuple_element_t<0,decltype(tuple)>; + using SystemMatrix = std::remove_cv_t>::element_type>; + if constexpr (std::is_same_v) + p_builder_and_solver->ProjectGrid(rReport.grid_level, + *std::get<0>(tuple), + *std::get<1>(tuple), + *std::get<2>(tuple)); + else (void)(rReport); + }, mSystem); + }, mpBuilderAndSolver); + KRATOS_CATCH("") + + const std::filesystem::path directory("builder_and_solver"); + if (!std::filesystem::is_directory(directory)) + std::filesystem::create_directories(directory); + + const std::string file_name = mpModelPart->Name() + "_" + std::to_string(mVtuPaths.size()); + mpMaybeVtuOutput.value()->PrintOutput((directory / file_name).string()); + mVtuPaths.emplace_back(file_name + ".vtu"); + + // Restore Dof values and their corresponding reactions. + std::visit([&state, &reaction](auto* p_builder_and_solver){ + auto& r_dof_set = p_builder_and_solver->GetDofSet(); + IndexPartition(r_dof_set.size()).for_each([&r_dof_set, &state, &reaction](std::size_t i_dof){ + Dof& r_dof = *(r_dof_set.begin() + i_dof); + r_dof.GetSolutionStepValue() = state[i_dof]; + r_dof.GetSolutionStepReactionValue() = reaction[i_dof]; + }); // for i_dof in range(len(r_dof_set)) + }, mpBuilderAndSolver); + } + + std::ostream* mpStream; + + ModelPart* mpModelPart; + + using DoubleSpace = TUblasSparseSpace; + using SingleSpace = TUblasSparseSpace; + std::variant< + std::tuple, + std::tuple + > mSystem; + + std::variant mpBuilderAndSolver; + + std::optional> mpMaybeVtuOutput; + + std::vector mVtuPaths; + + std::optional mMaybeLastReport; + + bool mUseAnsiColors; +}; // struct PMGStatusStream::Impl + + +template +PMGStatusStream::PMGStatusStream(std::ostream& rStream, + PMultigridBuilderAndSolver& rBuilderAndSolver, + ModelPart& rModelPart, + const TMatrix& rRootLhs, + const TVector& rRootSolution, + const TVector& rRootRhs, + bool UseAnsiColors) + : mpImpl(new Impl {/*mpStream=*/ &rStream, + /*mpModelPart=*/ &rModelPart, + /*mSystem=*/ std::make_tuple(&rRootLhs, &rRootSolution, &rRootRhs), + /*mpBuilderAndSolver=*/ &rBuilderAndSolver, + /*mpMaybeVtuOutput=*/ {}, + /*mVtuPaths=*/ {}, + /*mMaybeLastReport=*/ {}, + /*mUseAnsiColors=*/ UseAnsiColors}) +{ +} + + +PMGStatusStream::PMGStatusStream(PMGStatusStream&&) noexcept = default; + + +PMGStatusStream::~PMGStatusStream() = default; + + +PMGStatusStream& PMGStatusStream::operator=(PMGStatusStream&&) noexcept = default; + + +void PMGStatusStream::Submit(const TaggedReport& rReport, int SubmitterVerbosity) +{ + const auto& [report, report_verbosity] = rReport; + if (report_verbosity <= SubmitterVerbosity && !mpImpl->DuplicateFilter(report)) { + (*this) << report; + + if (5 <= SubmitterVerbosity) + mpImpl->WriteIntermediateState(report); + } +} + + +PMGStatusStream& PMGStatusStream::operator<<(const Report& rReport) +{ + mpImpl->PrintReport(rReport); + return *this; +} + + +template PMGStatusStream::PMGStatusStream(std::ostream&, + DoubleBnS&, + ModelPart&, + const typename TUblasSparseSpace::MatrixType&, + const typename TUblasSparseSpace::VectorType&, + const typename TUblasSparseSpace::VectorType&, + bool); + + +template PMGStatusStream::PMGStatusStream(std::ostream&, + SingleBnS&, + ModelPart&, + const typename TUblasSparseSpace::MatrixType&, + const typename TUblasSparseSpace::VectorType&, + const typename TUblasSparseSpace::VectorType&, + bool); + + +} // namespace Kratos diff --git a/kratos/solving_strategies/builder_and_solvers/p_multigrid/status_stream.hpp b/kratos/solving_strategies/builder_and_solvers/p_multigrid/status_stream.hpp new file mode 100644 index 000000000000..1cc39ed843df --- /dev/null +++ b/kratos/solving_strategies/builder_and_solvers/p_multigrid/status_stream.hpp @@ -0,0 +1,86 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Máté Kelemen +// + +#pragma once + +// System includes +#include // std::ostream +#include // std::unique_ptr +#include // std::optional +#include // std::pair + + +namespace Kratos { + + +class ModelPart; + +template +class PMultigridBuilderAndSolver; + + +/// @internal +class PMGStatusStream final { +public: + struct Report { + std::size_t grid_level; + bool multigrid_converged; + std::size_t multigrid_iteration; + std::optional maybe_multigrid_residual; + bool constraints_converged; + std::size_t constraint_iteration; + std::optional maybe_constraint_residual; + + std::pair Tag(int ReportVerbosity) + { + return std::make_pair(*this, ReportVerbosity); + } + }; // struct Report + + using TaggedReport = std::pair; + + template + PMGStatusStream(std::ostream& rStream, + PMultigridBuilderAndSolver& rBuilderAndSolver, + ModelPart& rModelPart, + const TMatrix& rRootLhs, + const TVector& rRootSolution, + const TVector& rRootRhs, + bool UseAnsiColors = true); + + PMGStatusStream(PMGStatusStream&&) noexcept; + + ~PMGStatusStream(); + + PMGStatusStream& operator=(PMGStatusStream&&) noexcept; + + /// @brief Submit a @ref Report to log. + /// @param rReport Report to submit. + /// @param SubmitterVerbosity Verbosity of the class that submitted the report. + void Submit(const TaggedReport& rReport, int SubmitterVerbosity); + +private: + PMGStatusStream& operator<<(const Report& rReport); + + PMGStatusStream(const PMGStatusStream&) = delete; + + PMGStatusStream& operator=(const PMGStatusStream&) = delete; + + struct Impl; + std::unique_ptr mpImpl; +}; // class PMGStatusStream + + +} // namespace Kratos diff --git a/kratos/sources/kratos_application.cpp b/kratos/sources/kratos_application.cpp index 033400022995..0e5bf04b5887 100644 --- a/kratos/sources/kratos_application.cpp +++ b/kratos/sources/kratos_application.cpp @@ -63,6 +63,7 @@ KratosApplication::KratosApplication(const std::string& ApplicationName) // Master-Slave Constraint mMasterSlaveConstraint(), mLinearMasterSlaveConstraint(), + mLinearMultifreedomConstraint(), // Periodic conditions mPeriodicCondition( 0, GeometryType::Pointer(new Line2D2(GeometryType::PointsArrayType(2)))), @@ -170,11 +171,12 @@ void KratosApplication::RegisterKratosCore() { Serializer::Register("DofDouble", Dof()); Serializer::Register("MasterSlaveConstraint", MasterSlaveConstraint()); + Serializer::Register("MultifreedomConstraint", MultifreedomConstraint()); //Register specific conditions ( must be completed : conditions defined in kratos_application.h) //generic condition KRATOS_REGISTER_CONDITION("GenericCondition", mGenericCondition); - + // Point conditions KRATOS_REGISTER_CONDITION("PointCondition2D1N", mPointCondition2D1N); KRATOS_REGISTER_CONDITION("PointCondition3D1N", mPointCondition3D1N); @@ -199,6 +201,7 @@ void KratosApplication::RegisterKratosCore() { //Master-slave constraints KRATOS_REGISTER_CONSTRAINT("MasterSlaveConstraint",mMasterSlaveConstraint); KRATOS_REGISTER_CONSTRAINT("LinearMasterSlaveConstraint",mLinearMasterSlaveConstraint); + KRATOS_REGISTER_CONSTRAINT("LinearMultifreedomConstraint", mLinearMultifreedomConstraint); KRATOS_REGISTER_CONDITION("PeriodicCondition", mPeriodicCondition) KRATOS_REGISTER_CONDITION("PeriodicConditionEdge", mPeriodicConditionEdge)