From 8210bafaac00d105eca7550562c68b491b1bc2da Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Mon, 23 Sep 2024 14:19:53 +0200 Subject: [PATCH] Simplify solution evaluator --- include/aspect/solution_evaluator.h | 169 ++++++++++++++++++++++++++-- source/solution_evaluator.cc | 167 ++++----------------------- 2 files changed, 176 insertions(+), 160 deletions(-) diff --git a/include/aspect/solution_evaluator.h b/include/aspect/solution_evaluator.h index 49eb28956f9..cd5e79fa851 100644 --- a/include/aspect/solution_evaluator.h +++ b/include/aspect/solution_evaluator.h @@ -35,6 +35,103 @@ namespace aspect { using namespace dealii; + namespace internal + { + /** + * Wrapper around dealii::FEPointEvaluation to choose number of components dynamically. + * This class is only an abstract base class that can be used to implement derived classes. + */ + template + class DynamicFEPointEvaluation + { + public: + /** + * Constructor which allows to select at run time the @p first_component and the @p n_components. + */ + DynamicFEPointEvaluation(const unsigned int first_component, const unsigned int n_components) + : first_component (first_component), + n_components (n_components) + {} + + /** + * Destructor. + */ + virtual ~DynamicFEPointEvaluation() = default; + + /** + * Evaluate the solution at the given positions. + * + * @p solution_values contains the values of the degrees of freedom. + * @p flags controls which values should be computed. + */ + virtual void evaluate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags flags) = 0; + + /** + * Return the gradient of the solution at the given @p evaluation_point. + */ + virtual + small_vector> + get_gradient(const unsigned int evaluation_point) const = 0; + + /** + * Copy the gradient of the solution at the given @p evaluation_point + * into the array @p gradients. + */ + virtual + void + get_gradient(const unsigned int evaluation_point, + const ArrayView> &gradients) const = 0; + + /** + * Return the value of the solution at the given @p evaluation_point. + */ + virtual + small_vector + get_value(const unsigned int evaluation_point) const = 0; + + /** + * Copy the value of the solution at the given @p evaluation_point. + * into the array @p solution. + */ + virtual + void + get_value(const unsigned int evaluation_point, + const ArrayView &solution) const = 0; + + /** + * Return the first component of the solution vector. + */ + virtual + unsigned int + get_first_component() const final + { + return first_component; + } + + /** + * Return the number of components of the solution vector. + */ + virtual + unsigned int + get_n_components() const final + { + return n_components; + } + + private: + /** + * The first component of the solution vector. + */ + unsigned int first_component; + + /** + * The number of components of the solution vector. + */ + unsigned int n_components; + }; + } + /** * This class evaluates the solution vector at arbitrary positions inside a cell. * This base class only provides the interface for SolutionEvaluatorImplementation. @@ -45,9 +142,12 @@ namespace aspect { public: /** - * virtual Destructor. + * Constructor. Create the member variables given a simulator access @p simulator and a set of + * update flags @p update_flags. The @p update_flags control if only the solution or also the gradients + * should be evaluated. */ - virtual ~SolutionEvaluator() = default; + SolutionEvaluator(const SimulatorAccess &simulator, + const UpdateFlags update_flags); /** * Reinitialize all variables to evaluate the given solution for the given cell @@ -55,11 +155,11 @@ namespace aspect * also the gradients should be evaluated. * If other flags are set an assertion is triggered. */ - virtual void + void reinit(const typename DoFHandler::active_cell_iterator &cell, const ArrayView> &positions, const ArrayView &solution_values, - const UpdateFlags update_flags) = 0; + const UpdateFlags update_flags); /** * Fill @p solution with all solution components at the given @p evaluation_point. Note @@ -67,8 +167,8 @@ namespace aspect * because this function only returns the results of the computation that * happened in reinit(). */ - virtual void get_solution(const unsigned int evaluation_point, - const ArrayView &solution) = 0; + void get_solution(const unsigned int evaluation_point, + const ArrayView &solution); /** * Fill @p gradients with all solution gradients at the given @p evaluation_point. Note @@ -76,21 +176,66 @@ namespace aspect * because this function only returns the results of the computation that * happened in reinit(). */ - virtual void get_gradients(const unsigned int evaluation_point, - const ArrayView> &gradients) = 0; + void get_gradients(const unsigned int evaluation_point, + const ArrayView> &gradients); /** * Return the evaluator for velocity or fluid velocity. This is the only * information necessary for advecting particles. */ - virtual FEPointEvaluation & - get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity) = 0; + FEPointEvaluation & + get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity); /** * Return the cached mapping information. */ - virtual NonMatching::MappingInfo & - get_mapping_info() = 0; + NonMatching::MappingInfo & + get_mapping_info(); + + private: + /** + * MappingInfo object for the FEPointEvaluation objects + */ + NonMatching::MappingInfo mapping_info; + + /** + * FEPointEvaluation objects for all common + * components of ASPECT's finite element solution. + * These objects are used inside of the member functions of this class. + */ + FEPointEvaluation velocity; + std::unique_ptr> pressure; + FEPointEvaluation<1, dim> temperature; + + /** + * We group compositions by type (FiniteElement) and evaluate + * them together. + */ + std::vector>> compositions; + + /** + * Pointers to FEPointEvaluation objects for all melt + * components of ASPECT's finite element solution, which only + * point to valid objects in case we use melt transport. Other + * documentation like for the objects directly above. + */ + std::unique_ptr> fluid_velocity; + std::unique_ptr> compaction_pressure; + std::unique_ptr> fluid_pressure; + + /** + * The component indices for the three melt formulation + * variables fluid velocity, compaction pressure, and + * fluid pressure (in this order). They are cached + * to avoid repeated expensive lookups. + */ + std::array melt_component_indices; + + /** + * Reference to the active simulator access object. Provides + * access to the general simulation variables. + */ + const SimulatorAccess &simulator_access; }; /** diff --git a/source/solution_evaluator.cc b/source/solution_evaluator.cc index 85fe1b04dec..a63de8be6e1 100644 --- a/source/solution_evaluator.cc +++ b/source/solution_evaluator.cc @@ -28,48 +28,6 @@ namespace aspect { namespace internal { - - /** - * Wrapper around dealii::FEPointEvaluation to choose number of components dynamically. - */ - template - class DynamicFEPointEvaluation - { - public: - DynamicFEPointEvaluation(const unsigned int first_component, const unsigned int n_components) - : first_component (first_component), - n_components (n_components) - {} - - virtual ~DynamicFEPointEvaluation() = default; - - unsigned int first_component; - unsigned int n_components; - - virtual void evaluate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags flags) = 0; - - virtual - small_vector> - get_gradient(const unsigned int evaluation_point) const = 0; - - virtual - void - get_gradient(const unsigned int evaluation_point, - const ArrayView> &gradients) const = 0; - - virtual - small_vector - get_value(const unsigned int evaluation_point) const = 0; - - virtual - void - get_value(const unsigned int evaluation_point, - const ArrayView &solution) const = 0; - }; - - - /** * The functions in this namespace allow us to use scalar and vector-valued FEEvaluation objects in the same * way, as the return type of FEEvaluation is double for one component, but a Tensor for more than one @@ -116,7 +74,7 @@ namespace aspect /** * Implementation of the base class DynamicFEPointEvaluation that wraps - * an FEEvaluation object. + * an FEPointEvaluation object. */ template class DynamicFEPointEvaluationImpl: public DynamicFEPointEvaluation @@ -177,6 +135,7 @@ namespace aspect gradients[c] = x[c]; } + private: FEPointEvaluation evaluation; }; @@ -221,98 +180,10 @@ namespace aspect } - // This class evaluates the solution vector at arbitrary positions inside a cell. - // It uses the deal.II class FEPointEvaluation to do this efficiently. Because - // FEPointEvaluation only supports a single finite element, but ASPECT uses a FESystem with - // many components, this class creates several FEPointEvaluation objects that are used for - // the individual finite elements of our solution (pressure, velocity, temperature, and - // all other optional variables). Because FEPointEvaluation is templated based on the - // number of components, but ASPECT only knows the number of components at runtime - // we create this derived class with an additional template. This makes it possible - // to access the functionality through the base class, but create an object of this - // derived class with the correct number of components at runtime. - template - class SolutionEvaluatorImplementation: public SolutionEvaluator - { - public: - // Constructor. Create the member variables given a simulator and a set of - // update flags. The update flags control if only the solution or also the gradients - // should be evaluated. - SolutionEvaluatorImplementation(const SimulatorAccess &simulator, - const UpdateFlags update_flags); - - // Reinitialize all variables to evaluate the given solution for the given cell - // and the given positions. The update flags control if only the solution or - // also the gradients should be evaluated. - // If other flags are set an assertion is triggered. - void - reinit(const typename DoFHandler::active_cell_iterator &cell, - const ArrayView> &positions, - const ArrayView &solution_values, - const UpdateFlags update_flags) override; - - // Return the value of all solution components at the given evaluation point. Note - // that this function only works after a successful call to reinit(), - // because this function only returns the results of the computation that - // happened in reinit(). - void get_solution(const unsigned int evaluation_point, - const ArrayView &solution) override; - - // Return the value of all solution gradients at the given evaluation point. Note - // that this function only works after a successful call to reinit(), - // because this function only returns the results of the computation that - // happened in reinit(). - void get_gradients(const unsigned int evaluation_point, - const ArrayView> &gradients) override; - - // Return the evaluator for velocity or fluid velocity. This is the only - // information necessary for advecting particles. - FEPointEvaluation & - get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity) override; - - // Return the cached mapping information. - NonMatching::MappingInfo & - get_mapping_info() override; - - private: - // MappingInfo object for the FEPointEvaluation objects - NonMatching::MappingInfo mapping_info; - - // FEPointEvaluation objects for all common - // components of ASPECT's finite element solution. - // These objects are used inside of the member functions of this class. - FEPointEvaluation velocity; - std::unique_ptr> pressure; - FEPointEvaluation<1, dim> temperature; - - // We group compositions by type (FiniteElement) and evaluate - // them together. - std::vector>> compositions; - - // Pointers to FEPointEvaluation objects for all melt - // components of ASPECT's finite element solution, which only - // point to valid objects in case we use melt transport. Other - // documentation like for the objects directly above. - std::unique_ptr> fluid_velocity; - std::unique_ptr> compaction_pressure; - std::unique_ptr> fluid_pressure; - - // The component indices for the three melt formulation - // variables fluid velocity, compaction pressure, and - // fluid pressure (in this order). They are cached - // to avoid repeated expensive lookups. - std::array melt_component_indices; - - // Reference to the active simulator access object. Provides - // access to the general simulation variables. - const SimulatorAccess &simulator_access; - }; - - template - SolutionEvaluatorImplementation::SolutionEvaluatorImplementation(const SimulatorAccess &simulator, - const UpdateFlags update_flags) + SolutionEvaluator::SolutionEvaluator(const SimulatorAccess &simulator, + const UpdateFlags update_flags) : mapping_info(simulator.get_mapping(), update_flags), @@ -409,10 +280,10 @@ namespace aspect template void - SolutionEvaluatorImplementation::reinit(const typename DoFHandler::active_cell_iterator &cell, - const ArrayView> &positions, - const ArrayView &solution_values, - const UpdateFlags update_flags) + SolutionEvaluator::reinit(const typename DoFHandler::active_cell_iterator &cell, + const ArrayView> &positions, + const ArrayView &solution_values, + const UpdateFlags update_flags) { // FEPointEvaluation uses different evaluation flags than the common UpdateFlags. // Translate between the two. @@ -511,8 +382,8 @@ namespace aspect template void - SolutionEvaluatorImplementation::get_solution(const unsigned int evaluation_point, - const ArrayView &solution) + SolutionEvaluator::get_solution(const unsigned int evaluation_point, + const ArrayView &solution) { Assert(solution.size() == simulator_access.introspection().n_components, ExcDimensionMismatch(solution.size(), simulator_access.introspection().n_components)); @@ -528,8 +399,8 @@ namespace aspect for (const auto &eval : compositions) { - const unsigned int start_index = eval->first_component; - const unsigned int n_components = eval->n_components; + const unsigned int start_index = eval->get_first_component(); + const unsigned int n_components = eval->get_n_components(); eval->get_value(evaluation_point, {&solution[start_index],n_components}); } @@ -549,8 +420,8 @@ namespace aspect template void - SolutionEvaluatorImplementation::get_gradients(const unsigned int evaluation_point, - const ArrayView> &gradients) + SolutionEvaluator::get_gradients(const unsigned int evaluation_point, + const ArrayView> &gradients) { Assert(gradients.size() == simulator_access.introspection().n_components, ExcDimensionMismatch(gradients.size(), simulator_access.introspection().n_components)); @@ -566,8 +437,8 @@ namespace aspect for (const auto &eval : compositions) { - const unsigned int start_index = eval->first_component; - const unsigned int n_components = eval->n_components; + const unsigned int start_index = eval->get_first_component(); + const unsigned int n_components = eval->get_n_components(); eval->get_gradient(evaluation_point, {&gradients[start_index],n_components}); @@ -587,7 +458,7 @@ namespace aspect template FEPointEvaluation & - SolutionEvaluatorImplementation::get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity) + SolutionEvaluator::get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity) { if (use_fluid_velocity) return *fluid_velocity; @@ -598,7 +469,7 @@ namespace aspect } template NonMatching::MappingInfo & - SolutionEvaluatorImplementation::get_mapping_info() + SolutionEvaluator::get_mapping_info() { return mapping_info; } @@ -610,7 +481,7 @@ namespace aspect construct_solution_evaluator (const SimulatorAccess &simulator_access, const UpdateFlags update_flags) { - return std::make_unique>(simulator_access, update_flags); + return std::make_unique>(simulator_access, update_flags); } }