diff --git a/.clang-tidy b/.clang-tidy index a2ced517b..045a9f667 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -22,6 +22,16 @@ CheckOptions: - key: readability-function-cognitive-complexity.IgnoreMacros value: 'true' +# Ignore bool conversion + - key: readability-implicit-bool-conversion.AllowIntegerConditions + value: 'true' + - key: readability-implicit-bool-conversion.AllowPointerConditions + value: 'true' + +# Structs are meant to be public + - key: misc-non-private-member-variables-in-classes.IgnoreClassesWithAllMemberVariablesBeingPublic + value: 'true' + # Naming conventions - key: readability-identifier-naming.ClassCase value: CamelCase diff --git a/include/prismspf/core/constraint_manager.h b/include/prismspf/core/constraint_manager.h new file mode 100644 index 000000000..29fb249eb --- /dev/null +++ b/include/prismspf/core/constraint_manager.h @@ -0,0 +1,222 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include "prismspf/core/types.h" + +PRISMS_PF_BEGIN_NAMESPACE + +// TODO (fractalsbyx): The following snippet is from dealii. Using this (by +// pre-constructing the needed maps and functions) may be cleaner than how dirichlet are +// currently done. +/* template + void + interpolate_boundary_values( + const Mapping &mapping, + const DoFHandler &dof, + const std::map *> + &function_map, + AffineConstraints &constraints, + const ComponentMask &component_mask = {}); */ + +// TODO (fractalsbyx): This class seems to parallel DofManager quite a bit. Consider +// merging them? +/** + * @brief The class handles the generation and application of boundary conditions based on + * the user-inputs. + */ +template +class ConstraintManager +{ +public: + /** + * @brief Constructor. + */ + ConstraintManager(const std::vector &field_attributes, + const std::set &solve_groups, + const DofManager &_dof_manager, + const PDEOperator &_pde_operator); + + /** + * @brief Getter function for the constraints. + */ + [[nodiscard]] std::vector *> + get_constraints(const std::set &field_indices, + unsigned int relative_level = 0) const; + + /** + * @brief Getter function for the constraint of an index (constant reference). + */ + [[nodiscard]] const dealii::AffineConstraints & + get_constraint(Types::Index index, unsigned int relative_level = 0) const; + + /** + * @brief Getter function for the constraints. + */ + [[nodiscard]] std::vector *> + get_change_constraints(const std::set &field_indices, + unsigned int relative_level = 0) const; + + /** + * @brief Getter function for the constraint of an index (constant reference). + */ + [[nodiscard]] const dealii::AffineConstraints & + get_change_constraint(Types::Index index, unsigned int relative_level = 0) const; + + /** + * @brief Make constraints based on the inputs of the constructor. + */ + void + reinit(const dealii::Mapping &mapping); + + /** + * @brief Update time-dependent constraints. + * For now this only updates the non-uniform dirichlet constraints. + */ + void + update_time_dependent_constraints( + const dealii::Mapping &mapping, + const std::vector *> &dof_handlers); + +private: + /** + * @brief Create a component mask. + */ + static const std::array vector_component_mask; + static const dealii::ComponentMask scalar_empty_mask; + + /** + * @brief Add boundary conditions to a single constraint. + */ + void + make_bc_constraints(dealii::AffineConstraints &constraint, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + const std::map &boundary_condition, + FieldInfo::TensorRank tensor_rank, + Types::Index field_index, + bool for_change_term = false); + + /** + * @brief Apply constraints for common boundary conditions. + */ + void + make_one_boundary_constraint(dealii::AffineConstraints &_constraints, + unsigned int boundary_id, + unsigned int component, + BoundaryCondition::Type boundary_type, + number dirichlet_value, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + FieldInfo::TensorRank tensor_rank, + Types::Index field_index, + bool for_change_term = false) const; + + /** + * @brief Apply natural constraints. + */ + void + make_natural_constraints() const; + + /** + * @brief make dirichlet constraints. + */ + void + make_uniform_dirichlet_constraints(dealii::AffineConstraints &_constraints, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + const unsigned int &boundary_id, + const bool &is_vector_field, + const number &value, + + const dealii::ComponentMask &mask) const; + + /** + * @brief make nonuniform dirichlet constraints. + */ + void + make_nonuniform_dirichlet_constraints(dealii::AffineConstraints &_constraints, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + const unsigned int &boundary_id, + const unsigned int &field_index, + const bool &is_vector_field, + const dealii::ComponentMask &mask, + bool is_change_term = false) const; + + /** + * @brief make periodic constraints. + */ + void + make_periodic_constraints(dealii::AffineConstraints &_constraints, + const dealii::DoFHandler &dof_handler, + const unsigned int &boundary_id, + const dealii::ComponentMask &mask) const; + + /** + * @brief Set the dirichlet constraint for the pinned point. + */ + void + set_pinned_point(dealii::AffineConstraints &constraint, + const dealii::Point &target_point, + const std::array &value, + const dealii::DoFHandler &dof_handler, + FieldInfo::TensorRank tensor_rank, + bool is_change_term = false) const; + + /** + * @brief User-inputs. + */ + const UserInputParameters *user_inputs; + + /** + * @brief Dof manager pointer. + */ + std::shared_ptr> dof_manager; + + /** + * @brief PDE operator number. + */ + std::shared_ptr> pde_operator; + + /** + * @brief Whether we have multigrid. + */ + bool has_multigrid = false; + + /** + * @brief Global minimum level for multigrid. + */ + unsigned int global_min_level = 0; + + /** + * @brief Constraints. Outer vector is indexed by field index. Inner vector is indexed + * by relative mg level. + */ + std::vector>>> constraints; + /** + * @brief Constraints for Newton-Change solutions. Outer vector is indexed by field + * index. Inner vector is indexed by relative mg level. + */ + std::vector>>> + change_constraints; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/dependencies.h b/include/prismspf/core/dependencies.h new file mode 100644 index 000000000..4a27f21b0 --- /dev/null +++ b/include/prismspf/core/dependencies.h @@ -0,0 +1,213 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +namespace Dependencies +{ + static const std::array nothing_arr = []() + { + std::array arr {}; + arr.fill(EvalFlags::nothing); + return arr; + }(); + + // NOLINTBEGIN(misc-non-private-member-variables-in-classes, hicpp-explicit-conversions) + struct Dependency + { + using EvalFlags = dealii::EvaluationFlags::EvaluationFlags; + EvalFlags flag = EvalFlags::nothing; + EvalFlags change_flag = EvalFlags::nothing; + + std::array old_flags {}; + + // Intentional implicit conversions + /** + * @brief Construct a Dependency with given flags. + */ + Dependency( + EvalFlags _flag = EvalFlags::nothing, + EvalFlags _change_flag = EvalFlags::nothing, + std::array _old_flags = nothing_arr) + : flag(_flag) + , change_flag(_change_flag) + , old_flags(_old_flags) + {} + + Dependency + operator|(const Dependency &other) const + { + Dependency result(static_cast(static_cast(flag) | + static_cast(other.flag)), + static_cast( + static_cast(change_flag) | + static_cast(other.change_flag))); + for (unsigned int i = 0; i < Numbers::max_saved_increments; ++i) + { + result.old_flags.at(i) = + static_cast(static_cast(old_flags.at(i)) | + static_cast(other.old_flags.at(i))); + } + return result; + } + + Dependency + operator&(const Dependency &other) const + { + Dependency result(static_cast(static_cast(flag) & + static_cast(other.flag)), + static_cast( + static_cast(change_flag) & + static_cast(other.change_flag))); + for (unsigned int i = 0; i < Numbers::max_saved_increments; ++i) + { + result.old_flags.at(i) = + static_cast(static_cast(old_flags.at(i)) & + static_cast(other.old_flags.at(i))); + } + return result; + } + + Dependency & + operator|=(const Dependency &other) + { + flag = static_cast(static_cast(flag) | + static_cast(other.flag)); + change_flag = static_cast(static_cast(change_flag) | + static_cast(other.change_flag)); + for (unsigned int i = 0; i < Numbers::max_saved_increments; ++i) + { + old_flags.at(i) = + static_cast(static_cast(old_flags.at(i)) | + static_cast(other.old_flags.at(i))); + } + + return *this; + } + + Dependency & + operator&=(const Dependency &other) + { + flag = static_cast(static_cast(flag) & + static_cast(other.flag)); + change_flag = static_cast(static_cast(change_flag) & + static_cast(other.change_flag)); + for (unsigned int i = 0; i < Numbers::max_saved_increments; ++i) + { + old_flags.at(i) = + static_cast(static_cast(old_flags.at(i)) & + static_cast(other.old_flags.at(i))); + } + + return *this; + } + }; + + // NOLINTEND(misc-non-private-member-variables-in-classes, hicpp-explicit-conversions) + + static const Dependency value {EvalFlags::values}; + static const Dependency gradient {EvalFlags::gradients}; + static const Dependency value_and_gradient {EvalFlags::values | EvalFlags::gradients}; +} // namespace Dependencies + +using Dependencies::Dependency; + +using DependencySet = std::map; + +DependencySet +make_dependency_set(const std::vector &field_attributes, + std::set dependency_strings) +{ + static const std::set, EvalFlags>> + reg_delimiters = { + {{"", ""}, EvalFlags::values }, + {{"grad(", ")"}, EvalFlags::gradients}, + {{"hess(", ")"}, EvalFlags::hessians }, + {{"hessdiag(", ")"}, EvalFlags::hessians }, + {{"lap(", ")"}, EvalFlags::hessians }, + {{"div(", ")"}, EvalFlags::gradients}, + {{"symgrad(", ")"}, EvalFlags::gradients}, + {{"curl(", ")"}, EvalFlags::gradients}, + }; + + DependencySet result; + for (unsigned int i = 0; i < field_attributes.size(); ++i) + { + const auto &attr = field_attributes[i]; + for (const auto &[delimiter, flag] : reg_delimiters) + { + std::string potential_match = delimiter.first + attr.name + delimiter.second; + auto iter = dependency_strings.find(potential_match); + if (iter != dependency_strings.end()) + { + result[i].flag |= flag; + dependency_strings.erase(iter); + } + potential_match = + delimiter.first + "change(" + attr.name + ")" + delimiter.second; + iter = dependency_strings.find(potential_match); + if (iter != dependency_strings.end()) + { + result[i].change_flag |= flag; + dependency_strings.erase(iter); + } + for (unsigned int old_index = 0; old_index < Numbers::max_saved_increments; + ++old_index) + { + potential_match = delimiter.first + "old_" + std::to_string(old_index + 1) + + "(" + attr.name + ")" + delimiter.second; + iter = dependency_strings.find(potential_match); + if (iter != dependency_strings.end()) + { + result[i].old_flags.at(old_index) |= flag; + dependency_strings.erase(iter); + } + } + } + } + if (!dependency_strings.empty()) + { + std::string error_message = "The following dependencies were not recognized: "; + for (const auto &dep : dependency_strings) + { + error_message += dep + ", "; + } + error_message.pop_back(); + error_message.pop_back(); +#ifndef DEBUG + ConditionalOStreams::pout_base() << "Warning: " << error_message << std::endl; +#endif + Assert(false, dealii::ExcMessage(error_message)); + } + return result; +} + +DependencySet +make_dependency_set(const std::vector &field_attributes, + const std::string &dependency_string) +{ + auto dependency_strings = dealii::Utilities::split_string_list(dependency_string); + return make_dependency_set(field_attributes, + std::set(dependency_strings.begin(), + dependency_strings.end())); +} + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/include/prismspf/core/dependency_extents.h b/include/prismspf/core/dependency_extents.h new file mode 100644 index 000000000..0c0909c54 --- /dev/null +++ b/include/prismspf/core/dependency_extents.h @@ -0,0 +1,55 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +struct DependencyExtent +{ + Dependency dependency; + unsigned int mg_depth = 0; + + [[nodiscard]] unsigned int + age() const + { + for (unsigned int i = Numbers::max_saved_increments - 1; i >= 0; --i) + { + if (dependency.old_flags.at(i) != dealii::EvaluationFlags::nothing) + { + return i; + } + } + return 0; + } + + static std::vector + calculate_extents(const std::vector &attributes_list, + const std::set &solve_groups) + { + std::vector extents(attributes_list.size()); + for (const auto &solve_group : solve_groups) + { + for (const auto &field_index : solve_group.field_indices) + { + DependencyExtent &extent = extents[field_index]; + for (const auto &[index, dep] : solve_group.dependencies_rhs) + { + extent.dependency |= dep; + extent.mg_depth = std::max(extent.mg_depth, solve_group.mg_depth); + } + } + } + return extents; + }; +}; + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/include/prismspf/core/dof_manager.h b/include/prismspf/core/dof_manager.h new file mode 100644 index 000000000..a0901e483 --- /dev/null +++ b/include/prismspf/core/dof_manager.h @@ -0,0 +1,112 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief Class that manages the deal.II DoFHandlers + */ +template +class DofManager +{ +public: + /** + * @brief Constructor. + */ + explicit DofManager(const std::vector &field_attributes); + + /** + * @brief Initialize the DoFHandlers + */ + void + init(const TriangulationManager &triangulation_handler); + + /** + * @brief Reinitialize the DoFHandlers + */ + void + reinit(const TriangulationManager &triangulation_handler); + + /** + * @brief Getter function for the DoFHandlers (constant reference). + */ + [[nodiscard]] const std::vector *> & + get_field_dof_handlers(const std::set &field_indices, + unsigned int relative_level = 0) const; + + /** + * @brief Getter function for the DoFHandler (reference). + */ + [[nodiscard]] const dealii::DoFHandler & + get_field_dof_handler(Types::Index field_index, unsigned int relative_level = 0) const + { + return *dof_handlers[field_index][relative_level]; + } + + /** + * @brief Getter function for the scalar and vector DoFHandlers. + */ + const std::vector, 2>> & + get_dof_handlers() const + { + return level_dof_handlers; + } + + /** + * @brief Getter function for a specific scalar or vector DoFHandler. + */ + [[nodiscard]] const dealii::DoFHandler & + get_dof_handler(const unsigned int &rank, unsigned int relative_level = 0) const + { + return level_dof_handlers[relative_level][rank]; + } + + /** + * @brief Get the total DoFs excluding multigrid DoFs. + */ + [[nodiscard]] dealii::types::global_dof_index + get_total_dofs() const + { + dealii::types::global_dof_index n_dofs = 0; + for (const auto &dof_handler_set : dof_handlers) + { + n_dofs += dof_handler_set[0]->n_dofs(); + } + return n_dofs; + } + +private: + /** + * @brief Collection of the triangulation DoFs. The number of DoFHandlers should be + * equal to or less than the number of fields. Technically, there's a small + * optimization we can use when multiple fields have the same constraints and + * quadrature rule, allowing us to share the same DoFHandler. An example of this might + * be grain growth. + * Outer vector is indexed by field index. Inner vector is indexed by relative mg level. + */ + std::vector>>> dof_handlers; + + /** + * @brief A scalar and a vector dof handler for each level + */ + std::vector, 2>> level_dof_handlers; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/dst_container.h b/include/prismspf/core/dst_container.h new file mode 100644 index 000000000..6f38ef87d --- /dev/null +++ b/include/prismspf/core/dst_container.h @@ -0,0 +1,202 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "prismspf/core/field_attributes.h" +#include "prismspf/core/solve_group.h" +#include "prismspf/core/variable_attributes.h" + +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class ElementVolume; + +/** + * @brief This class permits the access of a subset of indexed fields and gives an error + * if any non-allowed fields are requested. + * + * @tparam dim The number of dimensions in the problem. + * @tparam degree The polynomial degree of the shape functions. + * @tparam number Datatype to use for `dealii::VectorizedArray`. Either + * double or float. + */ +template +class DSTContainer +{ +public: + using TensorRank = FieldInfo::TensorRank; + /** + * @brief Typedef for the basic vector that apply our operations to. + */ + using SolutionVector = SolutionIndexer::SolutionVector; + using MatrixFree = dealii::MatrixFree>; + + /** + * @brief Typedef for the basic value that the user manipulates. + */ + using ScalarValue = dealii::VectorizedArray; + using VectorValue = dealii::Tensor<1, dim, ScalarValue>; + using ScalarGradient = VectorValue; + using VectorGradient = dealii::Tensor<2, dim, ScalarValue>; + /** + * @brief Typedef for scalar evaluation objects. + */ + using ScalarFEEvaluation = dealii:: + FEEvaluation>; + /** + * @brief Typedef for vector evaluation objects. + */ + using VectorFEEvaluation = dealii:: + FEEvaluation>; + + template + using Value = std::conditional_t>; + template + using Gradient = dealii::Tensor; + template + using Hessian = dealii::Tensor; + template + using FEEval = + dealii::FEEvaluation::n_independent_components, + number, + ScalarValue>; + + template + struct GetRankHelper + { + static constexpr TensorRank rank_from_val = TensorRank(Type::rank); + static constexpr TensorRank rank_from_grad = TensorRank(Type::rank - 1); + }; + + template <> + struct GetRankHelper + { + static constexpr TensorRank rank_from_val = TensorRank::Scalar; + // TODO: make sure the enable_if is correct syntax + template > + static constexpr TensorRank rank_from_grad = TensorRank::Scalar; + }; + + template + using GetRankFromVal = GetRankHelper::rank_from_val; + template + using GetRankFromGrad = GetRankHelper::rank_from_grad; + + template + std::vector> + get_relevant_feeval_vector() + { + if constexpr (Rank == TensorRank::Scalar) + { + return feeval_scalar; + } + else if constexpr (Rank == TensorRank::Vector) + { + return feeval_vector; + } + } + + /** + * @brief Constructor. + */ + DSTContainer(const std::set &_field_indices, + const std::vector &_field_attributes, + const MatrixFree &matrix_free, + const std::vector &_field_to_block_index); + + /** + * @brief Initialize based on cell for all dependencies. + */ + void + reinit(unsigned int cell); + + /** + * @brief Set the current quadtrature point + */ + void + set_q_point(unsigned int quad) + { + q_point = quad; + } + + /** + * @brief Set the residual value of the specified scalar/vector field. + */ + template + void + set_value_term(Types::Index global_variable_index, const ValType &val) + { + get_relevant_feeval_vector>()[global_variable_index] + .submit_value(val, q_point); + integration_flags[global_variable_index] |= EvalFlags::values; + } + + /** + * @brief Set the residual gradient of the specified scalar/vector field. + */ + template + void + set_gradient_term(Types::Index global_variable_index, const GradType &val) + { + get_relevant_feeval_vector>()[global_variable_index] + .submit_gradient(val, q_point); + integration_flags[global_variable_index] |= EvalFlags::gradients; + } + + /** + * @brief Check that a value is valid for submission. + */ + void + submission_valid(Types::Index field_index, DependencyType dependency_type) const; + +private: + /** + * @brief Integrate a single field. + */ + void + integrate(unsigned int field_index); + + /** + * @brief Integrate a single field and update the solution vector. + */ + void + integrate_and_distribute(unsigned int field_index, SolutionVector &dst); + + const std::vector *field_attributes_ptr; + const std::set *field_indices; + std::vector> feeval_scalar; + std::vector> feeval_vector; + std::vector integration_flags; + + /** + * @brief The quadrature point index. + */ + unsigned int q_point = 0; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/field_attributes.h b/include/prismspf/core/field_attributes.h new file mode 100644 index 000000000..66ebd1cfe --- /dev/null +++ b/include/prismspf/core/field_attributes.h @@ -0,0 +1,120 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE +using EvalFlags = dealii::EvaluationFlags::EvaluationFlags; +using Rank = FieldInfo::TensorRank; + +// NOLINTBEGIN(misc-non-private-member-variables-in-classes) +/** + * @brief Structure to hold the attributes of a field. This includes things like + * the name, rank, and nucleation information.. + */ +struct FieldAttributes +{ + /** + * @brief Constructor + */ + explicit FieldAttributes( + std::string _name = "", + Rank _field_type = Rank::Undefined, + EvalFlags _eval_flags_rhs = EvalFlags::nothing, + EvalFlags _eval_flags_lhs = EvalFlags::nothing, + bool _is_nucleation_rate_variable = false, + std::vector _nucleating_field_indices = std::vector()) + : name(std::move(_name)) + , field_type(_field_type) + , eval_flags_rhs(_eval_flags_rhs) + , eval_flags_lhs(_eval_flags_lhs) + , is_nucleation_rate_variable(_is_nucleation_rate_variable) + , nucleating_field_indices(std::move(_nucleating_field_indices)) + {} + + /** + * @brief Field name. + */ + std::string name; + + /** + * @brief Field type (Scalar/Vector). + */ + Rank field_type = Rank::Undefined; + + /** + * @brief Evaluation flags for the types of residual the user is expected to submit to + * on the RHS. + */ + EvalFlags eval_flags_rhs = EvalFlags::nothing; + + /** + * @brief Evaluation flags for the types of residual the user is expected to submit to + * on the LHS. This is empty for Explicit fields. + */ + EvalFlags eval_flags_lhs = EvalFlags::nothing; + + /** + * @brief Is a nucleation rate + */ + bool is_nucleation_rate_variable = false; + + /** + * @brief If this is a nucleation rate, the indices of the nucleating fields + */ + std::vector nucleating_field_indices; +}; + +// NOLINTEND(misc-non-private-member-variables-in-classes) + +// TODO: Consider making these static members instead of prismspf:: + +/** + * @brief Make a map that maps field names to field indices. + */ +std::map +field_index_map(const std::vector &fields) +{ + std::map map; + for (unsigned int i = 0; i < fields.size(); ++i) + { + AssertThrow(map.find(fields[i].name) == map.end(), + "The names of the fields are not unique. This is not allowed."); + map[fields[i].name] = i; + } + return map; +} + +/** + * @brief Make a map that maps field names to field attributes. + */ +std::map +field_map(const std::vector &fields) +{ + std::map map; + for (const FieldAttributes &field : fields) + { + AssertThrow(map.find(field.name) == map.end(), + "The names of the fields are not unique. This is not allowed."); + map[field.name] = field; + } + return map; +} + +// TODO: Submit a PR/issue to dealii to make operator| constexpr. +// constexpr EvalFlags values_and_gradients = EvalFlags::values | EvalFlags::gradients; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/field_container.h b/include/prismspf/core/field_container.h new file mode 100644 index 000000000..0ae3f3cd1 --- /dev/null +++ b/include/prismspf/core/field_container.h @@ -0,0 +1,558 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "prismspf/core/dependencies.h" +#include "prismspf/core/field_attributes.h" +#include "prismspf/core/solve_group.h" +#include "prismspf/core/variable_attributes.h" + +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class ElementVolume; + +/** + * @brief This class permits the access of a subset of indexed fields and gives an error + * if any non-allowed fields are requested. + * + * @tparam dim The number of dimensions in the problem. + * @tparam degree The polynomial degree of the shape functions. + * @tparam number Datatype to use for `dealii::VectorizedArray`. Either + * double or float. + */ +template +class FieldContainer +{ +public: + using TensorRank = FieldInfo::TensorRank; + /** + * @brief Typedef for the basic vector that apply our operations to. + */ + using SolutionVector = SolutionIndexer::SolutionVector; + using MatrixFree = dealii::MatrixFree>; + + /** + * @brief Typedef for the basic value that the use manipulates. + */ + using ScalarValue = dealii::VectorizedArray; + using VectorValue = dealii::Tensor<1, dim, ScalarValue>; + using ScalarGradient = VectorValue; + using VectorGradient = dealii::Tensor<2, dim, ScalarValue>; + /** + * @brief Typedef for scalar evaluation objects. + */ + using ScalarFEEvaluation = dealii:: + FEEvaluation>; + /** + * @brief Typedef for vector evaluation objects. + */ + using VectorFEEvaluation = dealii:: + FEEvaluation>; + + template + using Value = std::conditional_t>; + template + using Gradient = dealii::Tensor; + template + using Hessian = dealii::Tensor; + template + using FEEval = + dealii::FEEvaluation::n_independent_components, + number, + ScalarValue>; + + template + struct GetRankHelper + { + static constexpr TensorRank rank_from_val = TensorRank(Type::rank); + static constexpr TensorRank rank_from_grad = TensorRank(Type::rank - 1); + }; + + template <> + struct GetRankHelper + { + static constexpr TensorRank rank_from_val = TensorRank::Scalar; + // TODO: make sure the enable_if is correct syntax + template > + static constexpr TensorRank rank_from_grad = TensorRank::Scalar; + }; + + template + using GetRankFromVal = GetRankHelper::rank_from_val; + template + using GetRankFromGrad = GetRankHelper::rank_from_grad; + + template + std::vector> + get_relevant_feeval_vector() + { + if constexpr (Rank == TensorRank::Scalar) + { + return feeval_deps_scalar; + } + else if constexpr (Rank == TensorRank::Vector) + { + return feeval_deps_vector; + } + } + + /** + * @brief Typedef for scalar diagonal matrix objects. + */ + using ScalarDiagonal = dealii::AlignedVector; + + /** + * @brief Typedef for vector diagonal matrix objects. + */ + using VectorDiagonal = dealii::AlignedVector; + + /** + * @brief Constructor. + */ + FieldContainer(const std::vector &_field_attributes, + SolutionIndexer &_solution_indexer, + unsigned int _relative_level, + const DependencySet &dependency_map); + + /** + * @brief Initialize based on cell for all dependencies. + */ + void + reinit(unsigned int cell); + + /** + * @brief Read solution vector, and evaluate based on + * dependency flags for all dependencies. + */ + void + eval(); + + /** + * @brief Initialize based on cell, read solution vector, and evaluate based on + * dependency flags for all dependencies. + */ + void + reinit_and_eval(unsigned int cell); + + /** + * @brief Set the current quadtrature point + */ + void + set_q_point(unsigned int quad) + { + q_point = quad; + } + + /** + * @brief Return the value of the specified field. + * + * @tparam T the return type. Must be either a `ScalarValue` or `dealii::Tensor<1, dim, + * ScalarValue>`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + template + [[nodiscard]] Value + get_value(Types::Index global_variable_index) const + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_value(q_point); + } + + /** + * @brief Return the gradient of the specified field. + * + * @tparam T the return type. Must be either a `dealii::Tensor<1, dim, + * ScalarValue>` or `dealii::Tensor<2, dim, ScalarValue>`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + template + [[nodiscard]] Gradient + get_gradient(Types::Index global_variable_index) const + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_gradient(q_point); + } + + /** + * @brief Return the hessian of the specified field. + * + * @tparam T the return type. Must be either a `dealii::Tensor<2, dim, + * ScalarValue>` or `dealii::Tensor<3, dim, ScalarValue>`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + template + [[nodiscard]] Hessian + get_hessian(Types::Index global_variable_index) const + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_hessian(q_point); + } + + /** + * @brief Return the diagonal of the hessian of the specified field. + * + * @tparam T the return type. Must be either a `dealii::Tensor<1, dim, + * ScalarValue>` or `dealii::Tensor<2, dim, ScalarValue>`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + template + [[nodiscard]] Gradient + get_hessian_diagonal(Types::Index global_variable_index) const + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_hessian_diagonal(q_point); + } + + /** + * @brief Return the laplacian of the specified field. + * + * @tparam T the return type. Must be either a `ScalarValue` or `dealii::Tensor<1, dim, + * ScalarValue>`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + template + [[nodiscard]] Value + get_laplacian(Types::Index global_variable_index) const + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_laplacian(q_point); + } + + /** + * @brief Return the divergence of the specified field. + * + * @tparam T the return type. Must be `ScalarValue`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + // TODO: FIgure out the assertion here. Dealii probly will have one, but we should too. + template + [[nodiscard]] auto + get_divergence(Types::Index global_variable_index) const + -> ScalarValue /* Value::value_type */ + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_divergence(q_point); + } + + /** + * @brief Return the symmetric gradient of the specified field. + * + * @tparam T the return type. Must be either a `dealii::Tensor<2, dim, + * ScalarValue>` or `dealii::SymmetricTensor<2, dim, ScalarValue>`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + // TODO: FIgure out the assertion here. Dealii probly will have one, but we should too. + template + [[nodiscard]] auto + get_symmetric_gradient(Types::Index global_variable_index) const + -> dealii::SymmetricTensor<2, dim, ScalarValue> + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_symmetric_gradient(q_point); + } + + /** + * @brief Return the curl of the specified field. + * + * Note that this is dealii::VectorizedArray type for 2D and dealii::Tensor<1, + * dim, dealii::VectorizedArray> type for 3D. + * + * @tparam T the return type. Must be either a `dealii::Tensor<1, 1, + * ScalarValue>` or `dealii::Tensor<1, dim, ScalarValue>`. + * @param global_variable_index The global index of the variable to access. + * @param dependency_type The dependency type of the variable to access. + */ + // TODO: FIgure out the assertion here. Dealii probly will have one, but we should too. + template + [[nodiscard]] auto + get_curl(Types::Index global_variable_index) const + -> dealii::Tensor<1, (dim == 2 ? 1 : dim), ScalarValue> + { + return get_relevant_feeval_vector()[global_variable_index] + .template get() + .get_curl(q_point); + } + + /** + * @brief Set the residual value of the specified scalar/vector field. + */ + template + void + set_value_term(Types::Index global_variable_index, const ValType &val) + { + get_relevant_feeval_vector>()[global_variable_index] + .template get() + .submit_value(val, q_point); + } + + /** + * @brief Set the residual gradient of the specified scalar/vector field. + */ + template + void + set_gradient_term(Types::Index global_variable_index, const GradType &val) + { + get_relevant_feeval_vector>()[global_variable_index] + .template get() + .submit_gradient(val, q_point); + } + +private: + /** + * @brief Check whether the entry for the FEEvaluation is within the bounds of the + * vector. + */ + void + feevaluation_size_valid(Types::Index field_index) const; + + /** + * @brief Check whether the entry for the FEEvaluation is within the bounds of the + * vector and not a nullptr. + */ + void + feevaluation_exists(Types::Index field_index, Types::Index dependency_index) const; + + /** + * @brief Check that a variable value/gradient/hessians was marked as needed and thus + * properly initialized. + */ + void + access_valid(Types::Index field_index, + DependencyType dependency_type, + dealii::EvaluationFlags::EvaluationFlags flag) const; + + /** + * @brief Check that a value is valid for submission. + */ + void + submission_valid(Types::Index field_index, DependencyType dependency_type) const; + + /** + * @brief Return the number of quadrature points. + */ + [[nodiscard]] unsigned int + get_n_q_points() const; + + /** + * @brief Return the quadrature point location. + */ + [[nodiscard]] dealii::Point + get_q_point_location() const; + + /** + * @brief Integrate the residuals. + */ + void + integrate(); + + /** + * @brief Integrate the residuals and distribute from local to global. + */ + void + integrate_and_distribute(); + + /** + * @brief Struct to hold feevaluation relevant for this solve. + */ + template + struct FEEValuationDeps + { + using FEEDepPairPtr = std::unique_ptr>; + FEEDepPairPtr fe_eval; + FEEDepPairPtr fe_eval_change; + std::array fe_eval_old; + /** + * @brief The solution group and the block index + * @note It would look nicer to just use the SolutionIndexer, but this way decreases + * indexing + */ + const SolutionLevel *solution_level = nullptr; + unsigned int block_index = -1; + + FEEValuationDeps( + const Dependencies::Dependency &dependency, + std::pair *, unsigned int> mf_id_pair) + : solution_level(mf_id_pair.first) + , block_index(mf_id_pair.second) + { + if (dependency.flag) + { + fe_eval = + std::make_unique(solution_level->matrix_free, block_index); + } + if (dependency.change_flag) + { + fe_eval_change = + std::make_unique(solution_level->matrix_free, block_index); + } + for (unsigned int age = 0; age < Numbers::max_saved_increments; ++age) + { + if (dependency.old_flags.at(age)) + { + fe_eval_old[age] = + std::make_unique(solution_level->matrix_free, + block_index); + } + } + } + + template + FEEvaluationType & + get() + { + if constexpr (type == DependencyType::Normal) + { + return *fe_eval; + } + else if constexpr (type == DependencyType::Change) + { + return *fe_eval_change; + } + else + { + return *fe_eval_old[int(type)]; + } + } + + void + reinit(unsigned int cell) + { + if (fe_eval) + { + fe_eval->first.reinit(cell); + } + if (fe_eval_change) + { + fe_eval_change->first.reinit(cell); + } + for (auto &old_fe_eval : fe_eval_old) + { + if (old_fe_eval) + { + old_fe_eval->first.reinit(cell); + } + } + } + + void + eval() + { + if (fe_eval) + { + fe_eval->first.read_dof_values_plain( + solution_level->solutions.block(block_index)); + fe_eval->first.evaluate(fe_eval->second); + } + if (fe_eval_change) + { + fe_eval_change->first.read_dof_values_plain( + solution_level->change_solutions.block(block_index)); + fe_eval_change->first.evaluate(fe_eval_change->second); + } + for (unsigned int age = 0; age < Numbers::max_saved_increments; ++age) + { + if (FEEDepPairPtr &old_fe_eval = fe_eval_old[age]) + { + old_fe_eval->first.read_dof_values_plain( + solution_level->old_solutions[age].block(block_index)); + old_fe_eval->first.evaluate(old_fe_eval->second); + } + } + } + + void + reinit_and_eval(unsigned int cell) + { + if (fe_eval) + { + fe_eval->first.reinit(cell); + fe_eval->first.read_dof_values_plain( + solution_level->solutions.block(block_index)); + fe_eval->first.evaluate(fe_eval->second); + } + if (fe_eval_change) + { + fe_eval_change->first.reinit(cell); + fe_eval_change->first.read_dof_values_plain( + solution_level->change_solutions.block(block_index)); + fe_eval_change->first.evaluate(fe_eval_change->second); + } + for (unsigned int age = 0; age < Numbers::max_saved_increments; ++age) + { + if (FEEDepPairPtr &old_fe_eval = fe_eval_old[age]) + { + old_fe_eval->first.reinit(cell); + old_fe_eval->first.read_dof_values_plain( + solution_level->old_solutions[age].block(block_index)); + old_fe_eval->first.evaluate(old_fe_eval->second); + } + } + } + }; + + const std::vector *field_attributes_ptr; + const SolveGroup *solve_group; + const SolutionIndexer *solution_indexer; + unsigned int relative_level; + std::vector> feeval_deps_scalar; + std::vector> feeval_deps_vector; + std::vector> dst_feeval_scalar; + std::vector> dst_feeval_vector; + + /** + * @brief The element volume container + */ + const ElementVolume *element_volume_handler; + + /** + * @brief The quadrature point index. + */ + unsigned int q_point = 0; + + /** + * @brief Number of DoFs per cell. + */ + unsigned int n_dofs_per_cell = 0; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/grid_refiner.h b/include/prismspf/core/grid_refiner.h index 3158196ec..04f2a1c96 100644 --- a/include/prismspf/core/grid_refiner.h +++ b/include/prismspf/core/grid_refiner.h @@ -445,11 +445,12 @@ class GridRefiner marker_functions.begin(), marker_functions.end(), [&](const std::shared_ptr> &marker_function) - { - return marker_function->flag(*cell, - grid_refinement_context.get_user_inputs() - .get_temporal_discretization()); - })) + { + return marker_function->flag(*cell, + grid_refinement_context + .get_user_inputs() + .get_temporal_discretization()); + })) { cell->set_user_flag(); cell->clear_coarsen_flag(); diff --git a/include/prismspf/core/group_solution_handler.h b/include/prismspf/core/group_solution_handler.h new file mode 100644 index 000000000..5abd81828 --- /dev/null +++ b/include/prismspf/core/group_solution_handler.h @@ -0,0 +1,282 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief The solution vectors. + */ +template +struct SolutionLevel +{ + using BlockVector = dealii::LinearAlgebra::distributed::BlockVector; + using SolutionVector = BlockVector::BlockType; + using MatrixFree = dealii::MatrixFree>; + BlockVector solutions; + std::array old_solutions; + BlockVector change_solutions; + MatrixFree matrix_free; +}; + +/** + * @brief Class that manages solution initialization and swapping with old solutions. + */ +template +class GroupSolutionHandler +{ +public: + using BlockVector = SolutionLevel::BlockVector; + using SolutionVector = SolutionLevel::SolutionVector; + using MatrixFree = SolutionLevel::MatrixFree; +#if DEAL_II_VERSION_MAJOR >= 9 && DEAL_II_VERSION_MINOR >= 7 + using SolutionTransfer = dealii::SolutionTransfer; +#else + using SolutionTransfer = + dealii::parallel::distributed::SolutionTransfer; +#endif + + /** + * @brief Constructor. + */ + GroupSolutionHandler(const SolveGroup &_solve_group, + const std::vector &_attributes_list); + + // TODO (fractalsbyx): add 'const T& get() const' versions. + /** + * @brief Get the solution vector set. This contains all the normal fields and is + * typically used for output. + */ + [[nodiscard]] BlockVector & + get_solution_full_vector(unsigned int relative_level = 0); + + /** + * @brief Get the const solution vector set. This contains all the normal fields and is + * typically used for output. + */ + [[nodiscard]] const BlockVector & + get_solution_full_vector(unsigned int relative_level = 0) const; + + /** + * @brief Get a solution vector of a given field index. + */ + [[nodiscard]] SolutionVector & + get_solution_vector(unsigned int global_index, unsigned int relative_level = 0); + + /** + * @brief Get a solution vector of a given field index. + */ + [[nodiscard]] const SolutionVector & + get_solution_vector(unsigned int global_index, unsigned int relative_level = 0) const; + + /** + * @brief Get the old solution vector set at a given age. + */ + [[nodiscard]] BlockVector & + get_old_solution_full_vector(unsigned int age, unsigned int relative_level = 0); + + /** + * @brief Get the old solution vector set at a given age. + */ + [[nodiscard]] const BlockVector & + get_old_solution_full_vector(unsigned int age, unsigned int relative_level = 0) const; + + /** + * @brief Get a solution vector of a given field index at a given age. + */ + [[nodiscard]] SolutionVector & + get_old_solution_vector(unsigned int age, + unsigned int global_index, + unsigned int relative_level = 0); + + /** + * @brief Get a solution vector of a given field index at a given age. + */ + [[nodiscard]] const SolutionVector & + get_old_solution_vector(unsigned int age, + unsigned int global_index, + unsigned int relative_level = 0) const; + + /** + * @brief Get the "new" solution vector set. + */ + [[nodiscard]] BlockVector & + get_new_solution_full_vector(unsigned int relative_level = 0); + + /** + * @brief Get the "new" solution vector set. + */ + [[nodiscard]] const BlockVector & + get_new_solution_full_vector(unsigned int relative_level = 0) const; + + /** + * @brief Get the "new" solution vector of a given field index. + */ + [[nodiscard]] SolutionVector & + get_new_solution_vector(unsigned int index, unsigned int relative_level = 0); + + /** + * @brief Get the "new" solution vector of a given field index. + */ + [[nodiscard]] const SolutionVector & + get_new_solution_vector(unsigned int index, unsigned int relative_level = 0) const; + + /** + * @brief Get the solutions object at a level. + */ + [[nodiscard]] SolutionLevel & + get_solution_level(unsigned int relative_level = 0); + + /** + * @brief Get the solutions object at a level. + */ + [[nodiscard]] const SolutionLevel & + get_solution_level(unsigned int relative_level = 0) const; + + /** + * @brief Get the matrix_free object at a level. + */ + [[nodiscard]] MatrixFree & + get_matrix_free(unsigned int relative_level = 0); + + /** + * @brief Get the matrix_free object at a level. + */ + [[nodiscard]] const MatrixFree & + get_matrix_free(unsigned int relative_level = 0) const; + + /** + * @brief Get the block index from the global index. + */ + [[nodiscard]] unsigned int + get_block_index(unsigned int global_index) const; + + /** + * @brief Get the underlying solve group object. + */ + [[nodiscard]] const SolveGroup & + get_solve_group() const; + + /** + * @brief Get the block index from the global index. + */ + [[nodiscard]] const std::vector & + get_global_to_block_index() const; + + /** + * @brief Initialize the solution set. + */ + template + void + init(const dealii::Mapping &mapping, + const DofManager &dof_manager, + const ConstraintManager &constraint_manager, + const dealii::Quadrature &quad); + + /** + * @brief Reinitialize the solution set. + */ + void + reinit(); + + /** + * @brief Update the ghost values. + */ + void + update_ghosts(unsigned int relative_level = 0) const; + + /** + * @brief Zero out the ghost values. + */ + void + zero_out_ghosts(unsigned int relative_level = 0) const; + + /** + * @brief Apply the given constraints to a solution vector of a given field index. + */ + void + apply_constraints(unsigned int relative_level = 0); + + /** + * @brief Apply intial condition to the old fields. For now, this simply copies the + * values in the normal field to the old. + */ + void + apply_initial_condition_for_old_fields(); + + /** + * @brief Update the `solutions` to the `new_solution` and propagate the old solutions. + */ + void + update(unsigned int relative_level = 0); + + /** + * @brief Reinit the solution transfer objections. + */ + void + init_solution_transfer(); + + /** + * @brief Prepare for solution transfer + */ + void + prepare_for_solution_transfer(); + + /** + * @brief Transfer solutions + */ + void + execute_solution_transfer(); + +private: + /** + * @brief Information about the solve group this handler is responsible for. + */ + SolveGroup solve_group; + + // TODO (fractalsbyx): do this better + /** + * @brief Oldest saved age. + */ + int oldest_saved = 0; + + /** + * @brief Mapping from block index to global field index. + */ + std::vector block_to_global_index; + + /** + * @brief Mapping from global field index to block index. + */ + std::vector global_to_block_index; + + // TODO (fractalsbyx): Consider switching to dealii::MGLevelObject + // Right now, I prefer using the relative level. + std::vector> solution_levels; + + /** + * @brief Utility for solution transfer to different mesh (for AMR). + */ + SolutionTransfer solution_transfer; + std::array old_solution_transfer; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/invm_manager.h b/include/prismspf/core/invm_manager.h new file mode 100644 index 000000000..ff43abe4f --- /dev/null +++ b/include/prismspf/core/invm_manager.h @@ -0,0 +1,180 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief A little class that computes the element volume for our triangulation. + */ +template +class InvMManager +{ +public: + using SolutionVector = GroupSolutionHandler::SolutionVector; + using MatrixFree = GroupSolutionHandler::MatrixFree; + using ScalarValue = dealii::VectorizedArray; + using VectorValue = dealii::Tensor<1, dim, ScalarValue>; + + explicit InvMManager(const DofManager &dof_manager, + bool _calculate_scalar, + bool _calculate_vector) + : calculate_scalar(_calculate_scalar) + , calculate_vector(_calculate_vector) + { + auto dof_handlers = dof_manager.get_dof_handlers(); + data.resize(dof_handlers.size()); + for (unsigned int i = 0; i < data.size(); ++i) + { + for (unsigned int rank = 0; rank < 2; ++rank) + { + data[i][rank].reinit(SystemWide::mapping, + dof_handlers[i][rank], + dealii::AffineConstraints(), // make member? + SystemWide::quadrature); + } + } + } + + /** + * @brief Initialize. + */ + void + initialize() + { + for (unsigned int i = 0; i < data.size(); ++i) + { + if (calculate_scalar) + { + data[i][0]->initialize_dof_vector(jxw_scalar); + data[i][0]->initialize_dof_vector(invm_scalar); + } + if (calculate_vector) + { + data[i][1]->initialize_dof_vector(jxw_vector); + data[i][1]->initialize_dof_vector(invm_vector); + } + } + } + + /** + * @brief Compute element volume for the triangulation + */ + void + compute_scalar_invm() + { + data->cell_loop(&InvMManager::compute_local_scalar, this, jxw_scalar, 0); + invert(invm_scalar, jxw_scalar); + } + + void + compute_vector_invm() + { + data->cell_loop(&InvMManager::compute_local_vector, this, jxw_vector, 0); + invert(invm_vector, jxw_vector); + } + + void + compute_local_scalar(const MatrixFree &_data, + SolutionVector &dst, + [[maybe_unused]] const int &src, + const std::pair &cell_range) const + { + dealii::FEEvaluation fe_eval(_data); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval->first.evaluate(dealii::update_JxW_values); + for (unsigned int quad = 0; quad < SystemWide::quadrature.size(); + ++quad) + { + fe_eval.submit_value(fe_eval.JxW(quad)); + } + fe_eval.integrate_scatter(dealii::EvaluationFlags::values, dst); + } + } + + void + compute_local_vector(const MatrixFree &_data, + SolutionVector &dst, + [[maybe_unused]] const int &src, + const std::pair &cell_range) const + { + dealii::FEEvaluation fe_eval(_data); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval->first.evaluate(dealii::update_JxW_values); + for (unsigned int quad = 0; quad < SystemWide::quadrature.size(); + ++quad) + { + fe_eval.submit_value(fe_eval.JxW(quad) * one); + } + fe_eval.integrate_scatter(dealii::EvaluationFlags::values, dst); + } + } + + void + compute_invm() + { + if (calculate_scalar) + { + compute_scalar_invm(); + } + if (calculate_vector) + { + compute_vector_invm(); + } + } + +private: + void + invert(SolutionVector &dst, const SolutionVector &src) const + { + Assert(dst.size() == src.size(), dealii::ExcInternalError()); + for (unsigned int i = 0; i < src.locally_owned_size(); ++i) + { + dst[i] = 1.0 / src[i]; + } + } + + /** + * @brief Matrix-free object. + */ + std::vector> data; + + bool calculate_scalar = false; + bool calculate_vector = false; + + /** + * @brief Vector that stores element volumes + */ + SolutionVector jxw_scalar; + SolutionVector jxw_vector; + SolutionVector invm_scalar; + SolutionVector invm_vector; + + inline static const VectorValue one = []() + { + VectorValue one1; + for (unsigned int i = 0; i < dim; ++i) + { + one1[i] = 1.0; + } + return one1; + }(); +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/problem.h b/include/prismspf/core/problem.h new file mode 100644 index 000000000..5b38e2ff0 --- /dev/null +++ b/include/prismspf/core/problem.h @@ -0,0 +1,134 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +#include "prismspf/core/field_attributes.h" + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief This is the main class that handles the construction and solving of + * user-specified PDEs. + */ +template +class Problem +{ +public: + /** + * @brief Constructor. + */ + Problem(const std::vector &field_attributes, + const std::vector &solve_groups, + const UserInputParameters &_user_inputs, + PhaseFieldTools &_pf_tools, + const std::shared_ptr> &_pde_operator); + + /** + * @brief Run initialization and solving steps of the given problem. + */ + void + run(); + +private: + /** + * @brief Main time-stepping loop that calls solve_increment, reinit_system, + * output_results, etc... + */ + void + solve(); + + /** + * @brief Solve a single increment of the given PDEs. + */ + void + solve_increment(SimulationTimer &sim_timer); + + /** + * @brief Initialize the system. + */ + void + init_system(); + + /** + * @brief Reinitialize the system. + */ + void + reinit_system(); + + /** + * @brief Field attributes. + */ + std::vector field_attributes; + + /** + * @brief Solve groups. + */ + std::vector solve_groups; + + /** + * @brief User-inputs. + */ + const UserInputParameters *user_inputs_ptr; + + /** + * @brief Phase field tools. + */ + PhaseFieldTools *pf_tools; + + /** + * @brief Triangulation handler. + */ + TriangulationManager triangulation_manager; + + /** + * @brief DoF manager. + */ + DofManager dof_manager; + + /** + * @brief Constraint handler. + */ + ConstraintManager constraint_manager; + + /** + * @brief Solvers. + */ + std::vector>> solvers; + + /** + * @brief Solution indexer + */ + SolutionIndexer solution_indexer; + + /** + * @brief Solver context. + */ + SolveContext solve_context; + + /** + * @brief Grid refiner. + */ + GridRefiner grid_refiner; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/refinement_manager.h b/include/prismspf/core/refinement_manager.h new file mode 100644 index 000000000..824d48c19 --- /dev/null +++ b/include/prismspf/core/refinement_manager.h @@ -0,0 +1,473 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include "prismspf/core/field_attributes.h" + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class RefinementManager +{ +public: + /** + * @brief Constructor. Init the flags for refinement. + */ + explicit RefinementManager(SolverContext &solver_context) + : solver_context(solver_context) + , fe_values_flags() + , num_quad_points(SystemWide::quadrature.size()) + , max_refinement(solver_context.get_user_inputs() + .get_spatial_discretization() + .get_max_refinement()) + , min_refinement(solver_context.get_user_inputs() + .get_spatial_discretization() + .get_min_refinement()) + , marker_functions() + { + fe_values_flags.fill(dealii::UpdateFlags::update_default); + for (const auto &criterion : solver_context.get_user_inputs() + .get_spatial_discretization() + .get_refinement_criteria()) + { + // Grab the index and field type + const Types::Index index = criterion.get_index(); + const FieldInfo::TensorRank rank = solver_context.get_user_inputs() + .get_variable_attributes() + .at(index) + .field_info.tensor_rank; + + if (criterion.get_criterion() & GridRefinement::RefinementFlags::Value) + { + fe_values_flags[int(rank)] |= dealii::UpdateFlags::update_values; + } + else if (criterion.get_criterion() & GridRefinement::RefinementFlags::Gradient) + { + fe_values_flags[int(rank)] |= dealii::UpdateFlags::update_gradients; + } + } + // Create the FEValues + for (const auto field_type : + {FieldInfo::TensorRank::Scalar, FieldInfo::TensorRank::Vector}) + { + fe_values[field_type] = + dealii::FEValues(SystemWide::fe_systems[field_type], + SystemWide::quadrature, + fe_values_flags.at(int(field_type))); + } + } + + /** + * @brief Destructor. + */ + ~RefinementManager() = default; + + /** + * @brief Copy constructor. + * + * Deleted so grid refiner instances aren't copied. + */ + RefinementManager(const RefinementManager &grid_refiner) = delete; + + /** + * @brief Copy assignment. + * + * Deleted so grid refiner instances aren't copied. + */ + RefinementManager & + operator=(const RefinementManager &grid_refiner) = delete; + + /** + * @brief Move constructor. + * + * Deleted so grid refiner instances aren't moved. + */ + RefinementManager(RefinementManager &&grid_refiner) noexcept = delete; + + /** + * @brief Move assignment. + * + * Deleted so grid refiner instances aren't moved. + */ + RefinementManager & + operator=(RefinementManager &&grid_refiner) noexcept = delete; + + /** + * @brief Do the adaptive refinement + * + * This function consists of a few steps. + * 1. Flag the cells on the mesh according to the refinement criterion. + * 2. Refine the grid and initialize the solution transfer object + * 3. Redistribute the DoFs + * 4. Transfer the solution from old to new + * 5. Recompute and reapply the constraints (this is done in the solvers) + * 6. Recompute invm & element volume (if applicable) + */ + void + do_adaptive_refinement() + { + // Return early if adaptive meshing is disabled + if (!solver_context.get_user_inputs() + .get_spatial_discretization() + .get_has_adaptivity()) + { + return; + } + + TriangulationManager &triangulation_manager = + solver_context.get_triangulation_manager(); + DofManager &dof_manager = solver_context.get_dof_manager(); + ConstraintManager &constraint_manager = + solver_context.get_constraint_manager(); + std::vector &field_attributes = + solver_context.get_field_attributes(); + + // Step 1 + mark_cells_for_refinement_and_coarsening(); + bool first_iteration = true; + while ( + dealii::Utilities::MPI::logical_or(mark_cells_for_refinement(), MPI_COMM_WORLD) || + first_iteration) + { + first_iteration = false; + + // Update ghosts of all fields. + for (int field_index = 0; field_index < field_attributes.size(); field_index++) + { + solver_context.get_solution_indexer() + .get_solution(field_index) + .update_ghost_values(); + } + + // Step 2 + refine_grid(); + + // Step 3 + triangulation_manager.reinit(); + dof_manager.reinit(triangulation_manager, field_attributes); + constraint_manager().make_constraints(SystemWide::mapping, + dof_manager.get_dof_handlers()); + + // Todo: reinit matrix free operators? + + // Todo: reinit solutions? + + // Step 4 + // Todo execute solution transfer + + // Step 6 + // Todo: Recompute invm & element volume + } + } + + void + add_refinement_marker(std::shared_ptr> marker) + { + marker_functions.push_back(marker); + } + + void + clear_refinement_markers() + { + marker_functions.clear(); + } + + const std::vector>> & + get_refinement_markers() const + { + return marker_functions; + } + +private: + /** + * @brief Mark cells for refinement and coarsening + */ + void + mark_cells_for_refinement_and_coarsening() + { + // Create the an object for the refinement criterion at each of the quad points. This + // will either contain the value for scalar fields, the magnitude for vector fields, + // or the magnitude of the gradient for both of the fields. + std::vector values(num_quad_points, 0.0); + + // Clear user flags + solver_context.get_triangulation_handler().clear_user_flags(); + + // Loop over the cells provided by the triangulation + for (const auto &cell : solver_context.get_triangulation_handler() + .get_triangulation() + .active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + // Whether we should refine the cell + bool should_refine = false; + + // TODO (landinjm): We can probably avoid checking some of the neighboring + // cells when coarsening them + for (const auto &criterion : solver_context.get_user_inputs() + .get_spatial_discretization() + .get_refinement_criteria()) + { + // Grab the index + const Types::Index index = criterion.get_index(); + + // Grab the field type + const FieldInfo::TensorRank local_field_type = + solver_context.get_user_inputs() + .get_variable_attributes() + .at(index) + .field_info.tensor_rank; + + // Grab the DoFHandler iterator + const auto dof_iterator = cell->as_dof_handler_iterator( + solver_context.get_dof_handler().get_dof_handler(index)); + + // Reinit the cell + fe_values.at(local_field_type).reinit(dof_iterator); + + if (criterion.get_criterion() & GridRefinement::RefinementFlags::Value) + { + if (local_field_type == FieldInfo::TensorRank::Scalar) + { + // Get the values for a scalar field + fe_values.at(local_field_type) + .get_function_values( + *solver_context.get_solution_handler() + .get_solution_vector(index, DependencyType::Normal), + values); + } + else + { + // Get the magnitude of the value for vector fields + // TODO (landinjm): Should be zeroing this out? + std::vector> vector_values( + num_quad_points, + dealii::Vector(dim)); + fe_values.at(local_field_type) + .get_function_values( + *solver_context.get_solution_handler() + .get_solution_vector(index, DependencyType::Normal), + vector_values); + for (unsigned int q_point = 0; q_point < num_quad_points; + ++q_point) + { + values[q_point] = vector_values[q_point].l2_norm(); + } + } + + // Check if any of the quadrature points meet the refinement criterion + for (unsigned int q_point = 0; q_point < num_quad_points; ++q_point) + { + if (criterion.value_in_open_range(values[q_point])) + { + should_refine = true; + break; + } + } + + // Exit if we've already determined that the cell has to be refined + if (should_refine) + { + break; + } + } + if (criterion.get_criterion() & GridRefinement::RefinementFlags::Gradient) + { + if (local_field_type == FieldInfo::TensorRank::Scalar) + { + // Get the magnitude of the gradient for a scalar field + // TODO (landinjm): Should be zeroing this out? + std::vector> scalar_gradients( + num_quad_points); + fe_values.at(local_field_type) + .get_function_gradients( + *solver_context.get_solution_handler() + .get_solution_vector(index, DependencyType::Normal), + scalar_gradients); + for (unsigned int q_point = 0; q_point < num_quad_points; + ++q_point) + { + values[q_point] = scalar_gradients[q_point].norm(); + } + } + else + { + // TODO (landinjm): Should be zeroing this out? + std::vector>> + vector_gradients(num_quad_points, + std::vector>( + dim)); + fe_values.at(local_field_type) + .get_function_gradients( + *solver_context.get_solution_handler() + .get_solution_vector(index, DependencyType::Normal), + vector_gradients); + for (unsigned int q_point = 0; q_point < num_quad_points; + ++q_point) + { + dealii::Vector vector_gradient_component_magnitude( + dim); + for (unsigned int dimension = 0; dimension < dim; dimension++) + { + vector_gradient_component_magnitude[dimension] = + vector_gradients[q_point][dimension].norm(); + } + values[q_point] = + vector_gradient_component_magnitude.l2_norm(); + } + } + + // Check if any of the quadrature points meet the refinement criterion + for (unsigned int q_point = 0; q_point < num_quad_points; ++q_point) + { + if (criterion.gradient_magnitude_above_threshold(values[q_point])) + { + should_refine = true; + break; + } + } + + // Exit if we've already determined that the cell has to be refined + if (should_refine) + { + break; + } + } + } + + Assert(cell->level() > 0, + dealii::ExcMessage("Cell refinement level is less than one, which " + "will lead to underflow.")); + const auto cell_refinement = static_cast(cell->level()); + if (should_refine && cell_refinement < max_refinement) + { + cell->set_user_flag(); + cell->clear_coarsen_flag(); + cell->set_refine_flag(); + } + if (should_refine) + { + cell->set_user_flag(); + cell->clear_coarsen_flag(); + } + if (!should_refine && cell_refinement > min_refinement && + !cell->user_flag_set()) + { + cell->set_coarsen_flag(); + } + } + } + } + + /** + * @brief Mark cells based on function. Note: cells are only marked for refinement but + * not coarsening. + * @param refinement_function A function that determines if a cell should be refined. + * @return True if any cell was marked for refinement, false otherwise. + */ + bool + mark_cells_for_refinement() + { + bool any_cell_marked = false; + for (const auto &cell : solver_context.get_triangulation_manager() + .get_triangulation() + .active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + const auto cell_refinement = static_cast(cell->level()); + if (std::any_of( + marker_functions.begin(), + marker_functions.end(), + [&](const std::shared_ptr> &marker_function) + { + return marker_function->flag( + *cell, + solver_context.get_user_inputs().get_temporal_discretization()); + })) + { + cell->set_user_flag(); + cell->clear_coarsen_flag(); + if (cell_refinement < max_refinement) + { + cell->set_refine_flag(); + any_cell_marked = true; + } + } + } + } + return any_cell_marked; + } + + /** + * @brief Refine the grid + */ + void + refine_grid() + { + // Prepare for grid refinement + solver_context.get_triangulation_manager().prepare_for_grid_refinement(); + + // Prepare the solution transfer objects + solver_context.get_solution_handler().prepare_for_solution_transfer(); + + // Execute grid refinement + solver_context.get_triangulation_manager().execute_grid_refinement(); + } + + /** + * @brief Grid refinement context. + */ + SolverContext solver_context; + + /** + * @brief Update flags for the FEValues determined by the grid refinement + * criterion. For now, we share one flag set for scalar fields and one for vector + * fields. + */ + std::array fe_values_flags; + + /** + * @brief Finite element values for scalar and vector fields. + */ + std::array, 2> fe_values; + + /** + * @brief Number of quadrature points. + */ + unsigned int num_quad_points = 0; + + /** + * @brief Maximum global refinement level. + */ + unsigned int max_refinement = 0; + + /** + * @brief Minimum global refinement level. + */ + unsigned int min_refinement = 0; + + /** + * @brief Marker functions. + */ + std::list>> marker_functions; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/simulation_timer.h b/include/prismspf/core/simulation_timer.h new file mode 100644 index 000000000..e155660e2 --- /dev/null +++ b/include/prismspf/core/simulation_timer.h @@ -0,0 +1,83 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +class SimulationTimer +{ +public: + SimulationTimer() = default; + + explicit SimulationTimer(double step_size) + : time_step_size(step_size) + {} + + [[nodiscard]] unsigned int + get_increment() const + { + return current_increment; + } + + [[nodiscard]] double + get_time() const + { + return current_time; + } + + [[nodiscard]] double + get_timestep() const + { + return time_step_size; + } + + void + increment(double step_size) + { + current_increment++; + current_time += step_size; + } + + void + increment() + { + increment(time_step_size); + } + + void + set_timestep(double step_size) + { + time_step_size = step_size; + } + + void + set_time(double time) + { + current_time = time; + } + + void + set_increment(unsigned int increment) + { + current_increment = increment; + } + + void + reset() + { + current_increment = 0; + current_time = 0.0; + } + +private: + unsigned int current_increment = 0; + double current_time = 0.0; + double time_step_size = 0.0; +}; + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/include/prismspf/core/solution_indexer.h b/include/prismspf/core/solution_indexer.h new file mode 100644 index 000000000..18b947007 --- /dev/null +++ b/include/prismspf/core/solution_indexer.h @@ -0,0 +1,107 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief Class that provides access to solution vectors spread across different + * groups. + */ +template +class SolutionIndexer +{ +public: + using BlockVector = GroupSolutionHandler::BlockVector; + using SolutionVector = BlockVector::BlockType; + using MatrixFree = GroupSolutionHandler::MatrixFree; + + /** + * @brief Constructor. + */ + SolutionIndexer(unsigned int num_fields, // attributes_list.size() + const std::set> &solution_handlers); + + /** + * @brief Get a solution vector of a given field index. + */ + [[nodiscard]] const SolutionVector & + get_solution_vector(unsigned int global_index, unsigned int relative_level = 0) const; + /** + * @brief Get a solution vector of a given field index. + */ + [[nodiscard]] SolutionVector & + get_solution_vector(unsigned int global_index, unsigned int relative_level = 0); + + /** + * @brief Get a solution vector of a given field index at a given age. + */ + [[nodiscard]] const SolutionVector & + get_old_solution_vector(unsigned int age, + unsigned int global_index, + unsigned int relative_level = 0) const; + /** + * @brief Get a solution vector of a given field index at a given age. + */ + [[nodiscard]] SolutionVector & + get_old_solution_vector(unsigned int age, + unsigned int global_index, + unsigned int relative_level = 0); + + /** + * @brief Get the "new" solution vector set. + */ + [[nodiscard]] const BlockVector & + get_new_solution_full_vector(unsigned int relative_level = 0) const; + /** + * @brief Get the "new" solution vector set. + */ + [[nodiscard]] BlockVector & + get_new_solution_full_vector(unsigned int relative_level = 0); + + /** + * @brief Get the "new" solution vector of a given field index. + */ + [[nodiscard]] const SolutionVector & + get_new_solution_vector(unsigned int index, unsigned int relative_level = 0) const; + /** + * @brief Get the "new" solution vector of a given field index. + */ + [[nodiscard]] SolutionVector & + get_new_solution_vector(unsigned int index, unsigned int relative_level = 0); + + /** + * @brief Get the matrixfree object of the group a given field index. + */ + [[nodiscard]] const MatrixFree & + get_matrix_free(unsigned int index, unsigned int relative_level = 0) const; + /** + * @brief Get the matrixfree object of the group a given field index. + */ + [[nodiscard]] MatrixFree & + get_matrix_free(unsigned int index, unsigned int relative_level = 0); + + /** + * @brief Get the matrixfree object of the group a given field index. + */ + [[nodiscard]] std::pair &, unsigned int> + get_solution_level_and_block_index(unsigned int index, unsigned int relative_level = 0); + + /** + * @brief Get the block index of a field within its solve group + */ + [[nodiscard]] unsigned int + get_block_index(unsigned int global_index) const; + +private: + std::vector *> solutions; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/solution_output.h b/include/prismspf/core/solution_output.h index 6630d8684..03af39611 100644 --- a/include/prismspf/core/solution_output.h +++ b/include/prismspf/core/solution_output.h @@ -5,6 +5,9 @@ #include +#include +#include + #include #include @@ -40,6 +43,16 @@ class SolutionOutput const unsigned int °ree, const std::string &name, const UserInputParameters &user_inputs); + + /** + * @brief Outputs all fields in the solution set. + */ + SolutionOutput(const std::vector &field_attributes, + const SolutionIndexer &solution_indexer, + const DofManager &dof_manager, + const unsigned int °ree, + const std::string &file_prefix, + const UserInputParameters &user_inputs); }; PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/solve_group.h b/include/prismspf/core/solve_group.h new file mode 100644 index 000000000..5e9db0f6a --- /dev/null +++ b/include/prismspf/core/solve_group.h @@ -0,0 +1,134 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +// NOLINTBEGIN(misc-non-private-member-variables-in-classes, hicpp-explicit-conversions) +// readability-simplify-boolean-expr + +/** + * @brief Structure to hold the attributes of a solve-group. + */ +class SolveGroup +{ +public: + using EvalFlags = dealii::EvaluationFlags::EvaluationFlags; + using FieldType = FieldInfo::TensorRank; + + explicit SolveGroup(int _id = -1, + PDEType _pde_type = PDEType::Explicit, + bool _is_time_dependent = true, + std::set _field_indices = {}, + DependencySet _dependencies_rhs = {}, + DependencySet _dependencies_lhs = {}) + : id(_id) + , pde_type(_pde_type) + , is_time_dependent(_is_time_dependent) + , field_indices(std::move(_field_indices)) + , dependencies_rhs(std::move(_dependencies_rhs)) + , dependencies_lhs(std::move(_dependencies_lhs)) + {} + + /** + * @brief Unique identifier. Use this in 'if' statements or switch cases in equations + * lhs and rhs. + */ + int id; + + /** + * @brief PDE type (Constant | Explicit | Linear | Newton). + */ + PDEType pde_type; + + /** + * @brief Whether the solve is time dependent. This is used to determine whether to + * initialize the solution vector with the initial conditions or just solve. + */ + bool is_time_dependent; + + /** + * @brief Indices of the fields to be solved in this group. + */ + std::set field_indices; + + /** + * @brief Dependencies for the rhs equation(s) + */ + DependencySet dependencies_rhs; + /** + * @brief Dependencies for the lhs equation(s) + */ + DependencySet dependencies_lhs; + +private: + /** + * @brief A std::vector used to manage dynamic memory for holding an auxiliary solve + * group + * @remark Privately managed + */ + std::vector aux_solve_container; + +public: + void + set_auxiliary_solve(const SolveGroup &aux_solve = SolveGroup()) + { + aux_solve_container.clear(); + aux_solve_container.push_back(aux_solve); + } + + void + remove_auxiliary_solve() // unset_auxiliary_solve() + { + aux_solve_container.clear(); + } + + [[nodiscard]] bool + has_auxiliary_solve() const + { + return !aux_solve_container.empty(); + } + + SolveGroup & + auxiliary_solve() // get_auxiliary_solve() + { + AssertThrow(has_auxiliary_solve(), + dealii::ExcMessage("Cannot access auxiliary solve when none is set.\n " + "Attempted access from solve id: " + + std::to_string(id))); + return aux_solve_container[0]; + } + + bool + operator<(const SolveGroup &other) const + { + // TODO: Ensure that the pde_type enum is defined in the proper order. + // Swap Explicit and ImplicitTimeDependent? + // Constant | ExplicitTimeDependent | Explicit | ImplicitTimeDependent | Implicit | + // Postprocess + if (pde_type < other.pde_type) + { + return true; + } + return id < other.id; + } +}; + +// NOLINTEND(misc-non-private-member-variables-in-classes, hicpp-explicit-conversions) + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/system_wide.h b/include/prismspf/core/system_wide.h new file mode 100644 index 000000000..4652c682f --- /dev/null +++ b/include/prismspf/core/system_wide.h @@ -0,0 +1,45 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief This is the main class that handles the construction and solving of + * user-specified PDEs. + */ +template +class SystemWide +{ +public: + /** + * @brief Scalar and Vector FE systems. + */ + inline static const std::array, 2> fe_systems = { + dealii::FESystem(dealii::FE_Q(dealii::QGaussLobatto<1>(degree + 1)), 1), + dealii::FESystem(dealii::FE_Q(dealii::QGaussLobatto<1>(degree + 1)), dim)}; + + /** + * @brief Mappings to and from reference cell. + */ + inline static const dealii::MappingQ1 mapping; + + /** + * @brief Quadrature rule + */ + inline static const dealii::QGaussLobatto quadrature = + dealii::QGaussLobatto(degree + 1); +}; + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/include/prismspf/core/triangulation_manager.h b/include/prismspf/core/triangulation_manager.h new file mode 100644 index 000000000..f6f6baab0 --- /dev/null +++ b/include/prismspf/core/triangulation_manager.h @@ -0,0 +1,168 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include + +#include +#include +#include + +#include + +#include + + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief This class handlers the generation and manipulation of triangulations. + */ +template +class TriangulationManager +{ +public: + using Triangulation = + std::conditional_t, + dealii::parallel::distributed::Triangulation>; + + /** + * @brief Constructor. + */ + explicit TriangulationManager(bool _has_multigrid); + + /** + * @brief Set whether multigrid triangulations will be generated. + */ + void + set_using_multigrid(bool _has_multigrid) + { + has_multigrid = _has_multigrid; + } + + /** + * @brief Reinitialize the triangulation handler. + * This is used for AMR with multigrid so the coarsened meshes can be reinitialized. + */ + void + reinit() + { + // Check that the initial global refinement matches the maximum adaptive refinement + /* Assert(user_inputs->get_spatial_discretization().get_global_refinement() == + user_inputs->get_spatial_discretization().get_max_refinement(), + dealii::ExcMessage( + "Currently, we don't allow the initial refinement to be lower than the " + "maximum adpative refinement level when using multigrid. This is because we " + "have to create a sequence of coarser meshes.")); */ + if (has_multigrid) + { + coarsened_triangulations = + dealii::MGTransferGlobalCoarseningTools::create_geometric_coarsening_sequence( + triangulation); + std::reverse(coarsened_triangulations.begin(), coarsened_triangulations.end()); + } + else + { + coarsened_triangulations.clear(); + coarsened_triangulations.push_back(&triangulation); + } + // TODO (landinjm): p-multigrid + }; + + /** + * @brief Getter function for triangulation (constant reference). + */ + [[nodiscard]] const Triangulation & + get_triangulation() const; + + /** + * @brief Getter function for the multigrid triangulation (constant reference). + */ + [[nodiscard]] const std::vector>> & + get_mg_triangulation() const; + + /** + * @brief Getter function for a level in the globally coarsening multigrid triangulation + * (constant reference). + */ + [[nodiscard]] const dealii::Triangulation & + get_triangulation(unsigned int relative_level) const; + + /** + * @brief Generate mesh based on the inputs provided by the user. + */ + void + generate_mesh(const UserInputParameters &user_inputs); + + /** + * @brief Export triangulation to vtk. This is done for debugging purposes when dealing + * with unusual meshes (e.g., circular domains). + */ + void + export_triangulation_as_vtk(const std::string &filename) const; + + /** + * @brief Prepare the triangulation for grid refinement. + */ + void + prepare_for_grid_refinement() + { + triangulation.prepare_coarsening_and_refinement(); + } + + /** + * @brief Execute grid refinement on the triangulation. + */ + void + execute_grid_refinement() + { + triangulation.execute_coarsening_and_refinement(); + } + + /** + * @brief Clear all user flags. + */ + void + clear_user_flags() + { + triangulation.clear_user_flags(); + } + +private: + /** + * @brief Mark the domain ids on the triangulation to get the proper mapping of + * specified boundary conditions. + * + * TODO (landinjm): When the user has different meshs (or custom for that matter), we + * should let them manually set the boundary ids. + */ + void + mark_boundaries(const SpatialDiscretization &discretization_params) const; + + /** + * @brief Mark certain faces of the triangulation periodic. + */ + void + mark_periodic(const UserInputParameters &user_inputs); + + /** + * @brief Whether we have multigrid. + */ + bool has_multigrid = false; + + /** + * @brief Main triangulation. + */ + dealii::Triangulation triangulation; + + /** + * @brief Coarsened triangulations for multigrid. + */ + std::vector>> coarsened_triangulations; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/core/type_enums.h b/include/prismspf/core/type_enums.h index eb5f52979..d87f9728b 100644 --- a/include/prismspf/core/type_enums.h +++ b/include/prismspf/core/type_enums.h @@ -19,7 +19,10 @@ enum PDEType : std::uint8_t ImplicitTimeDependent, TimeIndependent, Auxiliary, - Constant + Constant, + Explicit, + Linear, + Newton }; /** @@ -69,9 +72,9 @@ enum FieldSolveType : std::uint8_t /** * @brief Internal classification for types of variable dependencies. */ -enum DependencyType : std::uint8_t +enum DependencyType : int { - Normal, + Normal = -1, Change, OldOne, OldTwo, diff --git a/include/prismspf/core/types.h b/include/prismspf/core/types.h index 7c7b283c3..d17b45992 100644 --- a/include/prismspf/core/types.h +++ b/include/prismspf/core/types.h @@ -30,12 +30,17 @@ namespace Numbers /** * @brief Max element degree. */ - static const unsigned int max_element_degree = 6; + static constexpr unsigned int max_element_degree = 6; + + /** + * @brief How many increments of field history prismspf can store. + */ + static constexpr unsigned int max_saved_increments = 4; /** * @brief Max number of subsections. */ - static const unsigned int max_subsections = 16; + static constexpr unsigned int max_subsections = 16; /** * @brief Invalid PDE type. diff --git a/include/prismspf/core/variable_attributes.h b/include/prismspf/core/variable_attributes.h index 28acc2dfb..8cd087718 100644 --- a/include/prismspf/core/variable_attributes.h +++ b/include/prismspf/core/variable_attributes.h @@ -83,9 +83,9 @@ struct FieldInfo */ enum class TensorRank : std::uint8_t { - Undefined, - Scalar, - Vector + Undefined = static_cast(-1), + Scalar = 0, + Vector = 1 }; /** diff --git a/include/prismspf/solvers/group_explicit_solver.h b/include/prismspf/solvers/group_explicit_solver.h new file mode 100644 index 000000000..4a6bb7a27 --- /dev/null +++ b/include/prismspf/solvers/group_explicit_solver.h @@ -0,0 +1,90 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include + +#include +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class SolverContext; + +/** + * @brief This class handles the explicit solves of all explicit fields + */ +template +class ExplicitSolver : public GroupSolverBase +{ +public: + using GroupSolverBase = GroupSolverBase; + using GroupSolverBase::rhs_operators; + using GroupSolverBase::solutions; + using GroupSolverBase::solve_group; + using GroupSolverBase::solver_context; + + /** + * @brief Constructor. + */ + ExplicitSolver(SolveGroup _solve_group, + const SolverContext &_solver_context) + : GroupSolverBase(_solve_group, _solver_context) + {} + + /** + * @brief Solve for a single update step. + */ + void + solve_level(unsigned int relative_level) override + { + // Zero out the ghosts + Timer::start_section("Zero ghosts"); + solutions.zero_out_ghosts(relative_level); + Timer::end_section("Zero ghosts"); + + rhs_operators[relative_level].compute_operator( + solutions.get_new_solution_full_vector(relative_level)); + + // TODO: if it's used for all solvers, define invm in solution handler. originally + // this was done as a loop over fields. Scale the update by the respective + // (Scalar/Vector) invm. + solutions.get_new_solution_full_vector(relative_level).scale(block_invm); + + // Update the solutions + solutions.update(relative_level); + + // Apply constraints + solutions.apply_constraints(relative_level); + + // Update the ghosts + Timer::start_section("Update ghosts"); + solutions.update_ghosts(relative_level); + Timer::end_section("Update ghosts"); + } + + /** + * @brief Solve for a single update step. + */ + void + solve() override; + + /** + * @brief Print information about the solver to summary.log. + */ + void + print() override; + +private: + /** + * @brief Matrix free operators for each level + */ + // std::vector> rhs_operators; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/solvers/group_linear_solver.h b/include/prismspf/solvers/group_linear_solver.h new file mode 100644 index 000000000..2968c988d --- /dev/null +++ b/include/prismspf/solvers/group_linear_solver.h @@ -0,0 +1,132 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include + +#include +#include + +#include +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class SolverContext; + +/** + * @brief This class handles the explicit solves of all explicit fields + */ +template +class LinearSolver : public GroupSolverBase +{ + using GroupSolverBase = GroupSolverBase; + using GroupSolverBase::rhs_operators; + using GroupSolverBase::solutions; + using GroupSolverBase::solve_group; + using GroupSolverBase::solver_context; + using BlockVector = SolutionHandler::BlockVector; + +public: + /** + * @brief Constructor. + */ + LinearSolver(SolveGroup _solve_group, + const SolverContext &_solver_context) + : GroupSolverBase(_solve_group, _solver_context) + {} + + /** + * @brief Initialize the solver. + */ + void + init() override + { + GroupSolverBase::init(); + // Initialize lhs_operators + for (unsigned int relative_level = 0; relative_level < lhs_operators.size(); + ++relative_level) + { + lhs_operators[relative_level] = MFOperator( + solver_context->pde_operator, + PDEOperator::compute_explicit_lhs, + solver_context->field_attributes, + solver_context->solution_indexer, + relative_level, + solve_group.dependencies_lhs); + lhs_operators[relative_level].initialize(solve_group, + solutions.get_matrix_free( + relative_level), + solutions.get_global_to_block_index()); + } + } + + /** + * @brief Solve for a single update step. + */ + void + solve_level(unsigned int relative_level) override + { + // Update the solutions + solutions.update(relative_level); + + // Zero out the ghosts + Timer::start_section("Zero ghosts"); + solutions.zero_out_ghosts(relative_level); + Timer::end_section("Zero ghosts"); + + // Set up linear solver + do_linear_solve(rhs_operators[relative_level], + solutions.get_new_solution_full_vector(relative_level), + lhs_operators[relative_level], + solutions.get_solution_full_vector(relative_level)); + + // Apply constraints + solutions.apply_constraints(relative_level); + + // Update the ghosts + Timer::start_section("Update ghosts"); + solutions.update_ghosts(relative_level); + Timer::end_section("Update ghosts"); + } + + void + do_linear_solve(MFOperator &rhs_operator, + BlockVector &b_vector, + MFOperator &lhs_operator, + BlockVector &x_vector) + { + // Compute rhs + rhs_operator.compute_operator(b_vector); + // Linear solve + try + { + // TODO: Make member? + dealii::SolverCG cg_solver(linear_solver_control); + cg_solver.solve(lhs_operator, x_vector, b_vector); + } + catch (...) // TODO: more specific catch + { + ConditionalOStreams::pout_base() + << "Warning: linear solver did not converge as per set tolerances.\n"; + } + } + +protected: + /** + * @brief Matrix free operators for each level + */ + std::vector> lhs_operators; + +private: + /** + * @brief Solver control. Contains max iterations and tolerance. + */ + dealii::SolverControl linear_solver_control; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/solvers/group_mg_solver.h b/include/prismspf/solvers/group_mg_solver.h new file mode 100644 index 000000000..7b8fff8c3 --- /dev/null +++ b/include/prismspf/solvers/group_mg_solver.h @@ -0,0 +1,159 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class SolverContext; + +/** + * @brief This class handles the explicit solves of all explicit fields + */ +template +class MGSolver : public LinearSolver +{ + using GroupSolverBase = GroupSolverBase; + using LinearSolver = LinearSolver; + using GroupSolverBase::rhs_operators; + using GroupSolverBase::solutions; + using GroupSolverBase::solve_group; + using GroupSolverBase::solver_context; + using LinearSolver::do_linear_solve; + using LinearSolver::lhs_operators; + using BlockVector = SolutionHandler::BlockVector; + +public: + /** + * @brief Constructor. + */ + MGSolver(SolveGroup _solve_group, + const SolverContext &_solver_context) + : LinearSolver(_solve_group, _solver_context) + {} + + /** + * @brief Solve for a single update step. + */ + void + solve_level(unsigned int relative_level) override + {} + + void + mg_solve() + { + int min_level = 0; + int max_level = 0; + // 1. Level operators + dealii::MGLevelObject> mg_lhs_operators(min_level, + max_level); + // TODO: Initialize them + + // 2. Solution and rhs vectors (rhs and solution here are used as scratch space by mg, + // but are not actually the same as the top level vectors) + dealii::MGLevelObject mg_solution_vectors(min_level, max_level); + dealii::MGLevelObject mg_rhs_vectors(min_level, max_level); + + // TODO: Initialize their shapes + + // 3. MG transfer + dealii::MGTransferBlockMatrixFree mg_transfer; // Constraints? + // NOTE: dof_handler.distribute_mg_dofs() must have been called + mg_transfer.build( + solver_context->dof_manager.get_dof_handlers(solve_group.field_indices, 0)); + + // 4. MG Smoother (takes in operators) This is similar to a solver, but is + // conceptually different. + // Preconditioner for smoother. + // TODO: use PreconditionBlockJacobi or other block preconditioner instead + using SmootherPrecond = + dealii::PreconditionChebyshev, BlockVector>; + dealii::MGLevelObject smoother_data( + min_level, + max_level); + for (unsigned int level = min_level; level <= max_level; ++level) + { + // smoother_data[level].smoothing_range; + // smoother_data[level].degree; + // smoother_data[level].eig_cg_n_iterations; + // smoother_data[level].preconditioner; + // smoother_data[level].constraints; + } + // Wrapper around a generic preconditioner to be used as a smoother + using Smoother = dealii::MGSmootherPrecondition, + SmootherPrecond, + BlockVector>; + Smoother mg_smoother; + mg_smoother.initialize(mg_lhs_operators, smoother_data); + + // The MG method requires the choice of a solve on the coarsest level. There are a + // couple of main options: + // A. Direct solve // we're not doing this + // B. Iterative solve (using a Krylov method, (like CG)) + /* dealii::MGCoarseGridIterativeSolver, + MFOperator, + dealii::PreconditionIdentity> mg_coarse_solver;*/ + // C. Smoothing operation (don't actually solve, just use the same smoother as on + // other levels) + /* dealii::MGCoarseGridApplySmoother mg_coarse_solver; + mg_coarse_solver.initialize(mg_smoother); */ + + // 5. Coarse grid solver + dealii::MGCoarseGridApplySmoother mg_coarse_solver; + mg_coarse_solver.initialize(mg_smoother); + + // 6. Multigrid object + dealii::Multigrid multigrid( + mg_lhs_operators, + mg_coarse_solver, + mg_transfer, + mg_smoother, + mg_smoother, + min_level, + max_level, + dealii::Multigrid::Cycle::v_cycle); + + // 7. Turn MG into a preconditioner object + dealii:: + PreconditionMG> + preconditioner_mg( + solver_context->dof_manager.get_dof_handlers(solve_group.field_indices, 0), + multigrid, + mg_transfer); + + // 8. Solve the system with CG + MG preconditioner + dealii::SolverControl cg_solver_control(/* TODO */); + dealii::SolverCG cg_solver(cg_solver_control); + } + +private: + /** + * @brief Solver control. Contains max iterations and tolerance. + */ + dealii::SolverControl linear_solver_control; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/solvers/group_newton_solver.h b/include/prismspf/solvers/group_newton_solver.h new file mode 100644 index 000000000..16c5a3646 --- /dev/null +++ b/include/prismspf/solvers/group_newton_solver.h @@ -0,0 +1,105 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include + +#include + +#include +#include +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class SolverContext; + +/** + * @brief This class handles the explicit solves of all explicit fields + */ +template +class NewtonSolver : public LinearSolver +{ + using GroupSolverBase = GroupSolverBase; + using LinearSolver = LinearSolver; + using GroupSolverBase::rhs_operators; + using GroupSolverBase::solutions; + using GroupSolverBase::solve_group; + using GroupSolverBase::solver_context; + using LinearSolver::do_linear_solve; + using LinearSolver::lhs_operators; + +public: + /** + * @brief Constructor. + */ + NewtonSolver(SolveGroup _solve_group, + const SolverContext &_solver_context) + : LinearSolver(_solve_group, _solver_context) + {} + + /** + * @brief Solve for a single update step. + */ + void + solve_level(unsigned int relative_level) override + { + number newton_step_length = number(1.0); // TODO + number newton_tolerance = number(1.0); // TODO + unsigned int newton_max_iterations = 1; // TODO + bool newton_unconverged = true; + unsigned int iter = 0; + + // Update the solutions + solutions.update(relative_level); + + // TODO: setup initial guess for solution vector. Maybe use old_solution if available. + + // Newton iteration loop. + while (newton_unconverged && iter++ < newton_max_iterations) + { + // Zero out the ghosts + Timer::start_section("Zero ghosts"); + solutions.zero_out_ghosts(relative_level); + Timer::end_section("Zero ghosts"); + + // Solve for Newton update. + do_linear_solve(rhs_operators[relative_level], + solutions.get_new_solution_full_vector(relative_level), + lhs_operators[relative_level], + solutions.get_change_solution_full_vector(relative_level)); + + // TODO: Apply zero-constraints to 'change' vector + + // Perform Newton update. + solutions.get_solution_full_vector(relative_level) + .add(newton_step_length, + solutions.get_change_solution_full_vector(relative_level)); + + // Apply constraints to solution vector + solutions.apply_constraints(relative_level); + + // Update the ghosts + Timer::start_section("Update ghosts"); + solutions.update_ghosts(relative_level); + Timer::end_section("Update ghosts"); + + // Check convergence. + number l2_norm = + solutions.get_change_solution_full_vector(relative_level).l2_norm(); + newton_unconverged = l2_norm > newton_tolerance; + } + ConditionalOStreams::pout_verbose() << iter << " Newton iterations to converge.\n\n"; + if (iter >= newton_max_iterations) + { + ConditionalOStreams::pout_base() << "Warning: nonlinear solver did not " + "converge as per set tolerances.\n\n"; + } + } +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/solvers/group_solver_base.h b/include/prismspf/solvers/group_solver_base.h new file mode 100644 index 000000000..eb53122d2 --- /dev/null +++ b/include/prismspf/solvers/group_solver_base.h @@ -0,0 +1,287 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include + +#include + +#include // +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class GroupSolverBase +{ +public: + /** + * @brief Constructor. + */ + GroupSolverBase(SolveGroup _solve_group, + const SolverContext &_solver_context) + : solve_group(std::move(_solve_group)) + , solutions() + , solver_context( + std::make_shared>(_solver_context)) + {} + + /** + * @brief Destructor. + */ + virtual ~GroupSolverBase() = default; + + /** + * @brief Copy constructor. + * + * Deleted so solver instances aren't copied. + */ + GroupSolverBase(const GroupSolverBase &solver_base) = delete; + + /** + * @brief Copy assignment. + * + * Deleted so solver instances aren't copied. + */ + GroupSolverBase & + operator=(const GroupSolverBase &solver_base) = delete; + + /** + * @brief Move constructor. + * + * Deleted so solver instances aren't moved. + */ + GroupSolverBase(GroupSolverBase &&solver_base) noexcept = delete; + + /** + * @brief Move assignment. + * + * Deleted so solver instances aren't moved. + */ + GroupSolverBase & + operator=(GroupSolverBase &&solver_base) noexcept = delete; + + /** + * @brief Initialize the solver. + */ + virtual void + init() + { + // Initialize vectors + solutions = + GroupSolutionHandler(solve_group, solver_context->field_attributes); + solutions.init(solver_context->mapping, + solver_context->dof_manager, + solver_context->constraint_manager, + solver_context->quadrature); + + // Set the initial condition + set_initial_condition(); + + // Apply constraints. + solutions.apply_constraints(); + + // Initialize rhs_operators + for (unsigned int relative_level = 0; relative_level < rhs_operators.size(); + ++relative_level) + { + rhs_operators[relative_level] = MFOperator( + solver_context->pde_operator, + PDEOperator::compute_explicit_rhs, + solver_context->field_attributes, + solver_context->solution_indexer, + relative_level, + solve_group.dependencies_rhs); + rhs_operators[relative_level].initialize(solve_group, + solutions.get_matrix_free( + relative_level), + solutions.get_global_to_block_index()); + } + } + + /** + * @brief Reinitialize the solution vectors & apply constraints. + */ + virtual void + reinit() + { + // Apply constraints. + solutions.reinit(); + solutions.apply_constraints(); + } + + /** + * @brief Solve one level. + */ + virtual void + solve_level([[maybe_unused]] unsigned int relative_level = 0) + {} + + /** + * @brief Solve for a single update step. + */ + virtual void + solve() + { + this->solve_level(0); + } + + /** + * @brief Update the fields. This is separate from solve because some derived solvers + * call solve methods from other solvers. + */ + virtual void + update() + {} + + /** + * @brief Update the fields. This is separate from solve because some derived solvers + * call solve methods from other solvers. + */ + virtual void + update_ghosts() + { + solutions.update_ghosts(); + } + + /** + * @brief Print information about the solver to summary.log. + */ + virtual void + print() + {} + + /** + * @brief Set the initial conditions. + */ + void + set_initial_condition() + { + for (const auto &global_index : solve_group.field_indices) + { + if (solver_context->get_user_inputs() + .get_load_initial_condition_parameters() + .get_read_initial_conditions_from_file()) + { + const auto &initial_condition_parameters = + solver_context->get_user_inputs().get_load_initial_condition_parameters(); + for (const auto &initial_condition_file : + initial_condition_parameters.get_initial_condition_files()) + { + auto name_it = + std::find(initial_condition_file.simulation_variable_names.begin(), + initial_condition_file.simulation_variable_names.end(), + attributes_list[global_index].name); + if (name_it != initial_condition_file.simulation_variable_names.end()) + { + dealii::VectorTools::interpolate( + solver_context->get_mapping(), + solver_context->get_dof_manager().get_dof_handler(global_index), + ReadInitialCondition( + *name_it, + attributes_list[global_index].field_type, + initial_condition_file, + solver_context->get_user_inputs().get_spatial_discretization()), + solutions.get_solution_vector(global_index)); + } + } + } + else + { + dealii::VectorTools::interpolate( + solver_context->get_mapping(), + solver_context->get_dof_manager().get_dof_handler(index), + InitialCondition( + index, + attributes_list[global_index].field_type, + solver_context->get_pde_operator()), + solutions.get_solution_vector(global_index)); + } + solutions.apply_initial_condition_for_old_fields(); + } + } + + /** + * @brief Get the invm handler. + */ + [[nodiscard]] const InvmHandler & + get_invm_handler() const + { + return solver_context->get_invm_handler(); + } + + /** + * @brief Get the solution handler. + */ + [[nodiscard]] GroupSolutionHandler & + get_solution_manager() const; + + /** + * @brief Get the element volume container. + */ + [[nodiscard]] const ElementVolumeContainer & + get_element_volume_container() const + { + return solver_context->get_element_volume_container(); + } + + /** + * @brief Get the pde operator. + */ + [[nodiscard]] const std::shared_ptr> & + get_pde_operator() const + { + return solver_context->get_pde_operator(); + } + + /** + * @brief Get the solver context. + */ + [[nodiscard]] const SolverContext & + get_solver_context() const + { + Assert(solver_context != nullptr, dealii::ExcNotInitialized()); + return *solver_context; + } + +protected: + /** + * @brief List of field attributes available. + */ + const std::vector attributes_list; + /** + * @brief Information about the solve group this handler is responsible for. + */ + SolveGroup solve_group; + /** + * @brief Solution vectors for fields handled by this solver. + */ + GroupSolutionHandler solutions; + + std::vector> rhs_operators; + + std::vector *> aux_solvers; + + /** + * @brief Solver context provides access to external information. + */ + const std::shared_ptr> solver_context; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/solvers/mf_operator.h b/include/prismspf/solvers/mf_operator.h new file mode 100644 index 000000000..f36741d82 --- /dev/null +++ b/include/prismspf/solvers/mf_operator.h @@ -0,0 +1,327 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include + +#if DEAL_II_VERSION_MAJOR >= 9 && DEAL_II_VERSION_MINOR >= 7 +# include +# define MATRIX_FREE_OPERATOR_BASE dealii::EnableObserverPointer +#else +# include +# define MATRIX_FREE_OPERATOR_BASE dealii::Subscriptor +#endif + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief This class exists to evaluate a single user-defined operator for the matrix-free + * implementation of some PDE. + * @note Information such as the pde operator are passed in at construction/initialization + * rather than being passed in during function calls because in certain contexts (eg. + * GMG,) the operator gets called by a dealii function, rather than by prismspf, so it + * needs to act as a standalone operator. + * + * @tparam dim The number of dimensions in the problem. + * @tparam degree The polynomial degree of the shape functions. + * @tparam number Datatype to use for `LinearAlgebra::distributed::Vector`. Either + * double or float. + */ +template +class MFOperator : public MATRIX_FREE_OPERATOR_BASE +{ +public: + using BlockVector = SolutionIndexer::BlockVector; + using SolutionVector = SolutionIndexer::SolutionVector; + using ScalarValue = dealii::VectorizedArray; + using VectorValue = dealii::Tensor<1, dim, ScalarValue>; + using MatrixFree = dealii::MatrixFree; + + using Operator = void (PDEOperator::*)( + FieldContainer &, /* variable_list */ + const dealii::Point> &, /* q_point_loc */ + const dealii::VectorizedArray &, /* element_volume */ + unsigned int /* solve_group_id */ + ) const; + + using TensorRank = FieldInfo::TensorRank; + template + using Value = std::conditional_t>; + + template + Value + identity() + { + static Value ident = []() + { + Value obj; + for (int i = 0; i < Value::n_independent_components; ++i) + { + obj[Value::unrolled_to_component_indices(i)] = 1.0; + } + }(); + return ident; + } + + template <> + ScalarValue + identity() + { + static ScalarValue ident(1.0); + return ident; + } + + template + Value + zero() + { + return Value(); + } + + template <> + ScalarValue + zero() + { + static ScalarValue zeroo(0.0); + return zeroo; + } + + /** + * @brief Constructor. + */ + explicit MFOperator(PDEOperator &operator_owner, + Operator oper, + const std::vector &_field_attributes, + const SolutionIndexer &_solution_indexer, + unsigned int _relative_level, + DependencySet _dependency_map) + : MATRIX_FREE_OPERATOR_BASE() + , pde_operator(&operator_owner) + , pde_op(oper) + , field_attributes(_field_attributes) + , solution_indexer(&_solution_indexer) + , relative_level(_relative_level) + , dependency_map(std::move(_dependency_map)) + {} + + /** + * @brief Initialize. + * @note This will look stylistically strange. We pass in the solution handler for the + * fields being solved here, but we don't actually use any of the field data from it. + * This is because the purpose of this MFOperator class is to provide an operator with + * the signature `vmult(VectorType &dst, const VectorType &src)`. Because we are + * using block vectors, for the MFOperator to work, it needs to have access to the index + * mapping in addition to the matrix_free object. Both are stored in the solution + * handler. + */ + void + initialize(const GroupSolutionHandler &dst_solution) + { + solve_group = dst_solution.get_solve_group(); + data = dst_solution.get_matrix_free(relative_level); + field_to_block_index = dst_solution.get_global_to_block_index(); + } + + // public: + /** + * @brief Calls cell_loop on function that calls user-defined operator + */ + void + compute_operator(BlockVector &dst, const BlockVector &src = BlockVector()) const; + +private: + /** + * @brief Calls user-defiend operator + */ + void + compute_local_operator(const MatrixFree &_data, + BlockVector &dst, + const BlockVector &src, + const std::pair &cell_range) const; + +public: + /** + * @brief Compute the diagonal of this operator. + */ + void + compute_diagonal(); + +private: + /** + * @brief Local computation of the diagonal of the operator. + */ + void + compute_local_diagonal(const MatrixFree &_data, + BlockVector &dst, + const unsigned int &dummy, + const std::pair &cell_range) const; + + template + dealii::AlignedVector> + compute_field_diagonal(FieldContainer &variable_list, + DSTContainer &dst_fields, + unsigned int field_index) const; + +public: + /** + * @brief Compute the element volume. (And store in this object? void?) + */ + void + compute_element_volume(); + +private: + /** + * @brief Local computation of the element volume. + */ + void + compute_local_element_volume( + const MatrixFree &_data, + SolutionVector &dst, + const unsigned int &dummy, + const std::pair &cell_range) const; + +public: + /** + * @brief Return the number of DoFs. + */ + dealii::types::global_dof_index + m() const; + + /** + * @brief Return the value of the matrix entry. This function is only valid when row == + * col and when the diagonal is initialized. Additionally, this is only used so that we + * may compile. Trying to use this function will throw an error. + */ + number + el(const unsigned int &row, const unsigned int &col) const; + + /** + * @brief Release all memory and return to state like having called the default + * constructor. + */ + void + clear(); + + /** + * @brief Initialize a given vector with the MatrixFree object that this object + * contains. + */ + void + initialize_dof_vector(SolutionVector &dst, unsigned int dof_handler_index = 0) const; + + /** + * @brief Set constrained entries to one. + */ + void + set_constrained_entries_to_one(SolutionVector &dst) const; + + /** + * @brief Get read access to the MatrixFree object stored with this operator. + */ + std::shared_ptr + get_matrix_free() const; + + /** + * @brief Get read access to the inverse diagonal of this operator. + */ + const std::shared_ptr> & + get_matrix_diagonal_inverse() const; + + /** + * @brief Matrix-vector multiplication. + */ + void + vmult(SolutionVector &dst, const SolutionVector &src) const; + + // NOLINTBEGIN(readability-identifier-naming) + + /** + * @brief Transpose matrix-vector multiplication. + */ + void + Tvmult(SolutionVector &dst, const SolutionVector &src) const; + + // NOLINTEND(readability-identifier-naming) + +private: + /** + * @brief The attribute list of the relevant variables. + */ + std::vector field_attributes; + + /** + * @brief The group being solved + */ + SolveGroup solve_group; + + /** + * @brief PDE operator object (owning class instance of pde_op) for user defined PDEs. + */ + const PDEOperator *pde_operator; + /** + * @brief The actual PDE operator function ptr (eg. compute_rhs) for user defined PDEs. + */ + Operator pde_op; + + /** + * @brief Read-access to fields. + */ + const SolutionIndexer *solution_indexer; + /** + * @brief Level so that correct fields are read from indexer. + */ + unsigned int relative_level; + + /** + * @brief Which fields should be available to the solve. + */ + DependencySet dependency_map; + + /** + * @brief Mapping from field index to block index (only for dst). + */ + std::vector field_to_block_index; + + /** + * @brief Matrix-free object. + */ + std::shared_ptr data; + + /** + * @brief Indices of DoFs on edge in case the operator is used in GMG context. + */ + std::vector> edge_constrained_indices; + + /** + * @brief The diagonal matrix. + */ + std::shared_ptr> diagonal_entries; + + /** + * @brief The inverse diagonal matrix. + */ + std::shared_ptr> inverse_diagonal_entries; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/solvers/solve_context.h b/include/prismspf/solvers/solve_context.h new file mode 100644 index 000000000..e5f633bf9 --- /dev/null +++ b/include/prismspf/solvers/solve_context.h @@ -0,0 +1,147 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +/** + * @brief This class provides context for a solver with ptrs to all the relevant + * dependencies. + * + * The context in this case refers to all the finite element machinery needed to solve the + * fields. For example, it contains the triangulation manager and pde operator that + * evaluates the user-specified PDEs. + */ +template +class SolveContext +{ +public: + /** + * @brief Constructor. + */ + SolveContext(const UserInputParameters &_user_inputs, + const TriangulationManager &_triangulation_manager, + const ConstraintManager &_constraint_manager, + const DofManager &_dof_manager, + SolutionIndexer &_solution_indexer, + std::shared_ptr> _pde_operator) + : user_inputs(&_user_inputs) + , triangulation_manager(&_triangulation_manager) + , constraint_manager(&_constraint_manager) + , dof_manager(&_dof_manager) + , solution_indexer(&_solution_indexer) + , pde_operator(_pde_operator) {}; + + /** + * @brief Get the user-inputs. + */ + [[nodiscard]] const UserInputParameters & + get_user_inputs() const + { + Assert(user_inputs != nullptr, dealii::ExcNotInitialized()); + return *user_inputs; + } + + /** + * @brief Get the triangulation manager. + */ + [[nodiscard]] const TriangulationManager & + get_triangulation_manager() const + { + Assert(triangulation_manager != nullptr, dealii::ExcNotInitialized()); + return *triangulation_manager; + } + + /** + * @brief Get the constraint manager. + */ + [[nodiscard]] const ConstraintManager & + get_constraint_manager() const + { + Assert(constraint_manager != nullptr, dealii::ExcNotInitialized()); + return *constraint_manager; + } + + /** + * @brief Get the dof manager. + */ + [[nodiscard]] const DofManager & + get_dof_manager() const + { + Assert(dof_manager != nullptr, dealii::ExcNotInitialized()); + return *dof_manager; + } + + /** + * @brief Get the solution manager. + */ + [[nodiscard]] SolutionIndexer & + get_solution_indexer() const + { + Assert(solution_indexer != nullptr, dealii::ExcNotInitialized()); + return *solution_indexer; + } + + /** + * @brief Get a shared pointer to the pde operator. + */ + [[nodiscard]] const std::shared_ptr> & + get_pde_operator() const + { + Assert(pde_operator != nullptr, dealii::ExcNotInitialized()); + return pde_operator; + } + +private: + /** + * @brief User-inputs. + */ + const UserInputParameters *user_inputs; + + /** + * @brief Triangulation manager. + */ + const TriangulationManager *triangulation_manager; + + /** + * @brief Constraint manager. + */ + const ConstraintManager *constraint_manager; + + /** + * @brief DoF manager. + */ + const DofManager *dof_manager; + + /** + * @brief Solution manager. + */ + SolutionIndexer *solution_indexer; + + /** + * @brief PDE operator. + */ + std::shared_ptr> pde_operator; + + /** + * @brief PDE operator for float precision. + */ + std::shared_ptr> pde_operator_float; +}; + +PRISMS_PF_END_NAMESPACE diff --git a/include/prismspf/solvers/solvers.h b/include/prismspf/solvers/solvers.h new file mode 100644 index 000000000..0f8cb1cab --- /dev/null +++ b/include/prismspf/solvers/solvers.h @@ -0,0 +1,8 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once + +#include +#include +#include diff --git a/include/prismspf/utilities/integrator.h b/include/prismspf/utilities/integrator.h index 10e525752..7aed30fa9 100644 --- a/include/prismspf/utilities/integrator.h +++ b/include/prismspf/utilities/integrator.h @@ -6,6 +6,8 @@ #include #include +#include + #include #include @@ -41,6 +43,63 @@ class Integrator compute_integral(std::vector &integral_value, const dealii::DoFHandler &dof_handler, const VectorType &vector) const; + + using TensorRank = FieldInfo::TensorRank; + template + using Value = std::conditional_t>; + + template + static Value + integrate(const dealii::DoFHandler &dof_handler, const VectorType &vector) + { + [[maybe_unused]] constexpr unsigned int expected_components = + dealii::Tensor::n_independent_components; + Assert(dof_handler.get_fe().n_components() == expected_components, + dealii::ExcMessage("The provided DoFHandler does not have the same number of " + "components as the expected ones. For scalar fields there " + "should be 1 component.")); + + // Update ghosts + vector.update_ghost_values(); + + // Set quadrature rule and FEValues to update the JxW values + const dealii::QGaussLobatto quadrature(degree + 1); + dealii::FEValues fe_values(dof_handler.get_fe(), + quadrature, + dealii::update_values | dealii::update_JxW_values); + + // Get the number of quadrature points + const unsigned int num_quad_points = quadrature.size(); + + // Create a value vector + std::vector> quad_values(num_quad_points); + + // Loop over the cells provided by the DoFHandler + Value value(0); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + // Reinitialize the cell + fe_values.reinit(cell); + + // Get the values + fe_values.get_function_values(vector, quad_values); + + // Sum up the product of the JxW and values at each quadrature point to + // compute the element integral. + for (unsigned int q_point = 0; q_point < num_quad_points; ++q_point) + { + value += quad_values[q_point] * fe_values.JxW(q_point); + } + } + } + + Value integral = dealii::Utilities::MPI::sum(value, MPI_COMM_WORLD); + return integral; + } }; PRISMS_PF_END_NAMESPACE diff --git a/src/core/constraint_manager.cc b/src/core/constraint_manager.cc new file mode 100644 index 000000000..132534093 --- /dev/null +++ b/src/core/constraint_manager.cc @@ -0,0 +1,466 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include "prismspf/core/dof_manager.h" + +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +ConstraintManager::ConstraintManager( + const std::vector &field_attributes, + const std::set &solve_groups, + const DofManager &_dof_manager, + const PDEOperator &_pde_operator) + : dof_manager(_dof_manager) + , pde_operator(_pde_operator) + , constraints(field_attributes.size()) +{ + for (const auto &solve_group : solve_groups) + { + /* TODO (fractalsbyx) figure out mg depth. Careful. Aux fields inherit this from + * primary fields, as well as any dependencies*/ + unsigned int num_mg_levels = 1; + for (const auto &field_index : solve_group.field_indices) + { + constraints[field_index].resize(num_mg_levels, nullptr); + for (unsigned int relative_level = 0; relative_level < num_mg_levels; + relative_level++) + { + const dealii::DoFHandler &dof_handler = + dof_manager->get_dof_handler(field_index, relative_level); + constraints[field_index][relative_level] = + std::make_shared>( + dof_handler.locally_owned_dofs(), + dealii::DoFTools::extract_locally_relevant_dofs(dof_handler)); + } + // TODO (fractalsbyx): construct change_constraints + } + } +} + +template +std::vector *> +ConstraintManager::get_constraints( + const std::set &field_indices, + unsigned int relative_level) const +{ + std::vector *> selected_constraints; + selected_constraints.reserve(field_indices.size()); + for (const auto index : field_indices) + { + selected_constraints.push_back(constraints[index][relative_level].get()); + } + return selected_constraints; +} + +template +const dealii::AffineConstraints & +ConstraintManager::get_constraint(Types::Index index, + unsigned int relative_level) const +{ + return *constraints[index][relative_level]; +} + +template +std::vector *> +ConstraintManager::get_change_constraints( + const std::set &field_indices, + unsigned int relative_level) const +{ + std::vector *> selected_constraints; + selected_constraints.reserve(field_indices.size()); + for (const auto index : field_indices) + { + selected_constraints.push_back(change_constraints[index][relative_level].get()); + } + return selected_constraints; +} + +template +const dealii::AffineConstraints & +ConstraintManager::get_change_constraint( + Types::Index index, + unsigned int relative_level) const +{ + return *change_constraints[index][relative_level]; +} + +template +void +ConstraintManager::reinit(const dealii::Mapping &mapping) +{ + for (unsigned int field_index = 0; field_index < constraints.size(); ++field_index) + { + for (unsigned int relative_level = 0; + relative_level < constraints[field_index].size(); + ++relative_level) + { + std::shared_ptr> &constraint = + constraints[field_index][relative_level]; + const dealii::DoFHandler &dof_handler = + dof_manager->get_dof_handler(field_index, relative_level); + constraint->reinit(dof_handler.locally_owned_dofs(), + dealii::DoFTools::extract_locally_relevant_dofs( + dof_handler)); + + // 1. Make hanging node constraints + dealii::DoFTools::make_hanging_node_constraints(dof_handler, *constraint); + // 2. Make boundary constraints + const BoundaryCondition &boundary_condition = + user_inputs->get_boundary_parameters().get_boundary_condition_list().at( + field_index); + make_bc_constraints(mapping, dof_handler, boundary_condition); + // 3. Make pinned point constraints + if (user_inputs->get_boundary_parameters().has_pinned_point(index)) + { + const auto &[value, target_point] = + user_inputs->get_boundary_parameters().get_pinned_point(index); + set_pinned_point(constraint, target_point, value, dof_handler, false); + } + constraint.close(); + } + for (unsigned int relative_level = 0; + relative_level < change_constraints[field_index].size(); + ++relative_level) + { + std::shared_ptr> &constraint = + change_constraints[field_index][relative_level]; + const dealii::DoFHandler &dof_handler = + dof_manager->get_dof_handler(field_index, relative_level); + constraint->reinit(dof_handler.locally_owned_dofs(), + dealii::DoFTools::extract_locally_relevant_dofs( + dof_handler)); + + // 1. Make hanging node constraints + dealii::DoFTools::make_hanging_node_constraints(dof_handler, *constraint); + // 2. Make boundary constraints + const std::map &boundary_condition = + user_inputs->get_boundary_parameters().get_boundary_condition_list().at( + field_index); + make_bc_constraints(mapping, dof_handler, boundary_condition, true); + constraint.close(); + } + } +} + +template +void +ConstraintManager::make_bc_constraints( + dealii::AffineConstraints &constraint, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + const std::map &boundary_condition, + const FieldInfo::TensorRank tensor_rank, + Types::Index field_index, + bool for_change_term) +{ + for (const auto &[component, condition] : boundary_condition) + { + for (const auto &[boundary_id, boundary_type] : + condition.get_boundary_condition_map()) + { + const number dirichlet_value = condition.get_dirichlet_value(boundary_id); + make_one_boundary_constraint(constraint, + boundary_id, + component, + boundary_type, + dirichlet_value, + mapping, + dof_handler, + tensor_rank, + field_index, + for_change_term); + } + } +} + +template +void +ConstraintManager::make_one_boundary_constraint( + dealii::AffineConstraints &_constraints, + const unsigned int boundary_id, + const unsigned int component, + BoundaryCondition::Type boundary_type, + const number dirichlet_value, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + FieldInfo::TensorRank tensor_rank, + Types::Index field_index, + bool for_change_term) const +{ + const bool is_vector_field = tensor_rank == FieldInfo::TensorRank::Vector; + const dealii::ComponentMask mask = + is_vector_field ? vector_component_mask[component] : scalar_empty_mask; + + // Apply the boundary conditions + switch (boundary_type) + { + case BoundaryCondition::Type::Natural: + { + make_natural_constraints(); + break; + } + case BoundaryCondition::Type::Dirichlet: + { + make_uniform_dirichlet_constraints(_constraints, + mapping, + dof_handler, + boundary_id, + is_vector_field, + for_change_term ? 0.0 : dirichlet_value, + + mask); + break; + } + case BoundaryCondition::Type::NonuniformDirichlet: + case BoundaryCondition::Type::TimeDependentNonuniformDirichlet: + { + make_nonuniform_dirichlet_constraints(_constraints, + mapping, + dof_handler, + boundary_id, + field_index, + is_vector_field, + + mask, + for_change_term); + break; + } + case BoundaryCondition::Type::Periodic: + { + // Skip boundary ids that are odd since those map to the even faces + if (boundary_id % 2 != 0) + { + break; + } + make_periodic_constraints(dof_handler, boundary_id, _constraints, mask); + break; + } + case BoundaryCondition::Type::Neumann: + { + Assert(false, FeatureNotImplemented("Neumann boundary conditions")); + break; + } + case BoundaryCondition::Type::NonuniformNeumann: + { + Assert(false, FeatureNotImplemented("Nonuniform neumann boundary conditions")); + break; + } + case BoundaryCondition::Type::TimeDependentNonuniformNeumann: + { + Assert(false, + FeatureNotImplemented( + "Time dependent nonuniform neumann boundary conditions")); + break; + } + default: + { + AssertThrow(false, UnreachableCode()); + } + } +} + +template +void +ConstraintManager::update_time_dependent_constraints( + const dealii::Mapping &mapping, + const std::vector *> &dof_handlers) +{ + for (const auto &[index, variable] : user_inputs->get_variable_attributes()) + { + // Check that we have time-dependent constraints before recreating the whole + // constraint set. + if (user_inputs->get_boundary_parameters().is_time_dependent(index)) + { + // TODO (landinjm): Is there a way to update the constraint set without + // recreating + // it? + make_constraint(mapping, *dof_handlers.at(index), index); + } + } +} + +template +const std::array + ConstraintManager::vector_component_mask = []() +{ + std::array masks {}; + for (unsigned int i = 0; i < dim; ++i) + { + dealii::ComponentMask temp_mask(dim, false); + temp_mask.set(i, true); + masks[i] = temp_mask; + } + return masks; +}(); + +template +const dealii::ComponentMask ConstraintManager::scalar_empty_mask {}; + +template +void +ConstraintManager::make_natural_constraints() const +{ + // Do nothing because they are naturally enforced. +} + +template +void +ConstraintManager::make_uniform_dirichlet_constraints( + dealii::AffineConstraints &_constraints, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + const unsigned int &boundary_id, + const bool &is_vector_field, + const number &value, + + const dealii::ComponentMask &mask) const +{ + dealii::VectorTools::interpolate_boundary_values( + mapping, + dof_handler, + boundary_id, + dealii::Functions::ConstantFunction(value, is_vector_field ? dim : 1), + _constraints, + mask); +} + +template +void +ConstraintManager::make_nonuniform_dirichlet_constraints( + dealii::AffineConstraints &_constraints, + const dealii::Mapping &mapping, + const dealii::DoFHandler &dof_handler, + const unsigned int &boundary_id, + const unsigned int &field_index, + const bool &is_vector_field, + + const dealii::ComponentMask &mask, + bool is_change_term) const +{ + if (!is_change_term) + { + dealii::VectorTools::interpolate_boundary_values( + mapping, + dof_handler, + boundary_id, + NonuniformDirichlet(field_index, + boundary_id, + pde_operator, + is_vector_field ? dim : 1), + _constraints, + mask); + } + else + { + dealii::VectorTools::interpolate_boundary_values( + mapping, + dof_handler, + boundary_id, + dealii::Functions::ZeroFunction(is_vector_field ? dim : 1), + _constraints, + mask); + } +} + +template +void +ConstraintManager::make_periodic_constraints( + dealii::AffineConstraints &_constraints, + const dealii::DoFHandler &dof_handler, + const unsigned int &boundary_id, + const dealii::ComponentMask &mask) const +{ + // Create a vector of matched pairs that we fill and enforce upon the + // constaints + std::vector< + dealii::GridTools::PeriodicFacePair::cell_iterator>> + periodicity_vector; + + // Determine the direction + const unsigned int direction = boundary_id / 2; + + // Collect the matched pairs on the coarsest level of the mesh + dealii::GridTools::collect_periodic_faces(dof_handler, + boundary_id, + boundary_id + 1, + direction, + periodicity_vector); + + // Set constraints + dealii::DoFTools::make_periodicity_constraints(periodicity_vector, + _constraints, + mask); +} + +template +void +ConstraintManager::set_pinned_point( + dealii::AffineConstraints &constraint, + const dealii::Point &target_point, + const std::array &value, + const dealii::DoFHandler &dof_handler, + const FieldInfo::TensorRank tensor_rank, + bool is_change_term) const +{ + const double tolerance = 1.0e-2; + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (!cell->is_locally_owned()) + { + continue; + } + + for (unsigned int vertex = 0; vertex < dealii::GeometryInfo::vertices_per_cell; + ++vertex) + { + const auto vertex_point = cell->vertex(vertex); + if (target_point.distance(vertex_point) >= tolerance * cell->diameter()) + { + continue; + } + + const unsigned int dof_index = cell->vertex_dof_index(vertex, 0); + const unsigned int dimension = + tensor_rank == FieldInfo::TensorRank::Vector ? dim : 1; + for (unsigned int component = 0; component < dimension; ++component) + { + constraint.add_line(dof_index + component); + constraint.set_inhomogeneity(dof_index + component, + is_change_term ? value[component] : 0.0); + } + } + } +} + +// #include "core/constraint_managser.inst" + +PRISMS_PF_END_NAMESPACE diff --git a/src/core/dof_manager.cc b/src/core/dof_manager.cc new file mode 100644 index 000000000..6055a2363 --- /dev/null +++ b/src/core/dof_manager.cc @@ -0,0 +1,70 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +#include +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +DofManager::DofManager(const std::vector &field_attributes) +{ + unsigned int num_mg_levels = 1; // Todo: upgrade to multigrid + level_dof_handlers.resize(num_mg_levels); + dof_handlers.resize(field_attributes.size()); + for (unsigned int field_index = 0; field_index < field_attributes.size(); ++field_index) + { + dof_handlers[field_index].resize(num_mg_levels, nullptr); + for (unsigned int relative_level = 0; relative_level < num_mg_levels; + ++relative_level) + { + dof_handlers[field_index][relative_level] = + &level_dof_handlers[relative_level][field_attributes[field_index].field_type]; + } + } +} + +template +void +DofManager::init(const TriangulationManager &triangulation_handler) +{ + for (unsigned int relative_level = 0; relative_level < level_dof_handlers.size(); + ++relative_level) + { + for (unsigned int rank = 0; rank < 2; ++rank) + { + dealii::DoFHandler &dof_handler = level_dof_handlers[relative_level][rank]; + dof_handler.reinit(triangulation_handler.get_triangulation(relative_level)); + dof_handler.distribute_dofs(SystemWide::fe_systems.at(rank)); + } + } +} + +template +void +DofManager::reinit(const TriangulationManager &triangulation_handler) +{ + init(triangulation_handler); +} + +template +const std::vector *> & +DofManager::get_dof_handlers(const std::set &field_indices, + unsigned int relative_level) const +{ + std::vector *> selected_dof_handlers; + selected_dof_handlers.reserve(field_indices.size()); + for (const auto index : field_indices) + { + selected_dof_handlers.push_back(dof_handlers[index][relative_level].get()); + } + return selected_dof_handlers; +} + +// #include "core/dof_manager.inst" + +PRISMS_PF_END_NAMESPACE diff --git a/src/core/dst_container.cc b/src/core/dst_container.cc new file mode 100644 index 000000000..2a38853cd --- /dev/null +++ b/src/core/dst_container.cc @@ -0,0 +1,105 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +DSTContainer::DSTContainer( + const std::set &_field_indices, + const std::vector &_field_attributes, + const MatrixFree &matrix_free, + const std::vector &_field_to_block_index) + : field_attributes_ptr(&_field_attributes) + , field_indices(&_field_indices) +{ + const std::vector field_attributes; + // Initialize the feeval vectors + feeval_scalar.clear(); + feeval_vector.clear(); + feeval_scalar.resize(field_attributes.size()); + feeval_vector.resize(field_attributes.size()); + integration_flags = std::vector(field_attributes.size(), EvalFlags::nothing); + for (const auto &field_index : *field_indices) + { + const unsigned int block_index = _field_to_block_index[field_index]; + if (field_attributes[field_index].field_type == FieldInfo::TensorRank::Scalar) + { + feeval_scalar[field_index] = std::make_unique(matrix_free, block_index); + } + else + { + feeval_vector[field_index] = std::make_unique(matrix_free, block_index); + } + } +} + +template +void +DSTContainer::reinit(unsigned int cell) +{ + for (const unsigned int &field_index : *field_indices) + { + if ((*field_attributes_ptr)[field_index].field_type == TensorRank::Scalar) + { + feeval_scalar[field_index].reinit(cell); + } + else + { + feeval_vector[field_index].reinit(cell); + } + } +} + +template +void +DSTContainer::integrate(unsigned int field_index) +{ + if ((*field_attributes_ptr)[field_index].field_type == TensorRank::Scalar) + { + feeval_scalar[field_index]->integrate(integration_flags[field_index]); + } + else /* vector */ + { + feeval_vector[field_index]->integrate(integration_flags[field_index]); + } +} + +template +void +DSTContainer::integrate_and_distribute(unsigned int field_index, + SolutionVector &dst) +{ + const std::vector field_attributes = *field_attributes_ptr; + + if (field_attributes[field_index].field_type == TensorRank::Scalar) + { + feeval_scalar[field_index]->integrate_scatter(integration_flags[field_index], dst); + } + else /* vector */ + { + feeval_vector[field_index]->integrate_scatter(integration_flags[field_index], dst); + } +} + +// #include "core/field_container.inst" + +PRISMS_PF_END_NAMESPACE diff --git a/src/core/field_container.cc b/src/core/field_container.cc new file mode 100644 index 000000000..3645beaf5 --- /dev/null +++ b/src/core/field_container.cc @@ -0,0 +1,183 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include "prismspf/core/dependencies.h" +#include "prismspf/core/field_attributes.h" +#include "prismspf/core/solution_indexer.h" +#include "prismspf/core/solve_group.h" + +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +FieldContainer::FieldContainer( + const std::vector &_field_attributes, + SolutionIndexer &_solution_indexer, + unsigned int _relative_level, + const DependencySet &dependency_map) + : field_attributes_ptr(&_field_attributes) + , relative_level(_relative_level) +{ + const std::vector field_attributes; + // Initialize the feeval vectors + feeval_deps_scalar.clear(); + feeval_deps_vector.clear(); + feeval_deps_scalar.resize(field_attributes.size()); + feeval_deps_vector.resize(field_attributes.size()); + + for (const auto &[field_index, dependency] : dependency_map) + { + const auto mf_id_pair = + solution_indexer->get_solution_level_and_block_index(field_index, relative_level); + if (field_attributes[field_index].field_type == FieldInfo::TensorRank::Scalar) + { + feeval_deps_scalar[field_index] = {dependency, mf_id_pair}; + } + else if (field_attributes[field_index].field_type == FieldInfo::TensorRank::Vector) + { + feeval_deps_vector[field_index] = {dependency, mf_id_pair}; + } + // else Assert unreachable + } + MatrixFree matrix_free; + unsigned int block_index = 0; + for (const auto &field_index : solve_group->field_indices) + { + if (field_attributes[field_index].field_type == FieldInfo::TensorRank::Scalar) + { + dst_feeval_scalar[field_index] = std::make_unique(matrix_free, block_index); + } + else if (field_attributes[field_index].field_type == FieldInfo::TensorRank::Vector) + { + dst_feeval_vector[field_index] = std::make_unique(matrix_free, block_index); + } + } +} + +template +unsigned int +FieldContainer::get_n_q_points() const +{} + +template +dealii::Point::ScalarValue> +FieldContainer::get_q_point_location() const +{ + Types::Index field_index = *(solve_group->field_indices.begin()); + const std::vector &field_attributes = *field_attributes_ptr; + if (field_attributes[field_index].field_type == TensorRank::Scalar) + { + return dst_feeval_scalar[field_index]->quadrature_point(q_point); + } + /*else if vector */ + { + return dst_feeval_vector[field_index]->quadrature_point(q_point); + } + /* unreachable */ + return dealii::Point(); +} + +template +void +FieldContainer::reinit(unsigned int cell) +{ + const DependencySet &dependency_map; // rhs or lhs + for (const auto &[field_index, dependency] : dependency_map) + { + feeval_deps_scalar[field_index].reinit(cell); + } +} + +template +void +FieldContainer::eval() +{ + const DependencySet &dependency_map; // rhs or lhs + for (const auto &[field_index, dependency] : dependency_map) + { + feeval_deps_scalar[field_index].eval(); + } +} + +template +void +FieldContainer::reinit_and_eval(unsigned int cell) +{ + const DependencySet &dependency_map; // rhs or lhs + for (const auto &[field_index, dependency] : dependency_map) + { + feeval_deps_scalar[field_index].reinit_and_eval(cell); + } +} + +template +void +FieldContainer::integrate() +{ + const std::vector field_attributes = *field_attributes_ptr; + for (const Types::Index &field_index : solve_group->field_indices) + { + const EvalFlags submission_flags = equation_type == EquationType::RHS + ? field_attributes[field_index].eval_flags_rhs + : field_attributes[field_index].eval_flags_lhs; + if (field_attributes[field_index].field_type == TensorRank::Scalar) + { + dst_feeval_scalar[field_index]->integrate(submission_flags); + } + else /* vector */ + { + dst_feeval_vector[field_index]->integrate(submission_flags); + } + } + /* AssertThrow(false, + dealii::ExcMessage( + "Integrate called for a solve type that is not NonexplicitLHS.")); */ +} + +template +void +FieldContainer::integrate_and_distribute() +{ + const std::vector field_attributes = *field_attributes_ptr; + for (const Types::Index &field_index : solve_group->field_indices) + { + const EvalFlags submission_flags = equation_type == EquationType::RHS + ? field_attributes[field_index].eval_flags_rhs + : field_attributes[field_index].eval_flags_lhs; + SolutionVector &dst = equation_type == EquationType::RHS + ? solutions.new_solutions(field_index, relative_level) + : solutions.change_solutions(field_index, relative_level); + ; + if (field_attributes[field_index].field_type == TensorRank::Scalar) + { + dst_feeval_scalar[field_index]->integrate_scatter(submission_flags, dst); + } + else /* vector */ + { + dst_feeval_vector[field_index]->integrate_scatter(submission_flags, dst); + } + } +} + +// #include "core/field_container.inst" + +PRISMS_PF_END_NAMESPACE diff --git a/src/core/group_solution_handler.cc b/src/core/group_solution_handler.cc new file mode 100644 index 000000000..fe8e68207 --- /dev/null +++ b/src/core/group_solution_handler.cc @@ -0,0 +1,384 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +#include +#include +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +GroupSolutionHandler::GroupSolutionHandler( + const SolveGroup &_solve_group, + const std::vector &_attributes_list) +{ + block_to_global_index.assign(_solve_group.field_indices.begin(), + _solve_group.field_indices.end()); + global_to_block_index = + std::vector(_attributes_list.size(), Numbers::invalid_index); + for (unsigned int i = 0; i < block_to_global_index.size(); ++i) + { + global_to_block_index[block_to_global_index[i]] = i; + } +} + +template +auto +GroupSolutionHandler::get_solution_full_vector(unsigned int relative_level) + -> BlockVector & +{ + return solution_levels[relative_level].solutions; +} + +template +auto +GroupSolutionHandler::get_solution_vector(unsigned int global_index, + unsigned int relative_level) + -> SolutionVector & +{ + return solution_levels[relative_level].solutions.block( + global_to_block_index[global_index]); +} + +template +auto +GroupSolutionHandler::get_old_solution_full_vector( + unsigned int age, + unsigned int relative_level) -> BlockVector & +{ + return solution_levels[relative_level].old_solutions[age]; +} + +template +auto +GroupSolutionHandler::get_old_solution_vector(unsigned int age, + unsigned int global_index, + unsigned int relative_level) + -> SolutionVector & +{ + return solution_levels[relative_level].old_solutions[age].block( + global_to_block_index[global_index]); +} + +template +auto +GroupSolutionHandler::get_new_solution_full_vector( + unsigned int relative_level) -> BlockVector & +{ + return solution_levels[relative_level].new_solutions; +} + +template +auto +GroupSolutionHandler::get_new_solution_vector(unsigned int global_index, + unsigned int relative_level) + -> SolutionVector & +{ + return solution_levels[relative_level].new_solutions.block( + global_to_block_index[global_index]); +} + +template +SolutionLevel & +GroupSolutionHandler::get_solution_level(unsigned int relative_level) +{ + return solution_levels[relative_level]; +} + +template +auto +GroupSolutionHandler::get_matrix_free(unsigned int relative_level) + -> MatrixFree & +{ + return solution_levels[relative_level].matrix_free; +} + +template +auto +GroupSolutionHandler::get_solution_full_vector( + unsigned int relative_level) const -> const BlockVector & +{ + return solution_levels[relative_level].solutions; +} + +template +auto +GroupSolutionHandler::get_solution_vector(unsigned int global_index, + unsigned int relative_level) const + -> const SolutionVector & +{ + return solution_levels[relative_level].solutions.block( + global_to_block_index[global_index]); +} + +template +auto +GroupSolutionHandler::get_old_solution_full_vector( + unsigned int age, + unsigned int relative_level) const -> const BlockVector & +{ + return solution_levels[relative_level].old_solutions[age]; +} + +template +auto +GroupSolutionHandler::get_old_solution_vector( + unsigned int age, + unsigned int global_index, + unsigned int relative_level) const -> const SolutionVector & +{ + return solution_levels[relative_level].old_solutions[age].block( + global_to_block_index[global_index]); +} + +template +auto +GroupSolutionHandler::get_new_solution_full_vector( + unsigned int relative_level) const -> const BlockVector & +{ + return solution_levels[relative_level].new_solutions; +} + +template +auto +GroupSolutionHandler::get_new_solution_vector( + unsigned int global_index, + unsigned int relative_level) const -> const SolutionVector & +{ + return solution_levels[relative_level].new_solutions.block( + global_to_block_index[global_index]); +} + +template +const SolutionLevel & +GroupSolutionHandler::get_solution_level(unsigned int relative_level) const +{ + return solution_levels[relative_level]; +} + +template +auto +GroupSolutionHandler::get_matrix_free(unsigned int relative_level) const + -> const MatrixFree & +{ + return solution_levels[relative_level].matrix_free; +} + +template +unsigned int +GroupSolutionHandler::get_block_index(unsigned int global_index) const +{ + return global_to_block_index[global_index]; +} + +template +auto +GroupSolutionHandler::get_solve_group() const -> const SolveGroup & +{ + return solve_group; +} + +template +const std::vector & +GroupSolutionHandler::get_global_to_block_index() const +{ + return global_to_block_index; +} + +// TODO (fractalsbyx): This might all need to go in reinit(). Check if dof_handler and +// constraint ptrs change. +template +template +void +GroupSolutionHandler::init( + const dealii::Mapping &mapping, + const DofManager &dof_manager, + const ConstraintManager &constraint_manager, + const dealii::Quadrature &quad) +{ + ConditionalOStreams::pout_base() + << "Initializing solution set for solver " << solve_group.id << "...\n" + << std::flush; + Timer::start_section("reinitialize solution set"); + + // Initialize matrixfree objects + // TODO + for (unsigned int relative_level = 0; relative_level < solution_levels.size(); + ++relative_level) + { + SolutionLevel &solution_level = solution_levels[relative_level]; + MatrixFree &matrix_free = solution_level.matrix_free; + matrix_free.reinit( + mapping, + dof_manager.get_dof_handlers(solve_group.field_indices, relative_level), + constraint_manager.get_constraints(solve_group.field_indices, relative_level), + quad); + } + // Initialize solution vectors + reinit(); + + Timer::end_section("reinitialize solution set"); +} + +template +void +GroupSolutionHandler::reinit() +{ + for (auto &solution_level : solution_levels) + { + BlockVector &solutions = solution_level.solutions; + BlockVector &new_solutions = solution_level.new_solutions; + std::array &old_solutions = + solution_level.old_solutions; + MatrixFree &matrix_free = solution_level.matrix_free; + + // These partitioners basically just provide the number of elements in a distributed + // way + std::vector> + partitioners; + partitioners.reserve(solve_group.field_indices.size()); + for (unsigned int block_index = 0; block_index < solve_group.field_indices.size(); + ++block_index) + { + partitioners.push_back(matrix_free.get_vector_partitioner(block_index)); + } + // TODO (fractalsbyx): Check that the default MPI communicator is correct here + solutions.reinit(partitioners); + new_solutions.reinit(partitioners); + for (unsigned int i = 0; i < Numbers::max_saved_increments; ++i) + { + if (i < oldest_saved) + { + old_solutions[i].reinit(partitioners); + } + } + } +} + +// TODO (fractalsbyx): Check if this is necessary to repeat. Check if dof_handler ptrs +// change. +template +void +GroupSolutionHandler::init_solution_transfer() +{ + const dealii::DoFHandler &dof_handler = + solution_levels[0].matrix_free.get_dof_handler(); + solution_transfer = SolutionTransfer(dof_handler); + for (unsigned int i = 0; i < oldest_saved; ++i) + { + old_solution_transfer[i] = SolutionTransfer(dof_handler); + } +} + +template +void +GroupSolutionHandler::prepare_for_solution_transfer() +{ + solution_transfer.prepare_for_coarsening_and_refinement(solution_levels[0].solutions); + for (unsigned int i = 0; i < oldest_saved; ++i) + { + old_solution_transfer[i].prepare_for_coarsening_and_refinement( + solution_levels[0].old_solutions[i]); + } +} + +template +void +GroupSolutionHandler::execute_solution_transfer() +{ + solution_transfer.interpolate(solution_levels[0].solutions); + for (unsigned int i = 0; i < oldest_saved; ++i) + { + old_solution_transfer[i].interpolate(solution_levels[0].old_solutions[i]); + } +} + +// TODO (fractalsbyx): Check if this is necessary for all solutions +template +void +GroupSolutionHandler::update_ghosts(unsigned int relative_level) const +{ + solution_levels[relative_level].solutions.update_ghost_values(); + solution_levels[relative_level].new_solutions.update_ghost_values(); + for (unsigned int i = 0; i < oldest_saved; ++i) + { + solution_levels[relative_level].old_solutions[i].update_ghost_values(); + } +} + +// TODO (fractalsbyx): Check if this is necessary for all solutions +template +void +GroupSolutionHandler::zero_out_ghosts(unsigned int relative_level) const +{ + solution_levels[relative_level].solutions.zero_out_ghost_values(); + solution_levels[relative_level].new_solutions.zero_out_ghost_values(); + for (unsigned int i = 0; i < oldest_saved; ++i) + { + solution_levels[relative_level].old_solutions[i].zero_out_ghost_values(); + } +} + +template +void +GroupSolutionHandler::apply_constraints(unsigned int relative_level) +{ + BlockVector &solutions = solution_levels[relative_level].solutions; + BlockVector &new_solutions = solution_levels[relative_level].new_solutions; + MatrixFree &matrix_free = solution_levels[relative_level].matrix_free; + + for (unsigned int i = 0; i < matrix_free.get_affine_constraints().size(); ++i) + { + matrix_free.get_affine_constraints()[i].distribute(solutions.block(i)); + matrix_free.get_affine_constraints()[i].distribute(new_solutions.block(i)); + // TODO: Check if this needs to be done to both or just one of the above + } +} + +// TODO (fractalsbyx): Replace into initial_conditions module +template +void +GroupSolutionHandler::apply_initial_condition_for_old_fields() +{ + for (auto &solution_level : solution_levels) + { + BlockVector &solutions = solution_level.solutions; + std::array &old_solutions = + solution_level.old_solutions; + + for (unsigned int i = 0; i < oldest_saved; ++i) + { + old_solutions[i] = solutions; + } + } +} + +template +void +GroupSolutionHandler::update(unsigned int relative_level) +{ + // bubble-swap method. bubble the discarded solution up to 'new_solution' + SolutionLevel &solution_level = solution_levels[relative_level]; + for (int age = oldest_saved - 1; age >= 0; --age) + { + if (age > 0) + { + solution_level.old_solutions[age].swap(solution_level.old_solutions[age - 1]); + } + else + { + solution_level.old_solutions[age].swap(solution_level.solutions); + } + } + solution_level.solutions.swap(solution_level.new_solutions); +} + +// #include "core/group_solution_handler.inst" + +PRISMS_PF_END_NAMESPACE diff --git a/src/core/problem.cc b/src/core/problem.cc new file mode 100644 index 000000000..65fcbe304 --- /dev/null +++ b/src/core/problem.cc @@ -0,0 +1,352 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#pragma once +#include +#include + +#include "prismspf/core/simulation_timer.h" +#include "prismspf/user_inputs/temporal_discretization.h" +#include "prismspf/user_inputs/user_input_parameters.h" + +PRISMS_PF_BEGIN_NAMESPACE + +template +std::set> +get_solution_managers_from_solvers( + const std::vector>> &solvers) +{ + // Todo: upgrade to recursive for aux solvers + std::set> solution_managers; + for (const auto &solver : solvers) + { + solution_managers.insert(solver->get_solution_manager()); + } + return solution_managers; +} + +template +Problem::Problem( + const std::vector &_field_attributes, + const std::vector &_solve_groups, + const UserInputParameters &_user_inputs, + PhaseFieldTools &_pf_tools, + const std::shared_ptr> &_pde_operator) + : field_attributes(_field_attributes) + , solve_groups(_solve_groups) + , user_inputs_ptr(&_user_inputs) + , pf_tools(&_pf_tools) + , triangulation_manager(false) + , dof_manager(field_attributes) + , constraint_manager(field_attributes, solve_groups, dof_manager, _pde_operator) + , solvers(solve_groups.size(), nullptr) + , solution_indexer(get_solution_managers_from_solvers(solvers)) + , solve_context(_user_inputs, + triangulation_manager, + constraint_manager, + dof_manager, + solution_indexer, + _pde_operator) + , grid_refiner(_user_inputs, triangulation_manager) +{} + +template +void +Problem::init_system() +{ + const UserInputParameters &user_inputs = *user_inputs_ptr; + const unsigned int n_proc = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + const unsigned int n_vect_doubles = dealii::VectorizedArray::size(); + const unsigned int n_vect_bits = 8 * sizeof(number) * n_vect_doubles; + + ConditionalOStreams::pout_base() << "number of processes: " << n_proc << "\n" + << std::flush; + ConditionalOStreams::pout_base() + << "vectorization over " << n_vect_doubles << " doubles = " << n_vect_bits + << " bits (" << dealii::Utilities::System::get_current_vectorization_level() << ')' + << "\n" + << std::flush; + + // Create the mesh + ConditionalOStreams::pout_base() << "creating triangulation...\n" << std::flush; + Timer::start_section("Generate mesh"); + triangulation_manager.generate_mesh(); + Timer::end_section("Generate mesh"); + + // Create the dof handlers. + ConditionalOStreams::pout_base() << "creating DoFHandlers...\n" << std::flush; + Timer::start_section("reinitialize DoFHandlers"); + dof_manager.init(triangulation_manager, SystemWide::fe_systems); + Timer::end_section("reinitialize DoFHandlers"); + + // Create the constraints + ConditionalOStreams::pout_base() << "creating constraints...\n" << std::flush; + Timer::start_section("Create constraints"); + for (const auto &solve_group : solve_groups) + { + // TODO: Loop over levels, pass in current time + constraint_manager.make_constraints(SystemWide::mapping, + dof_manager.get_dof_handlers( + solve_group.field_indices)); + } + Timer::end_section("Create constraints"); + + // Initialize the solvers + Timer::start_section("Initialize Solvers"); + solvers.reserve(solve_groups.size()); + for (const auto &solve_group : solve_groups) + { + std::shared_ptr> solver; + switch (solve_group.pde_type) + { + case PDEType::Explicit: + solver = std::make_shared>(solve_group, + solve_context); + break; + case PDEType::Linear: + solver = std::make_shared>(solve_group, + solve_context); + break; + case PDEType::Newton: + solver = std::make_shared>(solve_group, + solve_context); + break; + default: + AssertThrow(false, dealii::ExcMessage("Unknown solver type")); + } + solvers.push_back(std::move(solver)); + } + Timer::end_section("Initialize Solvers"); + + // TODO: InvM and element volume + + // Update the ghosts + Timer::start_section("Update ghosts"); + for (auto solver : solvers) + { + solver->update_ghosts(); + } + Timer::end_section("Update ghosts"); + + // Perform the initial grid refinement. For this one, we have to do a loop to sufficient + // coarsen cells to the minimum level + ConditionalOStreams::pout_base() << "initializing grid refiner..." << std::flush; + grid_refiner.init(SystemWide::fe_systems); + grid_refiner.add_refinement_marker(std::make_shared>( + user_inputs.get_nucleation_parameters(), + pf_tools->nuclei_list)); + dealii::types::global_dof_index old_dofs = dof_manager.get_total_dofs(); + dealii::types::global_dof_index new_dofs = 0; + for (unsigned int remesh_index = 0; + remesh_index < (user_inputs.get_spatial_discretization().get_max_refinement() - + user_inputs.get_spatial_discretization().get_min_refinement()); + remesh_index++) + { + // Perform grid refinement + ConditionalOStreams::pout_base() << "performing grid refinement...\n" << std::flush; + Timer::start_section("Grid refinement"); + grid_refiner.do_adaptive_refinement(); + Timer::end_section("Grid refinement"); + + // Reinitialize the solvers + for (auto solver : solvers) + { + solver->reinit(); + } + // Update the ghosts + Timer::start_section("Update ghosts"); + for (auto solver : solvers) + { + solver->update_ghosts(); + } + Timer::end_section("Update ghosts"); + + // Recalculate the total DoFs + new_dofs = dof_manager.get_total_dofs(); + + // Check for convergence + if (old_dofs == new_dofs) + { + break; + } + old_dofs = new_dofs; + } +} + +template +void +Problem::solve() +{ + // Print a warning if running in DEBUG mode + ConditionalOStreams::pout_verbose() + << "\n\n" + << "================================================\n" + " Warning: running in DEBUG mode \n" + << "================================================\n\n\n" + << std::flush; + ConditionalOStreams::pout_summary() + << "================================================\n" + " Initialization\n" + << "================================================\n" + << std::flush; + + Timer::start_section("Initialization"); + init_system(); + Timer::end_section("Initialization"); + + ConditionalOStreams::pout_base() << "\n"; + + ConditionalOStreams::pout_summary() + << "================================================\n" + " Solve\n" + << "================================================\n" + << std::flush; + + const UserInputParameters &user_inputs = *user_inputs_ptr; + const TemporalDiscretization &time_info = user_inputs.get_temporal_discretization(); + SimulationTimer sim_timer(time_info.get_timestep()); + // Main time-stepping loop + while (sim_timer.get_increment() < time_info.get_total_increments()) + { + // Solve a single increment + // Includes nucleation, refinement, constraints, solve, output, and update + solve_increment(sim_timer); + // Update time + sim_timer.increment(); + } + + // Print summary of nuclei seeded during the simulation + ConditionalOStreams::pout_summary() + << "================================================\n" + " Nuclei Seeded\n" + << "================================================\n" + << std::to_string(pf_tools->nuclei_list.size()) << " total nuclei seeded.\n"; + for (const Nucleus nucleus : pf_tools->nuclei_list) + { + ConditionalOStreams::pout_summary() << nucleus << "\n"; + } + ConditionalOStreams::pout_summary() << "\n" << std::flush; + + // Print timer summary + Timer::print_summary(); +} + +template +void +Problem::solve_increment(SimulationTimer &sim_timer) +{ + const UserInputParameters &user_inputs = *user_inputs_ptr; + unsigned int increment = sim_timer.get_increment(); + // Check for stochastic nucleation + Timer::start_section("Check for nucleation"); + bool any_nucleation_occurred = false; + if (user_inputs.get_nucleation_parameters().should_attempt_nucleation(increment)) + { + any_nucleation_occurred = + NucleationHandler::attempt_nucleation(solve_context, + pf_tools->nuclei_list); + } + Timer::end_section("Check for nucleation"); + + // Perform grid refinement if necessary + if (user_inputs.get_spatial_discretization().get_has_adaptivity() && + (user_inputs.get_spatial_discretization().should_refine_mesh(increment) || + any_nucleation_occurred)) + { + // Perform grid refinement + ConditionalOStreams::pout_base() + << "[Increment " << sim_timer.get_increment() << "] : Grid Refinement\n"; + Timer::start_section("Grid refinement"); + grid_refiner.do_adaptive_refinement(); + Timer::end_section("Grid refinement"); + ConditionalOStreams::pout_base() << "\n" << std::flush; + + // Reinitialize the solvers + for (auto solver : solvers) + { + solver->reinit(); + } + + // Update the ghosts + Timer::start_section("Update ghosts"); + for (auto solver : solvers) + { + solver->update_ghosts(); + } + Timer::end_section("Update ghosts"); + } + + // Update the time-dependent constraints + Timer::start_section("Update time-dependent constraints"); + for (const auto &solve_group : solve_groups) + { + // TODO: Loop over levels, pass in current time + constraint_manager.update_time_dependent_constraints( + SystemWide::mapping, + dof_manager.get_dof_handlers(solve_group.field_indices)); + } + Timer::end_section("Update time-dependent constraints"); + + // Solve a single increment + Timer::start_section("Solve Increment"); + for (auto solver : solvers) + { + solver->solve(); + } + Timer::end_section("Solve Increment"); + + // Output results if needed + if (user_inputs.get_output_parameters().should_output(increment)) + { + Timer::start_section("Output"); + SolutionOutput(field_attributes, + solve_context->solution_indexer, + dof_manager, + degree, + "solution", + user_inputs); + + // Print the l2-norms and integrals of each solution + ConditionalOStreams::pout_base() + << "Iteration: " << sim_timer.get_increment() << "\n"; + for (unsigned int index = 0; index < field_attributes.size(); ++index) + { + const auto &solution = + solve_context->solution_indexer.get_solution_vector(index); + ConditionalOStreams::pout_base() + << " Solution index " << index << " l2-norm: " << solution.l2_norm() + << " integrated value: "; + + const auto tensor_rank = field_attributes[index].field_type; + + if (tensor_rank == FieldInfo::TensorRank::Vector) + { + ConditionalOStreams::pout_base() + << Integrator::template integrate< + FieldInfo::TensorRank::Vector>(dof_manager.get_dof_handler(index), + solution) + << "\n"; + } + else if (tensor_rank == FieldInfo::TensorRank::Scalar) + { + ConditionalOStreams::pout_base() + << Integrator::template integrate< + FieldInfo::TensorRank::Scalar>(dof_manager.get_dof_handler(index), + solution) + << "\n"; + } + } + ConditionalOStreams::pout_base() << "\n" << std::flush; + Timer::end_section("Output"); + } + + // Update the field labels in preparation for next increment (c_n -> c_n-1) + for (auto solver : solvers) + { + solver.update(); + } +} + +// #include "core/problem.inst" + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/src/core/solution_indexer.cc b/src/core/solution_indexer.cc new file mode 100644 index 000000000..77dfb25d7 --- /dev/null +++ b/src/core/solution_indexer.cc @@ -0,0 +1,128 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +#include +#include + +#include + + +PRISMS_PF_BEGIN_NAMESPACE + +template +SolutionIndexer::SolutionIndexer( + unsigned int num_fields, // attributes_list.size() + const std::set> &solution_handlers) +{ + solutions.clear(); + solutions.resize(num_fields, nullptr); + for (const GroupSolutionHandler &solution_handler : solution_handlers) + { + for (unsigned int field_index : solution_handler.get_solve_group().field_indices) + { + solutions[field_index] = &solution_handler; + } + } +} + +template +auto +SolutionIndexer::get_solution_vector(unsigned int global_index, + unsigned int relative_level) const + -> const SolutionVector & +{ + return solutions[global_index]->get_solution_vector(global_index, relative_level); +} + +template +auto +SolutionIndexer::get_solution_vector(unsigned int global_index, + unsigned int relative_level) + -> SolutionVector & +{ + return solutions[global_index]->get_solution_vector(global_index, relative_level); +} + +template +auto +SolutionIndexer::get_old_solution_vector(unsigned int age, + unsigned int global_index, + unsigned int relative_level) const + -> const SolutionVector & +{ + return solutions[global_index]->get_old_solution_vector(age, + global_index, + relative_level); +} + +template +auto +SolutionIndexer::get_old_solution_vector(unsigned int age, + unsigned int global_index, + unsigned int relative_level) + -> SolutionVector & +{ + return solutions[global_index]->get_old_solution_vector(age, + global_index, + relative_level); +} + +template +auto +SolutionIndexer::get_new_solution_vector(unsigned int global_index, + unsigned int relative_level) const + -> const SolutionVector & +{ + return solutions[global_index]->get_new_solution_vector(global_index, relative_level); +} + +template +auto +SolutionIndexer::get_new_solution_vector(unsigned int global_index, + unsigned int relative_level) + -> SolutionVector & +{ + return solutions[global_index]->get_new_solution_vector(global_index, relative_level); +} + +template +auto +SolutionIndexer::get_matrix_free(unsigned int global_index, + unsigned int relative_level) const + -> const MatrixFree & +{ + return solutions[global_index]->get_matrix_free(relative_level); +} + +template +auto +SolutionIndexer::get_matrix_free(unsigned int global_index, + unsigned int relative_level) -> MatrixFree & +{ + return solutions[global_index]->get_matrix_free(relative_level); +} + +template +auto +SolutionIndexer::get_solution_level_and_block_index( + unsigned int global_index, + unsigned int relative_level) + -> std::pair &, unsigned int> +{ + auto &sol = solutions[global_index]; + return {sol->get_solution_level(relative_level), sol->get_block_index(global_index)}; +} + +template +unsigned int +SolutionIndexer::get_block_index(unsigned int global_index) const + +{ + return solutions[global_index]->get_block_index(global_index); +} + +// #include "core/solution_indexer.inst" + +PRISMS_PF_END_NAMESPACE diff --git a/src/core/solution_output.cc b/src/core/solution_output.cc index 2b33c2fdd..47bbd44c5 100644 --- a/src/core/solution_output.cc +++ b/src/core/solution_output.cc @@ -14,6 +14,10 @@ #include +#include "prismspf/core/dof_manager.h" +#include "prismspf/core/field_attributes.h" +#include "prismspf/core/solution_indexer.h" + #include #include #include @@ -216,6 +220,115 @@ SolutionOutput::SolutionOutput( } } +template +SolutionOutput::SolutionOutput( + const std::vector &field_attributes, + const SolutionIndexer &solution_indexer, + const DofManager &dof_manager, + const unsigned int °ree, + const std::string &file_prefix, + const UserInputParameters &user_inputs) +{ + // Some stuff to determine the actual name of the output file. + const auto n_trailing_digits = static_cast( + std::floor( + std::log10(user_inputs.get_temporal_discretization().get_total_increments())) + + 1); + + // Init data out + dealii::DataOut data_out; + + // Add data vectors + for (unsigned int index = 0; index < field_attributes.size(); ++index) + { + const FieldAttributes &field = field_attributes[index]; + auto &solution = solution_indexer.get_solution_vector(index); + solution.update_ghost_values(); + + // Mark field as Scalar/Vector + const unsigned int n_components = + (field.field_type == FieldInfo::TensorRank::Scalar) ? 1 : dim; + + const std::vector + data_type(n_components, + (n_components == 1) + ? dealii::DataComponentInterpretation::component_is_scalar + : dealii::DataComponentInterpretation::component_is_part_of_vector); + + const std::vector names(n_components, field.name); + + data_out.add_data_vector(dof_manager.get_dof_handler(index), + solution, + names, + data_type); + + solution.zero_out_ghost_values(); + } + + // Build patches to linearly interpolate from higher order element degrees. Note that + // this essentially converts the element to an equal amount of subdivisions in the + // output. This does not make subdivisions and element degree equivalent in the + // simulation! + const unsigned int n_divisions = + user_inputs.get_output_parameters().get_patch_subdivisions() == 0 + ? degree + : user_inputs.get_output_parameters().get_patch_subdivisions(); + data_out.build_patches(n_divisions); + + // Set some flags for data output + dealii::DataOutBase::VtkFlags flags; + flags.time = user_inputs.get_temporal_discretization().get_time(); + flags.cycle = user_inputs.get_temporal_discretization().get_increment(); + flags.print_date_and_time = true; +#ifdef PRISMS_PF_WITH_ZLIB + flags.compression_level = dealii::DataOutBase::CompressionLevel::best_speed; +#endif + data_out.set_flags(flags); + + // Write to file based on the user input. + const std::string directory = "./"; + const unsigned int increment = + user_inputs.get_temporal_discretization().get_increment(); + + if (user_inputs.get_output_parameters().get_file_type() == "vtu") + { + std::ostringstream increment_stream; + increment_stream << std::setw(static_cast(n_trailing_digits)) + << std::setfill('0') << increment; + const std::string filename = + directory + file_prefix + "_" + increment_stream.str() + ".vtu"; + data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD); + } + else if (user_inputs.get_output_parameters().get_file_type() == "pvtu") + { + data_out.write_vtu_with_pvtu_record(directory, + file_prefix, + increment, + MPI_COMM_WORLD, + n_trailing_digits); + } + else if (user_inputs.get_output_parameters().get_file_type() == "vtk") + { + std::ostringstream increment_stream; + increment_stream << std::setw(static_cast(n_trailing_digits)) + << std::setfill('0') << increment; + const std::string filename = + directory + file_prefix + "_" + increment_stream.str() + ".vtk"; + std::ofstream vtk_output(filename); + data_out.write_vtk(vtk_output); + } + else + { + AssertThrow(false, UnreachableCode()); + } + + // Update the ghost values again to allow for read access + for (unsigned int index = 0; index < field_attributes.size(); ++index) + { + solution_indexer.get_solution_vector(index).update_ghost_values(); + } +} + #include "core/solution_output.inst" PRISMS_PF_END_NAMESPACE diff --git a/src/core/triangulation_manager.cc b/src/core/triangulation_manager.cc new file mode 100644 index 000000000..923565bf9 --- /dev/null +++ b/src/core/triangulation_manager.cc @@ -0,0 +1,231 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template <> +TriangulationManager<1U>::TriangulationManager(bool _has_multigrid) + : has_multigrid(_has_multigrid) + , triangulation(dealii::Triangulation<1U>::limit_level_difference_at_vertices) +{} + +template +TriangulationManager::TriangulationManager(bool _has_multigrid) + : has_multigrid(_has_multigrid) + , triangulation(MPI_COMM_WORLD, + dealii::Triangulation::limit_level_difference_at_vertices) +{} + +template +const typename TriangulationManager::Triangulation & +TriangulationManager::get_triangulation() const +{ + return triangulation; +} + +template +const std::vector>> & +TriangulationManager::get_mg_triangulation() const +{ + Assert(!coarsened_triangulations.empty(), dealii::ExcNotInitialized()); + return coarsened_triangulations; +} + +template +const dealii::Triangulation & +TriangulationManager::get_triangulation(unsigned int relative_level) const +{ + Assert(!coarsened_triangulations.empty(), dealii::ExcNotInitialized()); + Assert(coarsened_triangulations.size() >= relative_level, + dealii::ExcMessage( + "The coarse triangulation set does not contain that specified level")); + return *coarsened_triangulations[relative_level]; +} + +template +void +TriangulationManager::generate_mesh(const UserInputParameters &user_inputs) +{ + const SpatialDiscretization &discretization_params = + user_inputs.get_spatial_discretization(); + // TODO (landinjm): Add more generality in selecting mesh types + if (discretization_params.get_radius() != 0.0) + { + // TODO (landinjm): Adding assertion about periodic boundary conditions for spheres + // Generate a sphere + + // TODO (landinjm): Add assertion that the user cannot specify multiple boundary + // conditions for a hyper_ball geometry + dealii::GridGenerator::hyper_ball(triangulation, + dealii::Point(), + discretization_params.get_radius()); + } + else + { + // TODO (landinjm): Add assertions about periodic boundary conditions for + // Rectangular domains here. Not sure whether it is better to check for assertions + // here or when we parse user inputs. + + // Generate rectangle + dealii::GridGenerator::subdivided_hyper_rectangle( + triangulation, + discretization_params.get_subdivisions(), + discretization_params.get_lower_bound(), + discretization_params.get_upper_bound()); + + // Mark boundaries. This is done before global refinement to reduce the number of + // cells we have to loop through. + mark_boundaries(discretization_params); + + // Mark periodicity + mark_periodic(user_inputs); + } + + // TODO (landinjm): This be better as a combination of a parameter flag and debug + // mode. Output triangulation to vtk if in debug mode +#ifdef DEBUG +// if(user_inputs ... output triangulation) +#endif + { + export_triangulation_as_vtk("triangulation"); + } + // Global refinement + triangulation.refine_global(discretization_params.get_global_refinement()); + + // Create the triangulations for the coarser levels if we have multigrid for any of the + // fields + reinit(); +} + +template +void +TriangulationManager::export_triangulation_as_vtk(const std::string &filename) const +{ + const dealii::GridOut grid_out; + std::ofstream out(filename + ".vtk"); + grid_out.write_vtk(triangulation, out); + ConditionalOStreams::pout_base() << "Triangulation written to " << filename << ".vtk\n"; +} + +// TODO (fractalsbyx): See if this can be made more efficient +template +void +TriangulationManager::mark_boundaries( + const SpatialDiscretization &discretization_params) const +{ + const double tolerance = 1e-12; + + // Loop through the cells + for (const auto &cell : triangulation.active_cell_iterators()) + { + // Mark the faces (faces_per_cell = 2*dim) + for (unsigned int face_number = 0; + face_number < dealii::GeometryInfo::faces_per_cell; + ++face_number) + { + // Direction for quad and hex cells + auto direction = face_number / 2; + + // Lower bound and upper bound + const double lower_bound = discretization_params.get_lower_bound()[direction]; + const double upper_bound = discretization_params.get_upper_bound()[direction]; + + // Mark the boundary id for lower and upper bounds + if (std::fabs(cell->face(face_number)->center()(direction) - lower_bound) < + tolerance || + std::fabs(cell->face(face_number)->center()(direction) - upper_bound) < + tolerance) + { + cell->face(face_number)->set_boundary_id(face_number); + } + } + } +} + +template +void +TriangulationManager::mark_periodic(const UserInputParameters &user_inputs) +{ + // Create a little set of boundary ids we've already marked + std::set periodic_ids; + + // Add periodicity in the triangulation where specified in the boundary conditions. Note + // that if one field is periodic all others should be as well. + // TODO (fractalsbyx): Enforce above condition systematically + for (const auto &[index, boundary_condition] : + user_inputs.get_boundary_parameters().get_boundary_condition_list()) + { + for (const auto &[component, condition] : boundary_condition) + { + for (const auto &[boundary_id, boundary_type] : + condition.get_boundary_condition_map()) + { + if (boundary_type == BoundaryCondition::Type::Periodic) + { + // Skip boundary ids that are odd since those map to the even faces + if (boundary_id % 2 != 0) + { + continue; + } + + // Skip the id if we've already added periodic boundaries to id + if (!periodic_ids.insert(boundary_id).second) + { + continue; + } + + // Create a vector of matched pairs that we fill and enforce upon the + // constaints + std::vector> + periodicity_vector; + + // Determine the direction + const auto direction = boundary_id / 2; + + // Grab the offset vector from one vertex to another + dealii::Tensor<1, dim> offset; + offset[direction] = + user_inputs.get_spatial_discretization().get_size()[direction]; + + // Collect the matched pairs on the coarsest level of the mesh + dealii::GridTools::collect_periodic_faces(triangulation, + boundary_id, + boundary_id + 1, + direction, + periodicity_vector, + offset); + + // Set constraints + triangulation.add_periodicity(periodicity_vector); + } + } + } + } +} + +// #include "core/triangulation_manager.inst" + +PRISMS_PF_END_NAMESPACE diff --git a/src/solvers/group_explicit_solver.cc b/src/solvers/group_explicit_solver.cc new file mode 100644 index 000000000..e152a07c2 --- /dev/null +++ b/src/solvers/group_explicit_solver.cc @@ -0,0 +1,10 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +// #include "solvers/group_explicit_solver.inst" + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/src/solvers/group_linear_solver.cc b/src/solvers/group_linear_solver.cc new file mode 100644 index 000000000..eb348e7bf --- /dev/null +++ b/src/solvers/group_linear_solver.cc @@ -0,0 +1,10 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +// #include "solvers/group_explicit_solver.inst" + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/src/solvers/group_newton_solver.cc b/src/solvers/group_newton_solver.cc new file mode 100644 index 000000000..cd1e96c4a --- /dev/null +++ b/src/solvers/group_newton_solver.cc @@ -0,0 +1,10 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +// #include "solvers/group_explicit_solver.inst" + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/src/solvers/group_solver_base.cc b/src/solvers/group_solver_base.cc new file mode 100644 index 000000000..f2c76a049 --- /dev/null +++ b/src/solvers/group_solver_base.cc @@ -0,0 +1,10 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +// #include "solvers/group_solver_base.inst" + +PRISMS_PF_END_NAMESPACE \ No newline at end of file diff --git a/src/solvers/mf_operator.cc b/src/solvers/mf_operator.cc new file mode 100644 index 000000000..c862ac929 --- /dev/null +++ b/src/solvers/mf_operator.cc @@ -0,0 +1,344 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include + +#include + +#include "prismspf/core/group_solution_handler.h" +#include "prismspf/core/solve_group.h" + +#include +#include + +#if DEAL_II_VERSION_MAJOR >= 9 && DEAL_II_VERSION_MINOR >= 7 +# include +#else +# include +#endif + +PRISMS_PF_BEGIN_NAMESPACE + +template +void +MFOperator::compute_operator(BlockVector &dst, + const BlockVector &src) const +{ + data->cell_loop(&MFOperator::compute_local_operator, this, dst, src, true); +} + +template +void +MFOperator::compute_local_operator( + const MatrixFree &_data, + BlockVector &dst, + const BlockVector &src, + const std::pair &cell_range) const +{ + // Construct FEEvaluation objects + // The reason this is constructed here, rather than as a private member is because + // compute_local_rhs is called by cell_loop, which multithreads. There would be data + // races. + FieldContainer variable_list(1 /*args*/); + DSTContainer dst_fields(solve_group.field_indices, + field_attributes, + *data, + field_to_block_index); + + // Initialize, evaluate, and submit based on user function. + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + // Grab the element volume + // const ScalarValue element_volume = + // element_volume_handler->get_element_volume(cell); + + // Initialize, read DOFs, and set evaulation flags for each variable + variable_list.reinit_and_eval(cell); + dst_fields.reinit(cell); + + // Evaluate the user-defined pde at each quadrature point + for (unsigned int quad = 0; quad < variable_list.get_n_q_points(); ++quad) + { + variable_list.set_q_point(quad); + dst_fields.set_q_point(quad); + // Evaluate the function pointer (the user-defined pde) + pde_operator->*pde_op(variable_list, dst_fields); + } + + // Integrate and add to global vector dst + for (unsigned int field_index : solve_group.field_indices) + { + dst_fields.integrate_and_distribute(field_index, + dst.block( + field_to_block_index[field_index])); + } + } +} + +template +void +MFOperator::compute_diagonal() +{ + inverse_diagonal_entries.reset(new dealii::DiagonalMatrix()); + SolutionVector &inverse_diagonal = inverse_diagonal_entries->get_vector(); + data->initialize_dof_vector(inverse_diagonal, field_index); + const unsigned int dummy = 0; + data->cell_loop(&MFOperator::local_compute_diagonal, this, inverse_diagonal, dummy); + + set_constrained_entries_to_one(inverse_diagonal); + + for (unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i) + { + Assert(inverse_diagonal.local_element(i) > static_cast(0.0), + dealii::ExcMessage( + "No diagonal entry in a positive definite operator should be zero")); + inverse_diagonal.local_element(i) = number(1.0) / inverse_diagonal.local_element(i); + } +} + +template +void +MFOperator::compute_local_diagonal( + const dealii::MatrixFree> &_data, + BlockVector &diagonal, + [[maybe_unused]] const unsigned int &dummy, + const std::pair &cell_range) const +{ + // Construct FEEvaluation objects + // The reason this is constructed here, rather than as a private member is because + // compute_local_rhs is called by cell_loop, which multithreads. There would be data + // races. + FieldContainer variable_list(1 /*args*/); + DSTContainer dst_fields(solve_group.field_indices, + field_attributes, + *data, + field_to_block_index); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + // Reinit the cell for all the dependencies + variable_list.reinit_and_eval(cell); + dst_fields.reinit(cell); + + // To get the diagonal of the "matrix", repeatedly "multiply the matrix" by a test x + // vector (change vector) and store the solution in the i-th position in the y + // vector (the diagonal). + // First (i=1) x vector: I 0 0 0 + // Second(i=2) x vector: 0 I 0 0 ... + // Submit zeros for everyting except the diagonals + for (unsigned int field_index : solve_group.field_indices) + { + if (/* Scalar */) + { + dealii::AlignedVector cell_diagonal = + compute_field_diagonal(variable_list, + dst_fields, + field_index); + // Submit calculated diagonal values and distribute + for (unsigned int i = 0; i < n_dofs_per_cell; ++i) + { + feeval_ptr->submit_dof_value(cell_diagonal[i], i); + } + feeval_ptr->distribute_local_to_global(dst); + } + + // Submit calculated diagonal values and distribute + for (unsigned int i = 0; i < n_dofs_per_cell; ++i) + { + if (/* Scalar */) + { + feeval_ptr->submit_dof_value((*diagonal_ptr)[i], i); + } + else + { + feeval_ptr->submit_dof_value((*diagonal_ptr)[i][0], i); + } + } + feeval_ptr->distribute_local_to_global( + diagonal.block(field_to_block_index[field_index])); + } + } +} + +template +template ::TensorRank Rank> +auto +MFOperator::compute_field_diagonal( + FieldContainer &variable_list, + DSTContainer &dst_fields, + unsigned int field_index) const -> dealii::AlignedVector> +{ + unsigned int n_dofs_per_cell = variable_list.get_dofs_per_component(field_index); + dealii::AlignedVector> cell_diagonal(n_dofs_per_cell, zero()); + // vector_feeval_ptr->dofs_per_component; + // scalar_feeval_ptr->dofs_per_cell; + for (unsigned int i = 0; i < n_dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < n_dofs_per_cell; ++j) + { + dst_fields.set_dof_value(field_index, + i == j ? identity() : zero(), + j); + } + + // Evaluate the dependencies based on the flags + variable_list.eval(); + + // Evaluate at each quadrature point + for (unsigned int quad = 0; quad < variable_list.get_n_q_points(); ++quad) + { + variable_list.set_q_point(quad); + dst_fields.set_q_point(quad); + pde_operator->pde_op(variable_list, dst_fields); + } + + // Integrate the diagonal + dst_fields.integrate(field_index); + + // TODO: fix this + dst_fields.eval(); + cell_diagonal[i] = + dst_fields.get_dof_value(field_index, + i, + variable_list.get_dof_value(field_index, i)); + } + return cell_diagonal; +} + +template +void +MFOperator::compute_element_volume(SolutionVector &dst) +{ + int dummy = 0; + data->cell_loop(&MFOperator::compute_local_element_volume, this, dst, dummy); +} + +template +void +MFOperator::compute_local_element_volume( + const dealii::MatrixFree> &_data, + SolutionVector &dst, + [[maybe_unused]] const int &dummy, + const std::pair &cell_range) const +{ + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + } +} + +template +dealii::types::global_dof_index +MFOperator::m() const +{ + Assert(data.get() != nullptr, dealii::ExcNotInitialized()); + + const unsigned int total_size = + std::accumulate(selected_fields.begin(), + selected_fields.end(), + 0U, + [this](unsigned int sum, unsigned int field) + { + return sum + data->get_vector_partitioner(field)->size(); + }); + + return total_size; +} + +template +number +MFOperator::el([[maybe_unused]] const unsigned int &row, + [[maybe_unused]] const unsigned int &col) const +{ + AssertThrow(false, FeatureNotImplemented("el()")); + return 0.0; +} + +template +void +MFOperator::clear() +{ + data.reset(); + diagonal_entries.reset(); + inverse_diagonal_entries.reset(); +} + +template +void +MFOperator::initialize_dof_vector( + SolutionVector &dst, + unsigned int dof_handler_index) const +{ + data->initialize_dof_vector(dst, dof_handler_index); +} + +template +void +MFOperator::set_constrained_entries_to_one(SolutionVector &dst) const +{ + for (unsigned int j = 0; j < dealii::MFOperators::BlockHelper::n_blocks(dst); ++j) + { + const std::vector &constrained_dofs = + data->get_constrained_dofs(selected_fields[j]); + for (const auto constrained_dof : constrained_dofs) + { + dealii::MFOperators::BlockHelper::subblock(dst, j).local_element( + constrained_dof) = 1.0; + } + } +} + +template +std::shared_ptr>> +MFOperator::get_matrix_free() const +{ + return data; +} + +template +const std::shared_ptr< + dealii::DiagonalMatrix::SolutionVector>> & +MFOperator::get_matrix_diagonal_inverse() const +{ + Assert(inverse_diagonal_entries.get() != nullptr && inverse_diagonal_entries->m() > 0, + dealii::ExcNotInitialized()); + return inverse_diagonal_entries; +} + +template +void +MFOperator::vmult(SolutionVector &dst, + const SolutionVector &src) const +{ + compute_operator(dst, src); +} + +// NOLINTBEGIN(readability-identifier-naming) + +template +void +MFOperator::Tvmult(SolutionVector &dst, + const SolutionVector &src) const +{ + vmult(dst, src); +} + +// NOLINTEND(readability-identifier-naming) + +// #include "core/matrix_free_operator.inst" + +PRISMS_PF_END_NAMESPACE