From 4a31088f6ed64855d12f073158c965248e3dabed Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Tue, 24 Sep 2024 16:34:00 +0200 Subject: [PATCH 1/2] Rework particle update interface --- doc/modules/changes/20240925_gassmoeller | 8 + .../aspect/particle/property/composition.h | 9 +- .../particle/property/cpo_bingham_average.h | 31 +- .../particle/property/cpo_elastic_tensor.h | 31 +- .../property/crystal_preferred_orientation.h | 31 +- .../aspect/particle/property/elastic_stress.h | 8 +- .../property/elastic_tensor_decomposition.h | 31 +- include/aspect/particle/property/grain_size.h | 4 +- .../particle/property/integrated_strain.h | 9 +- .../property/integrated_strain_invariant.h | 9 +- include/aspect/particle/property/interface.h | 90 +++-- .../aspect/particle/property/melt_particle.h | 9 +- include/aspect/particle/property/pT_path.h | 9 +- include/aspect/particle/property/position.h | 9 +- .../particle/property/reference_position.h | 9 +- .../aspect/particle/property/strain_rate.h | 9 +- include/aspect/particle/property/velocity.h | 9 +- .../property/viscoplastic_strain_invariants.h | 9 +- include/aspect/particle/world.h | 11 +- include/aspect/solution_evaluator.h | 18 +- source/particle/property/composition.cc | 26 +- .../particle/property/cpo_bingham_average.cc | 42 +-- .../particle/property/cpo_elastic_tensor.cc | 24 +- .../property/crystal_preferred_orientation.cc | 309 +++++++++--------- source/particle/property/elastic_stress.cc | 43 +-- .../property/elastic_tensor_decomposition.cc | 104 +++--- source/particle/property/grain_size.cc | 42 ++- source/particle/property/integrated_strain.cc | 56 ++-- .../property/integrated_strain_invariant.cc | 44 +-- source/particle/property/interface.cc | 66 ++-- source/particle/property/melt_particle.cc | 19 +- source/particle/property/pT_path.cc | 15 +- source/particle/property/position.cc | 11 +- .../particle/property/reference_position.cc | 11 +- source/particle/property/strain_rate.cc | 29 +- source/particle/property/velocity.cc | 15 +- .../viscoplastic_strain_invariants.cc | 127 +++---- source/particle/world.cc | 52 +-- source/solution_evaluator.cc | 19 +- 39 files changed, 721 insertions(+), 686 deletions(-) create mode 100644 doc/modules/changes/20240925_gassmoeller diff --git a/doc/modules/changes/20240925_gassmoeller b/doc/modules/changes/20240925_gassmoeller new file mode 100644 index 00000000000..b955db31128 --- /dev/null +++ b/doc/modules/changes/20240925_gassmoeller @@ -0,0 +1,8 @@ +New: The interface for particle property plugins has been updated. +The new interface can update particles cell-wise and is more efficient. +In addition, the interface is more extensible for future improvements. +All particle property plugins in ASPECT have been updated, user plugins +will have to be updated accordingly before the support for the now deprecated +old interface is dropped. +
+(Rene Gassmoeller, 2024/09/25) diff --git a/include/aspect/particle/property/composition.h b/include/aspect/particle/property/composition.h index a115a811e3a..2fad04a302e 100644 --- a/include/aspect/particle/property/composition.h +++ b/include/aspect/particle/property/composition.h @@ -58,14 +58,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/cpo_bingham_average.h b/include/aspect/particle/property/cpo_bingham_average.h index a2ffd097350..ea2938dc681 100644 --- a/include/aspect/particle/property/cpo_bingham_average.h +++ b/include/aspect/particle/property/cpo_bingham_average.h @@ -99,36 +99,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * Update function. This function is called every time an update is - * request by need_update() for every particle for every property. - * It is obvious that - * this function is called a lot, so its code should be efficient. - * The interface provides a default implementation that does nothing, - * therefore derived plugins that do not require an update do not - * need to implement this function. - * - * @param [in] data_position An unsigned integer that denotes which - * component of the particle property vector is associated with the - * current property. For properties that own several components it - * denotes the first component of this property, all other components - * fill consecutive entries in the @p particle_properties vector. - * - * @param [in] solution The values of the solution variables at the - * current particle position. - * - * @param [in] gradients The gradients of the solution variables at - * the current particle position. - * - * @param [in,out] particle The particle that is updated within - * the call of this function. The particle location can be accessed - * using particle->get_location() and its properties using - * particle->get_properties(). + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/cpo_elastic_tensor.h b/include/aspect/particle/property/cpo_elastic_tensor.h index 1b7cc82a201..4ddaa2ba785 100644 --- a/include/aspect/particle/property/cpo_elastic_tensor.h +++ b/include/aspect/particle/property/cpo_elastic_tensor.h @@ -90,36 +90,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * Update function. This function is called every time an update is - * request by need_update() for every particle for every property. - * It is obvious that - * this function is called a lot, so its code should be efficient. - * The interface provides a default implementation that does nothing, - * therefore derived plugins that do not require an update do not - * need to implement this function. - * - * @param [in] data_position An unsigned integer that denotes which - * component of the particle property vector is associated with the - * current property. For properties that own several components it - * denotes the first component of this property, all other components - * fill consecutive entries in the @p particle_properties vector. - * - * @param [in] solution The values of the solution variables at the - * current particle position. - * - * @param [in] gradients The gradients of the solution variables at - * the current particle position. - * - * @param [in,out] particle The particle that is updated within - * the call of this function. The particle location can be accessed - * using particle->get_location() and its properties using - * particle->get_properties(). + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This function tells the particle manager that diff --git a/include/aspect/particle/property/crystal_preferred_orientation.h b/include/aspect/particle/property/crystal_preferred_orientation.h index 29c038f87c8..73346248798 100644 --- a/include/aspect/particle/property/crystal_preferred_orientation.h +++ b/include/aspect/particle/property/crystal_preferred_orientation.h @@ -158,36 +158,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * Update function. This function is called every time an update is - * request by need_update() for every particle for every property. - * It is obvious that - * this function is called a lot, so its code should be efficient. - * The interface provides a default implementation that does nothing, - * therefore derived plugins that do not require an update do not - * need to implement this function. - * - * @param [in] data_position An unsigned integer that denotes which - * component of the particle property vector is associated with the - * current property. For properties that own several components it - * denotes the first component of this property, all other components - * fill consecutive entries in the @p particle_properties vector. - * - * @param [in] solution The values of the solution variables at the - * current particle position. - * - * @param [in] gradients The gradients of the solution variables at - * the current particle position. - * - * @param [in,out] particle The particle that is updated within - * the call of this function. The particle location can be accessed - * using particle->get_location() and its properties using - * particle->get_properties(). + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/elastic_stress.h b/include/aspect/particle/property/elastic_stress.h index d4d70d7da7f..66f7f2579ec 100644 --- a/include/aspect/particle/property/elastic_stress.h +++ b/include/aspect/particle/property/elastic_stress.h @@ -58,13 +58,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * @copydoc aspect::Particle::Property::Interface::need_update() diff --git a/include/aspect/particle/property/elastic_tensor_decomposition.h b/include/aspect/particle/property/elastic_tensor_decomposition.h index 76ec5a314e0..faa78eef7b7 100644 --- a/include/aspect/particle/property/elastic_tensor_decomposition.h +++ b/include/aspect/particle/property/elastic_tensor_decomposition.h @@ -216,36 +216,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * Update function. This function is called every time an update is - * request by need_update() for every particle for every property. - * It is obvious that - * this function is called a lot, so its code should be efficient. - * The interface provides a default implementation that does nothing, - * therefore derived plugins that do not require an update do not - * need to implement this function. - * - * @param [in] data_position An unsigned integer that denotes which - * component of the particle property vector is associated with the - * current property. For properties that own several components it - * denotes the first component of this property, all other components - * fill consecutive entries in the @p particle_properties vector. - * - * @param [in] solution The values of the solution variables at the - * current particle position. - * - * @param [in] gradients The gradients of the solution variables at - * the current particle position. - * - * @param [in,out] particle The particle that is updated within - * the call of this function. The particle location can be accessed - * using particle->get_location() and its properties using - * particle->get_properties(). + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This function tells the particle manager that diff --git a/include/aspect/particle/property/grain_size.h b/include/aspect/particle/property/grain_size.h index bdfa42aac07..2c306be02b5 100644 --- a/include/aspect/particle/property/grain_size.h +++ b/include/aspect/particle/property/grain_size.h @@ -64,9 +64,7 @@ namespace aspect * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ void - update_particle_properties (const unsigned int data_position, - const std::vector> &solution, - const std::vector>> &gradients, + update_particle_properties (const ParticleUpdateInputs &inputs, typename ParticleHandler::particle_iterator_range &particles) const override; /** diff --git a/include/aspect/particle/property/integrated_strain.h b/include/aspect/particle/property/integrated_strain.h index e773dd9089f..41d10caca67 100644 --- a/include/aspect/particle/property/integrated_strain.h +++ b/include/aspect/particle/property/integrated_strain.h @@ -59,14 +59,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/integrated_strain_invariant.h b/include/aspect/particle/property/integrated_strain_invariant.h index 5ef0c8e78d9..9c969542961 100644 --- a/include/aspect/particle/property/integrated_strain_invariant.h +++ b/include/aspect/particle/property/integrated_strain_invariant.h @@ -58,14 +58,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/interface.h b/include/aspect/particle/property/interface.h index 7b6184847b1..46e734c7998 100644 --- a/include/aspect/particle/property/interface.h +++ b/include/aspect/particle/property/interface.h @@ -42,6 +42,36 @@ namespace aspect { using namespace dealii::Particles; + /** + * A data structure with all inputs for the + * Particle::Property::update_particle_properties() method. + */ + template + struct ParticleUpdateInputs + { + public: + /** + * The solution vector at each particle position. This vector is + * only filled if the update function requires the solution values. + */ + small_vector> solution; + + /** + * The solution gradients at each particle position. + * This vector is only filled if the update function requires the + * gradients of the solution values. + */ + small_vector>> gradients; + + /** + * Cell iterator of the cell that is currently being updated. + * This allows for evaluating additional properties at the cell vertices, + * or to query the cell for material ids, neighbors, or other + * information that is not available solely from the particles. + */ + typename DoFHandler::active_cell_iterator current_cell; + }; + /** * This class is used to store all the necessary information to translate * between the data structure of the particle properties (a flat vector of @@ -335,26 +365,18 @@ namespace aspect * therefore derived plugins that do not require an update do not * need to implement this function. * - * @param [in] data_position An unsigned integer that denotes which - * component of each particle property vector is associated with the - * current property. For properties that own several components it - * denotes the first component of this property, all other components - * fill consecutive entries in the properties vector. - * - * @param [in] solution A vector of values of the solution variables - * at the given particle positions. - * - * @param [in] gradients A vector of gradients of the solution - * variables at the given particle positions. + * @param [in] inputs A struct of type ParticleUpdateInputs that contains + * all necessary inputs to compute the particle updates. See + * the documentation of this struct in + * include/aspect/particle/property/interface.h for a list of all + * available inputs. * * @param [in,out] particles The particles that are to be updated * within this function. */ virtual void - update_particle_properties (const unsigned int data_position, - const std::vector> &solution, - const std::vector>> &gradients, + update_particle_properties (const ParticleUpdateInputs &inputs, typename ParticleHandler::particle_iterator_range &particles) const; /** @@ -383,7 +405,9 @@ namespace aspect * using particle->get_location() and its properties using * particle->get_properties(). * - * @deprecated Use update_particle_properties() instead. + * @deprecated This version of the function is deprecated. + * Use update_particle_properties() instead, which allows to + * update all particles of a cell in one function call. */ DEAL_II_DEPRECATED virtual @@ -456,6 +480,28 @@ namespace aspect virtual std::vector> get_property_information() const = 0; + + /** + * Set the position of this property in the particle property vector. + */ + virtual + void + set_data_position (const unsigned int data_position); + + /** + * Get the position of this property in the particle property vector. + */ + virtual + unsigned int + get_data_position () const; + + protected: + /** + * Store the position of the particle property in the particle property vector. + * If the property has multiple components, the first component is stored + * and all other components are stored consecutively after the first one. + */ + unsigned int data_position; }; /** @@ -547,17 +593,17 @@ namespace aspect * Update function for particle properties. This function is * called once every time step for every cell. * + * @param inputs A struct of type ParticleUpdateInputs that contains + * all necessary inputs to compute the particle updates. See + * the documentation of this struct in + * include/aspect/particle/property/interface.h for a list of all + * available inputs. * @param particles The particles that are to be updated within * this function. - * @param solution The values of the solution variables at the - * given particle positions. - * @param gradients The gradients of the solution variables at - * the given particle positions. */ void - update_particles (typename ParticleHandler::particle_iterator_range &particles, - const std::vector> &solution, - const std::vector>> &gradients) const; + update_particles (ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const; /** * Returns an enum, which denotes at what time this class needs to diff --git a/include/aspect/particle/property/melt_particle.h b/include/aspect/particle/property/melt_particle.h index 52a8ce8efbf..e9e4673f861 100644 --- a/include/aspect/particle/property/melt_particle.h +++ b/include/aspect/particle/property/melt_particle.h @@ -59,14 +59,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/pT_path.h b/include/aspect/particle/property/pT_path.h index 6c6331a883a..82f85c38957 100644 --- a/include/aspect/particle/property/pT_path.h +++ b/include/aspect/particle/property/pT_path.h @@ -70,14 +70,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/position.h b/include/aspect/particle/property/position.h index 67f24247df6..aa4a379108f 100644 --- a/include/aspect/particle/property/position.h +++ b/include/aspect/particle/property/position.h @@ -54,14 +54,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/reference_position.h b/include/aspect/particle/property/reference_position.h index 0de6a35f06b..e65c7c3cfcd 100644 --- a/include/aspect/particle/property/reference_position.h +++ b/include/aspect/particle/property/reference_position.h @@ -54,14 +54,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/strain_rate.h b/include/aspect/particle/property/strain_rate.h index d8c7460757d..98e42da387c 100644 --- a/include/aspect/particle/property/strain_rate.h +++ b/include/aspect/particle/property/strain_rate.h @@ -55,14 +55,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** diff --git a/include/aspect/particle/property/velocity.h b/include/aspect/particle/property/velocity.h index 70d710430bd..94272b035e8 100644 --- a/include/aspect/particle/property/velocity.h +++ b/include/aspect/particle/property/velocity.h @@ -55,14 +55,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * This implementation tells the particle manager that diff --git a/include/aspect/particle/property/viscoplastic_strain_invariants.h b/include/aspect/particle/property/viscoplastic_strain_invariants.h index d7a9c18da85..a0bc25986cf 100644 --- a/include/aspect/particle/property/viscoplastic_strain_invariants.h +++ b/include/aspect/particle/property/viscoplastic_strain_invariants.h @@ -61,14 +61,11 @@ namespace aspect std::vector &particle_properties) const override; /** - * @copydoc aspect::Particle::Property::Interface::update_particle_property() + * @copydoc aspect::Particle::Property::Interface::update_particle_properties() */ - virtual void - update_particle_property (const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const override; + update_particle_properties (const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const override; /** * @copydoc aspect::Particle::Property::Interface::need_update() diff --git a/include/aspect/particle/world.h b/include/aspect/particle/world.h index b4c797915fe..1a2bc6864cc 100644 --- a/include/aspect/particle/world.h +++ b/include/aspect/particle/world.h @@ -413,10 +413,17 @@ namespace aspect /** * Update the particle properties of one cell. + * + * @param inputs The input data required for the particle update. This + * function will fill this structure with the necessary data. + * @param positions The reference positions of the particles in the cell. + * This function will update these positions for the current cell. + * @param evaluator The solution evaluator that is used to update the particles. */ void - local_update_particles(const typename DoFHandler::active_cell_iterator &cell, - SolutionEvaluator &evaluators); + local_update_particles(Property::ParticleUpdateInputs &inputs, + small_vector> &positions, + SolutionEvaluator &evaluator); /** * Advect the particles of one cell. Performs only one step for diff --git a/include/aspect/solution_evaluator.h b/include/aspect/solution_evaluator.h index cd5e79fa851..4d5e64c2c07 100644 --- a/include/aspect/solution_evaluator.h +++ b/include/aspect/solution_evaluator.h @@ -166,18 +166,26 @@ namespace aspect * 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(). + * + * @param evaluation_point The index of the evaluation point in the positions array. + * @param solution The array to fill with the solution values. This array has to be + * of size n_components(). */ void get_solution(const unsigned int evaluation_point, - const ArrayView &solution); + const ArrayView &solution) const; /** * Fill @p gradients with all solution gradients at the given @p 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(). + * + * @param evaluation_point The index of the evaluation point in the positions array. + * @param gradients The array to fill with the solution gradients. This array has to be + * of size n_components(). */ void get_gradients(const unsigned int evaluation_point, - const ArrayView> &gradients); + const ArrayView> &gradients) const; /** * Return the evaluator for velocity or fluid velocity. This is the only @@ -192,6 +200,12 @@ namespace aspect NonMatching::MappingInfo & get_mapping_info(); + /** + * Return the number of components in the solution vector. + */ + unsigned int + n_components() const; + private: /** * MappingInfo object for the FEPointEvaluation objects diff --git a/source/particle/property/composition.cc b/source/particle/property/composition.cc index 8584d0f65ad..331a3f52a35 100644 --- a/source/particle/property/composition.cc +++ b/source/particle/property/composition.cc @@ -37,20 +37,27 @@ namespace aspect data.push_back(this->get_initial_composition_manager().initial_composition(position,i)); } + + template void - Composition::update_particle_property(const unsigned int data_position, - const Vector &solution, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + Composition::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - for (unsigned int i = 0; i < this->n_compositional_fields(); ++i) + unsigned int p = 0; + const auto &composition_components = this->introspection().component_indices.compositional_fields; + for (auto &particle: particles) { - const unsigned int solution_component = this->introspection().component_indices.compositional_fields[i]; - particle->get_properties()[data_position+i] = solution[solution_component]; + for (unsigned int j = 0; j < this->n_compositional_fields(); ++j) + { + particle.get_properties()[this->data_position+j] = inputs.solution[p][composition_components[j]]; + } + ++p; } } + + template UpdateTimeFlags Composition::need_update() const @@ -58,6 +65,8 @@ namespace aspect return update_time_step; } + + template UpdateFlags Composition::get_needed_update_flags () const @@ -65,6 +74,8 @@ namespace aspect return update_values; } + + template std::vector> Composition::get_property_information() const @@ -87,7 +98,6 @@ namespace aspect } return property_information; } - } } } diff --git a/source/particle/property/cpo_bingham_average.cc b/source/particle/property/cpo_bingham_average.cc index cfdfd423a32..956f683668e 100644 --- a/source/particle/property/cpo_bingham_average.cc +++ b/source/particle/property/cpo_bingham_average.cc @@ -81,31 +81,35 @@ namespace aspect template void - CpoBinghamAverage::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + CpoBinghamAverage::update_particle_properties(const ParticleUpdateInputs &/*inputs*/, + typename ParticleHandler::particle_iterator_range &particles) const { std::vector volume_fractions_grains(n_grains); std::vector> rotation_matrices_grains(n_grains); - ArrayView data = particle->get_properties(); - for (unsigned int mineral_i = 0; mineral_i < n_minerals; ++mineral_i) + + unsigned int p = 0; + for (auto &particle: particles) { - // create volume fractions and rotation matrix vectors in the order that it is stored in the data array - for (unsigned int grain_i = 0; grain_i < n_grains; ++grain_i) + ArrayView data = particle.get_properties(); + for (unsigned int mineral_i = 0; mineral_i < n_minerals; ++mineral_i) { - volume_fractions_grains[grain_i] = cpo_particle_property->get_volume_fractions_grains(cpo_data_position,data,mineral_i,grain_i); - rotation_matrices_grains[grain_i] = cpo_particle_property->get_rotation_matrix_grains(cpo_data_position,data,mineral_i,grain_i); + // create volume fractions and rotation matrix vectors in the order that it is stored in the data array + for (unsigned int grain_i = 0; grain_i < n_grains; ++grain_i) + { + volume_fractions_grains[grain_i] = cpo_particle_property->get_volume_fractions_grains(cpo_data_position,data,mineral_i,grain_i); + rotation_matrices_grains[grain_i] = cpo_particle_property->get_rotation_matrix_grains(cpo_data_position,data,mineral_i,grain_i); + } + + const std::vector> weighted_rotation_matrices = Utilities::rotation_matrices_random_draw_volume_weighting(volume_fractions_grains, rotation_matrices_grains, n_samples, this->random_number_generator); + std::array,3> bingham_average = compute_bingham_average(weighted_rotation_matrices); + + for (unsigned int i = 0; i < 3; ++i) + for (unsigned int j = 0; j < 6; ++j) + { + data[this->data_position + mineral_i*18 + i*6 + j] = bingham_average[i][j]; + } } - - const std::vector> weighted_rotation_matrices = Utilities::rotation_matrices_random_draw_volume_weighting(volume_fractions_grains, rotation_matrices_grains, n_samples, this->random_number_generator); - std::array,3> bingham_average = compute_bingham_average(weighted_rotation_matrices); - - for (unsigned int i = 0; i < 3; ++i) - for (unsigned int j = 0; j < 6; ++j) - { - data[data_position + mineral_i*18 + i*6 + j] = bingham_average[i][j]; - } + ++p; } } diff --git a/source/particle/property/cpo_elastic_tensor.cc b/source/particle/property/cpo_elastic_tensor.cc index 5a947da1d37..751b49639b2 100644 --- a/source/particle/property/cpo_elastic_tensor.cc +++ b/source/particle/property/cpo_elastic_tensor.cc @@ -140,25 +140,23 @@ namespace aspect template void - CpoElasticTensor::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + CpoElasticTensor::update_particle_properties(const ParticleUpdateInputs &/*inputs*/, + typename ParticleHandler::particle_iterator_range &particles) const { // Get a reference to the CPO particle property. const Particle::Property::CrystalPreferredOrientation &cpo_particle_property = this->get_particle_world(this->get_particle_world_index()).get_property_manager().template get_matching_active_plugin>(); + for (auto &particle: particles) + { + const SymmetricTensor<2,6> C_average = voigt_average_elastic_tensor(cpo_particle_property, + cpo_data_position, + particle.get_properties()); - const SymmetricTensor<2,6> C_average = voigt_average_elastic_tensor(cpo_particle_property, - cpo_data_position, - particle->get_properties()); - - Particle::Property::CpoElasticTensor::set_elastic_tensor(data_position, - particle->get_properties(), - C_average); - - + Particle::Property::CpoElasticTensor::set_elastic_tensor(this->data_position, + particle.get_properties(), + C_average); + } } diff --git a/source/particle/property/crystal_preferred_orientation.cc b/source/particle/property/crystal_preferred_orientation.cc index 58aab2dfdcb..73d011af993 100644 --- a/source/particle/property/crystal_preferred_orientation.cc +++ b/source/particle/property/crystal_preferred_orientation.cc @@ -224,179 +224,184 @@ namespace aspect template void - CrystalPreferredOrientation::update_particle_property(const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const + CrystalPreferredOrientation::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - // STEP 1: Load data and preprocess it. - - // need access to the pressure, viscosity, - // get velocity - Tensor<1,dim> velocity; - for (unsigned int i = 0; i < dim; ++i) - velocity[i] = solution[this->introspection().component_indices.velocities[i]]; - - // get velocity gradient tensor. - Tensor<2,dim> velocity_gradient; - for (unsigned int i = 0; i < dim; ++i) - velocity_gradient[i] = gradients[this->introspection().component_indices.velocities[i]]; - - // Calculate strain rate from velocity gradients - const SymmetricTensor<2,dim> strain_rate = symmetrize (velocity_gradient); - const SymmetricTensor<2,dim> deviatoric_strain_rate - = (this->get_material_model().is_compressible() - ? - strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor() - : - strain_rate); - - const double pressure = solution[this->introspection().component_indices.pressure]; - const double temperature = solution[this->introspection().component_indices.temperature]; - const double water_content = solution[this->introspection().component_indices.compositional_fields[water_index]]; - - // get the composition of the particle - std::vector compositions; - for (unsigned int i = 0; i < this->n_compositional_fields(); ++i) + const unsigned int data_position = this->data_position; + std::vector compositions(this->n_compositional_fields()); + + unsigned int p = 0; + for (auto &particle: particles) { - const unsigned int solution_component = this->introspection().component_indices.compositional_fields[i]; - compositions.push_back(solution[solution_component]); - } + // STEP 1: Load data and preprocess it. + + // need access to the pressure, viscosity, + // get velocity + Tensor<1,dim> velocity; + for (unsigned int i = 0; i < dim; ++i) + velocity[i] = inputs.solution[p][this->introspection().component_indices.velocities[i]]; + + // get velocity gradient tensor. + Tensor<2,dim> velocity_gradient; + for (unsigned int i = 0; i < dim; ++i) + velocity_gradient[i] = inputs.gradients[p][this->introspection().component_indices.velocities[i]]; + + // Calculate strain rate from velocity gradients + const SymmetricTensor<2,dim> strain_rate = symmetrize (velocity_gradient); + const SymmetricTensor<2,dim> deviatoric_strain_rate + = (this->get_material_model().is_compressible() + ? + strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor() + : + strain_rate); + + const double pressure = inputs.solution[p][this->introspection().component_indices.pressure]; + const double temperature = inputs.solution[p][this->introspection().component_indices.temperature]; + const double water_content = inputs.solution[p][this->introspection().component_indices.compositional_fields[water_index]]; + + // get the composition of the particle + for (unsigned int i = 0; i < this->n_compositional_fields(); ++i) + { + const unsigned int solution_component = this->introspection().component_indices.compositional_fields[i]; + compositions[i] = inputs.solution[p][solution_component]; + } - const double dt = this->get_timestep(); + const double dt = this->get_timestep(); - // even in 2d we need 3d strain-rates and velocity gradient tensors. So we make them 3d by - // adding an extra dimension which is zero. - SymmetricTensor<2,3> strain_rate_3d; - strain_rate_3d[0][0] = strain_rate[0][0]; - strain_rate_3d[0][1] = strain_rate[0][1]; - //sym: strain_rate_3d[1][0] = strain_rate[1][0]; - strain_rate_3d[1][1] = strain_rate[1][1]; + // even in 2d we need 3d strain-rates and velocity gradient tensors. So we make them 3d by + // adding an extra dimension which is zero. + SymmetricTensor<2,3> strain_rate_3d; + strain_rate_3d[0][0] = strain_rate[0][0]; + strain_rate_3d[0][1] = strain_rate[0][1]; + //sym: strain_rate_3d[1][0] = strain_rate[1][0]; + strain_rate_3d[1][1] = strain_rate[1][1]; - if (dim == 3) - { - strain_rate_3d[0][2] = strain_rate[0][2]; - strain_rate_3d[1][2] = strain_rate[1][2]; - //sym: strain_rate_3d[2][0] = strain_rate[0][2]; - //sym: strain_rate_3d[2][1] = strain_rate[1][2]; - strain_rate_3d[2][2] = strain_rate[2][2]; - } - Tensor<2,3> velocity_gradient_3d; - velocity_gradient_3d[0][0] = velocity_gradient[0][0]; - velocity_gradient_3d[0][1] = velocity_gradient[0][1]; - velocity_gradient_3d[1][0] = velocity_gradient[1][0]; - velocity_gradient_3d[1][1] = velocity_gradient[1][1]; - if (dim == 3) - { - velocity_gradient_3d[0][2] = velocity_gradient[0][2]; - velocity_gradient_3d[1][2] = velocity_gradient[1][2]; - velocity_gradient_3d[2][0] = velocity_gradient[2][0]; - velocity_gradient_3d[2][1] = velocity_gradient[2][1]; - velocity_gradient_3d[2][2] = velocity_gradient[2][2]; - } - - ArrayView data = particle->get_properties(); + if (dim == 3) + { + strain_rate_3d[0][2] = strain_rate[0][2]; + strain_rate_3d[1][2] = strain_rate[1][2]; + //sym: strain_rate_3d[2][0] = strain_rate[0][2]; + //sym: strain_rate_3d[2][1] = strain_rate[1][2]; + strain_rate_3d[2][2] = strain_rate[2][2]; + } + Tensor<2,3> velocity_gradient_3d; + velocity_gradient_3d[0][0] = velocity_gradient[0][0]; + velocity_gradient_3d[0][1] = velocity_gradient[0][1]; + velocity_gradient_3d[1][0] = velocity_gradient[1][0]; + velocity_gradient_3d[1][1] = velocity_gradient[1][1]; + if (dim == 3) + { + velocity_gradient_3d[0][2] = velocity_gradient[0][2]; + velocity_gradient_3d[1][2] = velocity_gradient[1][2]; + velocity_gradient_3d[2][0] = velocity_gradient[2][0]; + velocity_gradient_3d[2][1] = velocity_gradient[2][1]; + velocity_gradient_3d[2][2] = velocity_gradient[2][2]; + } - for (unsigned int mineral_i = 0; mineral_i < n_minerals; ++mineral_i) - { + ArrayView data = particle.get_properties(); - /** - * Now we have loaded all the data and can do the actual computation. - * The computation consists of two parts. The first part is computing - * the derivatives for the directions and grain sizes. Then those - * derivatives are used to advect the particle properties. - */ - double sum_volume_mineral = 0; - std::pair, std::vector>> - derivatives_grains = this->compute_derivatives(data_position, - data, - mineral_i, - strain_rate_3d, - velocity_gradient_3d, - particle->get_location(), - temperature, - pressure, - velocity, - compositions, - strain_rate, - deviatoric_strain_rate, - water_content); - - switch (advection_method) + for (unsigned int mineral_i = 0; mineral_i < n_minerals; ++mineral_i) { - case AdvectionMethod::forward_euler: - sum_volume_mineral = this->advect_forward_euler(data_position, - data, - mineral_i, - dt, - derivatives_grains); + /** + * Now we have loaded all the data and can do the actual computation. + * The computation consists of two parts. The first part is computing + * the derivatives for the directions and grain sizes. Then those + * derivatives are used to advect the particle properties. + */ + double sum_volume_mineral = 0; + std::pair, std::vector>> + derivatives_grains = this->compute_derivatives(data_position, + data, + mineral_i, + strain_rate_3d, + velocity_gradient_3d, + particle.get_location(), + temperature, + pressure, + velocity, + compositions, + strain_rate, + deviatoric_strain_rate, + water_content); + + switch (advection_method) + { + case AdvectionMethod::forward_euler: - break; + sum_volume_mineral = this->advect_forward_euler(data_position, + data, + mineral_i, + dt, + derivatives_grains); - case AdvectionMethod::backward_euler: - sum_volume_mineral = this->advect_backward_euler(data_position, - data, - mineral_i, - dt, - derivatives_grains); + break; - break; - } + case AdvectionMethod::backward_euler: + sum_volume_mineral = this->advect_backward_euler(data_position, + data, + mineral_i, + dt, + derivatives_grains); - // normalize the volume fractions back to a total of 1 for each mineral - const double inv_sum_volume_mineral = 1.0/sum_volume_mineral; + break; + } - Assert(std::isfinite(inv_sum_volume_mineral), - ExcMessage("inv_sum_volume_mineral is not finite. sum_volume_enstatite = " - + std::to_string(sum_volume_mineral))); + // normalize the volume fractions back to a total of 1 for each mineral + const double inv_sum_volume_mineral = 1.0/sum_volume_mineral; - for (unsigned int grain_i = 0; grain_i < n_grains; ++grain_i) - { - const double volume_fraction_grains = get_volume_fractions_grains(data_position,data,mineral_i,grain_i)*inv_sum_volume_mineral; - set_volume_fractions_grains(data_position,data,mineral_i,grain_i,volume_fraction_grains); - Assert(isfinite(get_volume_fractions_grains(data_position,data,mineral_i,grain_i)), - ExcMessage("volume_fractions_grains[mineral_i]" + std::to_string(grain_i) + "] is not finite: " - + std::to_string(get_volume_fractions_grains(data_position,data,mineral_i,grain_i)) + ", inv_sum_volume_mineral = " - + std::to_string(inv_sum_volume_mineral) + ".")); + Assert(std::isfinite(inv_sum_volume_mineral), + ExcMessage("inv_sum_volume_mineral is not finite. sum_volume_enstatite = " + + std::to_string(sum_volume_mineral))); - /** - * Correct direction rotation matrices numerical error (orthnormality) after integration - * Follows same method as in matlab version from Thissen (see https://github.com/cthissen/Drex-MATLAB/) - * of finding the nearest orthonormal matrix using the SVD - */ - Tensor<2,3> rotation_matrix = get_rotation_matrix_grains(data_position,data,mineral_i,grain_i); - for (size_t i = 0; i < 3; ++i) + for (unsigned int grain_i = 0; grain_i < n_grains; ++grain_i) { - for (size_t j = 0; j < 3; ++j) + const double volume_fraction_grains = get_volume_fractions_grains(data_position,data,mineral_i,grain_i)*inv_sum_volume_mineral; + set_volume_fractions_grains(data_position,data,mineral_i,grain_i,volume_fraction_grains); + Assert(isfinite(get_volume_fractions_grains(data_position,data,mineral_i,grain_i)), + ExcMessage("volume_fractions_grains[mineral_i]" + std::to_string(grain_i) + "] is not finite: " + + std::to_string(get_volume_fractions_grains(data_position,data,mineral_i,grain_i)) + ", inv_sum_volume_mineral = " + + std::to_string(inv_sum_volume_mineral) + ".")); + + /** + * Correct direction rotation matrices numerical error (orthnormality) after integration + * Follows same method as in matlab version from Thissen (see https://github.com/cthissen/Drex-MATLAB/) + * of finding the nearest orthonormal matrix using the SVD + */ + Tensor<2,3> rotation_matrix = get_rotation_matrix_grains(data_position,data,mineral_i,grain_i); + for (size_t i = 0; i < 3; ++i) { - Assert(!std::isnan(rotation_matrix[i][j]), ExcMessage("rotation_matrix is nan before orthogonalization.")); + for (size_t j = 0; j < 3; ++j) + { + Assert(!std::isnan(rotation_matrix[i][j]), ExcMessage("rotation_matrix is nan before orthogonalization.")); + } } - } - rotation_matrix = dealii::project_onto_orthogonal_tensors(rotation_matrix); - for (size_t i = 0; i < 3; ++i) - for (size_t j = 0; j < 3; ++j) - { - // I don't think this should happen with the projection, but D-Rex - // does not do the orthogonal projection, but just clamps the values - // to 1 and -1. - Assert(std::fabs(rotation_matrix[i][j]) <= 1.0, - ExcMessage("The rotation_matrix has a entry larger than 1.")); - - Assert(!std::isnan(rotation_matrix[i][j]), - ExcMessage("rotation_matrix is nan after orthoganalization: " - + std::to_string(rotation_matrix[i][j]))); - - Assert(abs(rotation_matrix[i][j]) <= 1.0, - ExcMessage("3. rotation_matrix[" + std::to_string(i) + "][" + std::to_string(j) + - "] is larger than one: " - + std::to_string(rotation_matrix[i][j]) + " (" + std::to_string(rotation_matrix[i][j]-1.0) + "). rotation_matrix = \n" - + std::to_string(rotation_matrix[0][0]) + " " + std::to_string(rotation_matrix[0][1]) + " " + std::to_string(rotation_matrix[0][2]) + "\n" - + std::to_string(rotation_matrix[1][0]) + " " + std::to_string(rotation_matrix[1][1]) + " " + std::to_string(rotation_matrix[1][2]) + "\n" - + std::to_string(rotation_matrix[2][0]) + " " + std::to_string(rotation_matrix[2][1]) + " " + std::to_string(rotation_matrix[2][2]))); - } + rotation_matrix = dealii::project_onto_orthogonal_tensors(rotation_matrix); + for (size_t i = 0; i < 3; ++i) + for (size_t j = 0; j < 3; ++j) + { + // I don't think this should happen with the projection, but D-Rex + // does not do the orthogonal projection, but just clamps the values + // to 1 and -1. + Assert(std::fabs(rotation_matrix[i][j]) <= 1.0, + ExcMessage("The rotation_matrix has a entry larger than 1.")); + + Assert(!std::isnan(rotation_matrix[i][j]), + ExcMessage("rotation_matrix is nan after orthoganalization: " + + std::to_string(rotation_matrix[i][j]))); + + Assert(abs(rotation_matrix[i][j]) <= 1.0, + ExcMessage("3. rotation_matrix[" + std::to_string(i) + "][" + std::to_string(j) + + "] is larger than one: " + + std::to_string(rotation_matrix[i][j]) + " (" + std::to_string(rotation_matrix[i][j]-1.0) + "). rotation_matrix = \n" + + std::to_string(rotation_matrix[0][0]) + " " + std::to_string(rotation_matrix[0][1]) + " " + std::to_string(rotation_matrix[0][2]) + "\n" + + std::to_string(rotation_matrix[1][0]) + " " + std::to_string(rotation_matrix[1][1]) + " " + std::to_string(rotation_matrix[1][2]) + "\n" + + std::to_string(rotation_matrix[2][0]) + " " + std::to_string(rotation_matrix[2][1]) + " " + std::to_string(rotation_matrix[2][2]))); + } + } } + ++p; } } diff --git a/source/particle/property/elastic_stress.cc b/source/particle/property/elastic_stress.cc index f0be9a72362..49d811ca9ef 100644 --- a/source/particle/property/elastic_stress.cc +++ b/source/particle/property/elastic_stress.cc @@ -89,36 +89,39 @@ namespace aspect template void - ElasticStress::update_particle_property(const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const + ElasticStress::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - material_inputs.position[0] = particle->get_location(); + unsigned int p = 0; + for (auto &particle: particles) + { + material_inputs.position[0] = particle.get_location(); + + material_inputs.current_cell = inputs.current_cell; - material_inputs.current_cell = typename DoFHandler::active_cell_iterator(*particle->get_surrounding_cell(), - &(this->get_dof_handler())); + material_inputs.temperature[0] = inputs.solution[p][this->introspection().component_indices.temperature]; - material_inputs.temperature[0] = solution[this->introspection().component_indices.temperature]; + material_inputs.pressure[0] = inputs.solution[p][this->introspection().component_indices.pressure]; - material_inputs.pressure[0] = solution[this->introspection().component_indices.pressure]; + for (unsigned int d = 0; d < dim; ++d) + material_inputs.velocity[0][d] = inputs.solution[p][this->introspection().component_indices.velocities[d]]; - for (unsigned int d = 0; d < dim; ++d) - material_inputs.velocity[0][d] = solution[this->introspection().component_indices.velocities[d]]; + for (unsigned int n = 0; n < this->n_compositional_fields(); ++n) + material_inputs.composition[0][n] = inputs.solution[p][this->introspection().component_indices.compositional_fields[n]]; - for (unsigned int n = 0; n < this->n_compositional_fields(); ++n) - material_inputs.composition[0][n] = solution[this->introspection().component_indices.compositional_fields[n]]; + Tensor<2,dim> grad_u; + for (unsigned int d=0; d grad_u; - for (unsigned int d=0; dget_material_model().evaluate (material_inputs,material_outputs); - this->get_material_model().evaluate (material_inputs,material_outputs); + for (unsigned int i = 0; i < SymmetricTensor<2,dim>::n_independent_components ; ++i) + particle.get_properties()[this->data_position + i] += material_outputs.reaction_terms[0][i]; - for (unsigned int i = 0; i < SymmetricTensor<2,dim>::n_independent_components ; ++i) - particle->get_properties()[data_position + i] += material_outputs.reaction_terms[0][i]; + ++p; + } } diff --git a/source/particle/property/elastic_tensor_decomposition.cc b/source/particle/property/elastic_tensor_decomposition.cc index 25308ff2623..b8466e93fb2 100644 --- a/source/particle/property/elastic_tensor_decomposition.cc +++ b/source/particle/property/elastic_tensor_decomposition.cc @@ -347,60 +347,62 @@ namespace aspect template void - ElasticTensorDecomposition::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + ElasticTensorDecomposition::update_particle_properties(const ParticleUpdateInputs &/*inputs*/, + typename ParticleHandler::particle_iterator_range &particles) const { - ArrayView data = particle->get_properties(); - const SymmetricTensor<2,6> elastic_matrix = Particle::Property::CpoElasticTensor::get_elastic_tensor(cpo_elastic_tensor_data_position, - data); - - const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = Utilities::compute_dilatation_stiffness_tensor(elastic_matrix); - const SymmetricTensor<2,3> voigt_stiffness_tensor_full = Utilities::compute_voigt_stiffness_tensor(elastic_matrix); - const Tensor<2,3> SCCS_full = Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); - - const std::array,7 > norms = Utilities::compute_elastic_tensor_SCCS_decompositions(SCCS_full, elastic_matrix); - - // get max hexagonal element index, which is the same as the permutation index - const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); - std::array permutation = Utilities::indexed_even_permutation(max_hexagonal_element_index); - // reorder the SCCS by the SCCS permutation which yields the largest hexagonal vector (percentage of anisotropy) - Tensor<2,3> hexagonal_permutated_SCCS; - for (size_t index = 0; index < 3; ++index) + const unsigned int data_position = this->data_position; + for (auto &particle: particles) { - hexagonal_permutated_SCCS[index] = SCCS_full[permutation[index]]; + ArrayView data = particle.get_properties(); + const SymmetricTensor<2,6> elastic_matrix = Particle::Property::CpoElasticTensor::get_elastic_tensor(cpo_elastic_tensor_data_position, + data); + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = Utilities::compute_dilatation_stiffness_tensor(elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = Utilities::compute_voigt_stiffness_tensor(elastic_matrix); + const Tensor<2,3> SCCS_full = Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + + const std::array,7 > norms = Utilities::compute_elastic_tensor_SCCS_decompositions(SCCS_full, elastic_matrix); + + // get max hexagonal element index, which is the same as the permutation index + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = Utilities::indexed_even_permutation(max_hexagonal_element_index); + // reorder the SCCS by the SCCS permutation which yields the largest hexagonal vector (percentage of anisotropy) + Tensor<2,3> hexagonal_permutated_SCCS; + for (size_t index = 0; index < 3; ++index) + { + hexagonal_permutated_SCCS[index] = SCCS_full[permutation[index]]; + } + + data[data_position] = SCCS_full[0][0]; + data[data_position+1] = SCCS_full[0][1]; + data[data_position+2] = SCCS_full[0][2]; + data[data_position+3] = SCCS_full[1][0]; + data[data_position+4] = SCCS_full[1][1]; + data[data_position+5] = SCCS_full[1][2]; + data[data_position+6] = SCCS_full[2][0]; + data[data_position+7] = SCCS_full[2][1]; + data[data_position+8] = SCCS_full[2][2]; + data[data_position+9] = hexagonal_permutated_SCCS[2][0]; + data[data_position+10] = hexagonal_permutated_SCCS[2][1]; + data[data_position+11] = hexagonal_permutated_SCCS[2][2]; + data[data_position+12] = norms[6][0]; + data[data_position+13] = norms[0][0]; // triclinic + data[data_position+14] = norms[0][1]; // triclinic + data[data_position+15] = norms[0][2]; // triclinic + data[data_position+16] = norms[1][0]; // monoclinic + data[data_position+17] = norms[1][1]; // monoclinic + data[data_position+18] = norms[1][2]; // monoclinic + data[data_position+19] = norms[2][0]; // orthorhomic + data[data_position+20] = norms[2][1]; // orthorhomic + data[data_position+21] = norms[2][2]; // orthorhomic + data[data_position+22] = norms[3][0]; // tetragonal + data[data_position+23] = norms[3][1]; // tetragonal + data[data_position+24] = norms[3][2]; // tetragonal + data[data_position+25] = norms[4][0]; // hexagonal + data[data_position+26] = norms[4][1]; // hexagonal + data[data_position+27] = norms[4][2]; // hexagonal + data[data_position+28] = norms[5][0]; // isotropic } - - data[data_position] = SCCS_full[0][0]; - data[data_position+1] = SCCS_full[0][1]; - data[data_position+2] = SCCS_full[0][2]; - data[data_position+3] = SCCS_full[1][0]; - data[data_position+4] = SCCS_full[1][1]; - data[data_position+5] = SCCS_full[1][2]; - data[data_position+6] = SCCS_full[2][0]; - data[data_position+7] = SCCS_full[2][1]; - data[data_position+8] = SCCS_full[2][2]; - data[data_position+9] = hexagonal_permutated_SCCS[2][0]; - data[data_position+10] = hexagonal_permutated_SCCS[2][1]; - data[data_position+11] = hexagonal_permutated_SCCS[2][2]; - data[data_position+12] = norms[6][0]; - data[data_position+13] = norms[0][0]; // triclinic - data[data_position+14] = norms[0][1]; // triclinic - data[data_position+15] = norms[0][2]; // triclinic - data[data_position+16] = norms[1][0]; // monoclinic - data[data_position+17] = norms[1][1]; // monoclinic - data[data_position+18] = norms[1][2]; // monoclinic - data[data_position+19] = norms[2][0]; // orthorhomic - data[data_position+20] = norms[2][1]; // orthorhomic - data[data_position+21] = norms[2][2]; // orthorhomic - data[data_position+22] = norms[3][0]; // tetragonal - data[data_position+23] = norms[3][1]; // tetragonal - data[data_position+24] = norms[3][2]; // tetragonal - data[data_position+25] = norms[4][0]; // hexagonal - data[data_position+26] = norms[4][1]; // hexagonal - data[data_position+27] = norms[4][2]; // hexagonal - data[data_position+28] = norms[5][0]; // isotropic } diff --git a/source/particle/property/grain_size.cc b/source/particle/property/grain_size.cc index 145d00f267c..eb82d821d36 100644 --- a/source/particle/property/grain_size.cc +++ b/source/particle/property/grain_size.cc @@ -66,52 +66,50 @@ namespace aspect template void - GrainSize::update_particle_properties(const unsigned int data_position, - const std::vector> &solution, - const std::vector>> &gradients, + GrainSize::update_particle_properties(const ParticleUpdateInputs &inputs, typename ParticleHandler::particle_iterator_range &particles) const { - material_inputs = MaterialModel::MaterialModelInputs(solution.size(), this->n_compositional_fields()); - material_outputs = MaterialModel::MaterialModelOutputs(solution.size(), this->n_compositional_fields()); + material_inputs = MaterialModel::MaterialModelInputs(inputs.solution.size(), this->n_compositional_fields()); + material_outputs = MaterialModel::MaterialModelOutputs(inputs.solution.size(), this->n_compositional_fields()); material_inputs.requested_properties = MaterialModel::MaterialProperties::reaction_terms; - material_inputs.current_cell = typename DoFHandler::active_cell_iterator(*particles.begin()->get_surrounding_cell(), - &(this->get_dof_handler())); + material_inputs.current_cell = inputs.current_cell; - unsigned int i = 0; + unsigned int p = 0; for (auto particle: particles) { // Make sure all particles are in the same cell - Assert(particle.get_surrounding_cell() == particles.begin()->get_surrounding_cell(), + Assert(particle.get_surrounding_cell() == inputs.current_cell, ExcMessage("All particles must be in the same cell.")); - material_inputs.position[i] = particle.get_location(); - material_inputs.temperature[i] = solution[i][this->introspection().component_indices.temperature]; - material_inputs.pressure[i] = solution[i][this->introspection().component_indices.pressure]; + material_inputs.position[p] = particle.get_location(); + material_inputs.temperature[p] = inputs.solution[p][this->introspection().component_indices.temperature]; + material_inputs.pressure[p] = inputs.solution[p][this->introspection().component_indices.pressure]; for (unsigned int d = 0; d < dim; ++d) - material_inputs.velocity[i][d] = solution[i][this->introspection().component_indices.velocities[d]]; + material_inputs.velocity[p][d] = inputs.solution[p][this->introspection().component_indices.velocities[d]]; for (unsigned int n = 0; n < this->n_compositional_fields(); ++n) - material_inputs.composition[i][n] = solution[i][this->introspection().component_indices.compositional_fields[n]]; + material_inputs.composition[p][n] = inputs.solution[p][this->introspection().component_indices.compositional_fields[n]]; - material_inputs.composition[i][grain_size_index] = particle.get_properties()[data_position]; + material_inputs.composition[p][grain_size_index] = particle.get_properties()[this->data_position]; Tensor<2,dim> grad_u; for (unsigned int d=0; dintrospection().component_indices.velocities[d]]; - ++i; + material_inputs.strain_rate[p] = symmetrize (grad_u); + + ++p; } this->get_material_model().evaluate(material_inputs, material_outputs); - i = 0; - for (auto particle: particles) + p = 0; + for (auto &particle: particles) { - particle.get_properties()[data_position] += material_outputs.reaction_terms[i][grain_size_index]; - ++i; + particle.get_properties()[this->data_position] += material_outputs.reaction_terms[p][grain_size_index]; + ++p; } } diff --git a/source/particle/property/integrated_strain.cc b/source/particle/property/integrated_strain.cc index 7afa788bad3..7f386a90adb 100644 --- a/source/particle/property/integrated_strain.cc +++ b/source/particle/property/integrated_strain.cc @@ -38,43 +38,45 @@ namespace aspect template void - IntegratedStrain::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const + IntegratedStrain::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - const Tensor<2,dim> old_strain(make_array_view(&particle->get_properties()[data_position], - &particle->get_properties()[data_position] + Tensor<2,dim>::n_independent_components)); + const double dt = this->get_timestep(); - Tensor<2,dim> grad_u; - for (unsigned int d=0; d old_strain(make_array_view(&particle.get_properties()[this->data_position], + &particle.get_properties()[this->data_position] + Tensor<2,dim>::n_independent_components)); - const double dt = this->get_timestep(); + Tensor<2,dim> grad_u; + for (unsigned int d=0; dintrospection().component_indices.velocities[d]]; - Tensor<2,dim> new_strain; + // here we integrate the equation + // new_deformation_gradient = velocity_gradient * old_deformation_gradient + // using a RK4 integration scheme. + const Tensor<2,dim> k1 = grad_u * old_strain * dt; + Tensor<2,dim> new_strain = old_strain + 0.5*k1; - // here we integrate the equation - // new_deformation_gradient = velocity_gradient * old_deformation_gradient - // using a RK4 integration scheme. - const Tensor<2,dim> k1 = grad_u * old_strain * dt; - new_strain = old_strain + 0.5*k1; + const Tensor<2,dim> k2 = grad_u * new_strain * dt; + new_strain = old_strain + 0.5*k2; - const Tensor<2,dim> k2 = grad_u * new_strain * dt; - new_strain = old_strain + 0.5*k2; + const Tensor<2,dim> k3 = grad_u * new_strain * dt; + new_strain = old_strain + k3; - const Tensor<2,dim> k3 = grad_u * new_strain * dt; - new_strain = old_strain + k3; + const Tensor<2,dim> k4 = grad_u * new_strain * dt; - const Tensor<2,dim> k4 = grad_u * new_strain * dt; + // the new strain is the rotated old strain plus the + // strain of the current time step + new_strain = old_strain + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0; - // the new strain is the rotated old strain plus the - // strain of the current time step - new_strain = old_strain + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0; + // unroll and store the new strain + new_strain.unroll(&particle.get_properties()[this->data_position], + &particle.get_properties()[this->data_position] + Tensor<2,dim>::n_independent_components); - // unroll and store the new strain - new_strain.unroll(&particle->get_properties()[data_position], - &particle->get_properties()[data_position] + Tensor<2,dim>::n_independent_components); + ++p; + } } template diff --git a/source/particle/property/integrated_strain_invariant.cc b/source/particle/property/integrated_strain_invariant.cc index a0d13bb0f82..9688f9a0d66 100644 --- a/source/particle/property/integrated_strain_invariant.cc +++ b/source/particle/property/integrated_strain_invariant.cc @@ -38,32 +38,36 @@ namespace aspect template void - IntegratedStrainInvariant::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const + IntegratedStrainInvariant::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - // Integrated strain invariant from prior time step - const auto data = particle->get_properties(); - double old_strain = data[data_position]; - // Current timestep - const double dt = this->get_timestep(); + unsigned int p = 0; + for (auto &particle: particles) + { + // Integrated strain invariant from prior time step + const auto data = particle.get_properties(); + const double old_strain = data[this->data_position]; - // Velocity gradients - Tensor<2,dim> grad_u; - for (unsigned int d=0; dget_timestep(); - // Calculate strain rate from velocity gradients - const SymmetricTensor<2,dim> strain_rate = symmetrize (grad_u); + // Velocity gradients + Tensor<2,dim> grad_u; + for (unsigned int d=0; d strain_rate = symmetrize (grad_u); - // New strain is the old strain plus dt*edot_ii - const double new_strain = old_strain + dt*edot_ii; - data[data_position] = new_strain; + // Calculate strain rate second invariant + const double edot_ii = std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)); + + // New strain is the old strain plus dt*edot_ii + const double new_strain = old_strain + dt*edot_ii; + data[this->data_position] = new_strain; + ++p; + } } diff --git a/source/particle/property/interface.cc b/source/particle/property/interface.cc index 564a7a0dfd1..0cac16f4bcf 100644 --- a/source/particle/property/interface.cc +++ b/source/particle/property/interface.cc @@ -209,25 +209,36 @@ namespace aspect DEAL_II_DISABLE_EXTRA_DIAGNOSTICS template void - Interface::update_particle_properties (const unsigned int data_position, - const std::vector> &solution, - const std::vector>> &gradients, + Interface::update_particle_properties (const ParticleUpdateInputs &inputs, typename ParticleHandler::particle_iterator_range &particles) const { - const Vector invalid_solution; - const std::vector> invalid_gradient; + Vector solution; + std::vector> gradient; unsigned int i = 0; - for (typename ParticleHandler::particle_iterator particle = particles.begin(); - particle != particles.end(); ++particle, ++i) + for (auto particle = particles.begin(); particle != particles.end(); ++particle) { + if (get_needed_update_flags() & update_values) + { + solution.reinit(inputs.solution[i].size()); + for (unsigned int j=0; j0) ? - solution[i] : invalid_solution, - (gradients.size()>0) ? - gradients[i] : invalid_gradient, + update_particle_property(this->data_position, + solution, + gradient, particle); + + ++i; } } DEAL_II_ENABLE_EXTRA_DIAGNOSTICS @@ -271,6 +282,24 @@ namespace aspect + template + void + Interface::set_data_position (const unsigned int index) + { + data_position = index; + } + + + + template + unsigned int + Interface::get_data_position () const + { + return data_position; + } + + + template void IntegratorProperties::initialize_one_particle_property(const Point &/*position*/, @@ -325,9 +354,12 @@ namespace aspect // Initialize our property information property_information = ParticlePropertyInformation(info); + unsigned int plugin_index = 0; for (const auto &p : this->plugin_objects) { + p->set_data_position(property_information.get_position_by_plugin_index(plugin_index)); p->initialize(); + ++plugin_index; } } @@ -518,18 +550,14 @@ namespace aspect template void - Manager::update_particles (typename ParticleHandler::particle_iterator_range &particles, - const std::vector> &solution, - const std::vector>> &gradients) const + Manager::update_particles (ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { unsigned int plugin_index = 0; for (typename std::list>>::const_iterator p = this->plugin_objects.begin(); p!=this->plugin_objects.end(); ++p,++plugin_index) { - (*p)->update_particle_properties(property_information.get_position_by_plugin_index(plugin_index), - solution, - gradients, - particles); + (*p)->update_particle_properties(inputs,particles); } } diff --git a/source/particle/property/melt_particle.cc b/source/particle/property/melt_particle.cc index fa80b32bee3..c99f8b27e6b 100644 --- a/source/particle/property/melt_particle.cc +++ b/source/particle/property/melt_particle.cc @@ -39,20 +39,23 @@ namespace aspect template void - MeltParticle::update_particle_property(const unsigned int data_position, - const Vector &solution, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + MeltParticle::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { AssertThrow(this->introspection().compositional_name_exists("porosity"), ExcMessage("Particle property melt particle only works if" "there is a compositional field called porosity.")); const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity"); - if (solution[this->introspection().component_indices.compositional_fields[porosity_idx]] > threshold_for_melt_presence) - particle->get_properties()[data_position] = 1.0; - else - particle->get_properties()[data_position] = 0.0; + unsigned int p = 0; + for (auto &particle: particles) + { + if (inputs.solution[p][this->introspection().component_indices.compositional_fields[porosity_idx]] > threshold_for_melt_presence) + particle.get_properties()[this->data_position] = 1.0; + else + particle.get_properties()[this->data_position] = 0.0; + ++p; + } } template diff --git a/source/particle/property/pT_path.cc b/source/particle/property/pT_path.cc index ffff2519a9e..88f6ba75cb2 100644 --- a/source/particle/property/pT_path.cc +++ b/source/particle/property/pT_path.cc @@ -65,13 +65,16 @@ namespace aspect template void - PTPath::update_particle_property(const unsigned int data_position, - const Vector &solution, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + PTPath::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - particle->get_properties()[data_position] = solution[this->introspection().component_indices.pressure]; - particle->get_properties()[data_position+1] = solution[this->introspection().component_indices.temperature]; + unsigned int p = 0; + for (auto &particle: particles) + { + particle.get_properties()[this->data_position] = inputs.solution[p][this->introspection().component_indices.pressure]; + particle.get_properties()[this->data_position+1] = inputs.solution[p][this->introspection().component_indices.temperature]; + ++p; + } } diff --git a/source/particle/property/position.cc b/source/particle/property/position.cc index a21a0837482..171194fb519 100644 --- a/source/particle/property/position.cc +++ b/source/particle/property/position.cc @@ -37,13 +37,12 @@ namespace aspect template void - Position::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + Position::update_particle_properties(const ParticleUpdateInputs &/*inputs*/, + typename ParticleHandler::particle_iterator_range &particles) const { - for (unsigned int i = 0; i < dim; ++i) - particle->get_properties()[data_position+i] = particle->get_location()[i]; + for (auto &particle: particles) + for (unsigned int i = 0; i < dim; ++i) + particle.get_properties()[this->data_position+i] = particle.get_location()[i]; } template diff --git a/source/particle/property/reference_position.cc b/source/particle/property/reference_position.cc index 6d5fa5f978e..a9d5614839f 100644 --- a/source/particle/property/reference_position.cc +++ b/source/particle/property/reference_position.cc @@ -37,13 +37,12 @@ namespace aspect template void - ReferencePosition::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + ReferencePosition::update_particle_properties(const ParticleUpdateInputs &/*inputs*/, + typename ParticleHandler::particle_iterator_range &particles) const { - for (unsigned int i = 0; i < dim; ++i) - particle->get_properties()[data_position+i] = particle->get_reference_location()[i]; + for (auto &particle: particles) + for (unsigned int i = 0; i < dim; ++i) + particle.get_properties()[this->data_position+i] = particle.get_reference_location()[i]; } template diff --git a/source/particle/property/strain_rate.cc b/source/particle/property/strain_rate.cc index be79a182387..78c32f2bf16 100644 --- a/source/particle/property/strain_rate.cc +++ b/source/particle/property/strain_rate.cc @@ -41,22 +41,25 @@ namespace aspect template void - StrainRate::update_particle_property(const unsigned int data_position, - const Vector &/*solution*/, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const + StrainRate::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - const auto data = particle->get_properties(); - // Velocity gradients - Tensor<2,dim> grad_u; - for (unsigned int d=0; d grad_u; + for (unsigned int d=0; d strain_rate = symmetrize (grad_u); + // Calculate strain rate from velocity gradients + const SymmetricTensor<2,dim> strain_rate = symmetrize (grad_u); - for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i) - data[data_position + i] = strain_rate[Tensor<2,dim>::unrolled_to_component_indices(i)]; + for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i) + data[this->data_position + i] = strain_rate[Tensor<2,dim>::unrolled_to_component_indices(i)]; + ++p; + } } diff --git a/source/particle/property/velocity.cc b/source/particle/property/velocity.cc index 82270e8db5b..329f6df140d 100644 --- a/source/particle/property/velocity.cc +++ b/source/particle/property/velocity.cc @@ -37,13 +37,16 @@ namespace aspect template void - Velocity::update_particle_property(const unsigned int data_position, - const Vector &solution, - const std::vector> &/*gradients*/, - typename ParticleHandler::particle_iterator &particle) const + Velocity::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { - for (unsigned int i = 0; i < dim; ++i) - particle->get_properties()[data_position+i] = solution[this->introspection().component_indices.velocities[i]]; + unsigned int p = 0; + for (auto &particle: particles) + { + for (unsigned int i = 0; i < dim; ++i) + particle.get_properties()[this->data_position+i] = inputs.solution[p][this->introspection().component_indices.velocities[i]]; + ++p; + } } template diff --git a/source/particle/property/viscoplastic_strain_invariants.cc b/source/particle/property/viscoplastic_strain_invariants.cc index 8ea89a2cda2..6df9095373a 100644 --- a/source/particle/property/viscoplastic_strain_invariants.cc +++ b/source/particle/property/viscoplastic_strain_invariants.cc @@ -97,85 +97,92 @@ namespace aspect template void - ViscoPlasticStrainInvariant::update_particle_property(const unsigned int data_position, - const Vector &solution, - const std::vector> &gradients, - typename ParticleHandler::particle_iterator &particle) const + ViscoPlasticStrainInvariant::update_particle_properties(const ParticleUpdateInputs &inputs, + typename ParticleHandler::particle_iterator_range &particles) const { + // Find out plastic yielding by calling function in material model. + const MaterialModel::ViscoPlastic &viscoplastic + = Plugins::get_plugin_as_type>(this->get_material_model()); + // Current timestep const double dt = this->get_timestep(); + const unsigned int data_position = this->data_position; - // Velocity gradients - Tensor<2,dim> grad_u; - for (unsigned int d=0; dintrospection().component_indices.pressure]; - material_inputs.temperature[0] = solution[this->introspection().component_indices.temperature]; - material_inputs.position[0] = particle->get_location(); + // Velocity gradients + Tensor<2,dim> grad_u; + for (unsigned int d=0; dintrospection().component_indices.pressure]; + material_inputs.temperature[0] = inputs.solution[p][this->introspection().component_indices.temperature]; + material_inputs.position[0] = particle.get_location(); - // Put compositional fields into single variable - for (unsigned int i = 0; i < this->n_compositional_fields(); ++i) - { - material_inputs.composition[0][i] = solution[this->introspection().component_indices.compositional_fields[i]]; - } + // Calculate strain rate from velocity gradients + material_inputs.strain_rate[0] = symmetrize (grad_u); - // Find out plastic yielding by calling function in material model. - const MaterialModel::ViscoPlastic &viscoplastic - = Plugins::get_plugin_as_type>(this->get_material_model()); + // Put compositional fields into single variable + for (unsigned int i = 0; i < this->n_compositional_fields(); ++i) + { + material_inputs.composition[0][i] = inputs.solution[p][this->introspection().component_indices.compositional_fields[i]]; + } - const bool plastic_yielding = viscoplastic.is_yielding(material_inputs); - // Next take the integrated strain invariant from the prior time step. - const auto data = particle->get_properties(); + const bool plastic_yielding = viscoplastic.is_yielding(material_inputs); - // Calculate strain rate second invariant - const double edot_ii = std::sqrt(std::max(-second_invariant(deviator(material_inputs.strain_rate[0])), 0.)); + // Next take the integrated strain invariant from the prior time step. + const auto data = particle.get_properties(); - // Calculate strain invariant magnitude over the last time step - const double strain_update = dt*edot_ii; + // Calculate strain rate second invariant + const double edot_ii = std::sqrt(std::max(-second_invariant(deviator(material_inputs.strain_rate[0])), 0.)); - /* Update the strain values that are used in the simulation, which use the following assumptions - * to identify the correct position in the data vector for each value: - * (1) Total strain cannot be used in combination with any other strain field - * (2) If plastic strain is tracked, it will always be in the first data position - * (3) If noninitial plastic strain is tracked, it will always be in the last data position - * (4) If noninitial plastic strain is tracked, plastic strain is also being tracked - * (5) If only viscous strain is tracked, it will be in the first data position. - * (6) If both viscous and plastic strain are tracked, viscous strain will be in the second data position - * If these assumptions change in the future, they will need to be updated. - * */ + // Calculate strain invariant magnitude over the last time step + const double strain_update = dt*edot_ii; - if (this->introspection().compositional_name_exists("plastic_strain") && plastic_yielding == true) - data[data_position] += strain_update; + /* Update the strain values that are used in the simulation, which use the following assumptions + * to identify the correct position in the data vector for each value: + * (1) Total strain cannot be used in combination with any other strain field + * (2) If plastic strain is tracked, it will always be in the first data position + * (3) If noninitial plastic strain is tracked, it will always be in the last data position + * (4) If noninitial plastic strain is tracked, plastic strain is also being tracked + * (5) If only viscous strain is tracked, it will be in the first data position. + * (6) If both viscous and plastic strain are tracked, viscous strain will be in the second data position + * If these assumptions change in the future, they will need to be updated. + * */ - if (this->introspection().compositional_name_exists("viscous_strain") && plastic_yielding == false) - { - // Not yielding and only one field, which tracks the viscous strain. - if (n_components == 1) + if (this->introspection().compositional_name_exists("plastic_strain") && plastic_yielding == true) data[data_position] += strain_update; - // Not yielding and either two or three fields are tracked. If two fields are tracked, - // they represent plastic strain (first data position) and viscous strain (second data - // data position, updated below). If three fields are tracked, they represent plastic - // strain (first data position), viscous strain (second data position, updated below), - // and noninitial plastic strain (third data position). In either case, the viscous - // strain is in the second data position, allowing us to use a single expression. - if (n_components > 1) - data[data_position+1] += strain_update; - } + if (this->introspection().compositional_name_exists("viscous_strain") && plastic_yielding == false) + { + // Not yielding and only one field, which tracks the viscous strain. + if (n_components == 1) + data[data_position] += strain_update; + + // Not yielding and either two or three fields are tracked. If two fields are tracked, + // they represent plastic strain (first data position) and viscous strain (second data + // data position, updated below). If three fields are tracked, they represent plastic + // strain (first data position), viscous strain (second data position, updated below), + // and noninitial plastic strain (third data position). In either case, the viscous + // strain is in the second data position, allowing us to use a single expression. + if (n_components > 1) + data[data_position+1] += strain_update; + } + + // Only one field, which tracks total strain and is updated regardless of whether the + // material is yielding or not. + if (this->introspection().compositional_name_exists("total_strain")) + data[data_position] += strain_update; - // Only one field, which tracks total strain and is updated regardless of whether the - // material is yielding or not. - if (this->introspection().compositional_name_exists("total_strain")) - data[data_position] += strain_update; + // Yielding, and noninitial plastic strain (last data position, updated below) is tracked. + if (this->introspection().compositional_name_exists("noninitial_plastic_strain") && plastic_yielding == true) + data[data_position+(n_components-1)] += strain_update; - // Yielding, and noninitial plastic strain (last data position, updated below) is tracked. - if (this->introspection().compositional_name_exists("noninitial_plastic_strain") && plastic_yielding == true) - data[data_position+(n_components-1)] += strain_update; + ++p; + } } diff --git a/source/particle/world.cc b/source/particle/world.cc index 09f0910b0bd..5126d4a31ec 100644 --- a/source/particle/world.cc +++ b/source/particle/world.cc @@ -443,51 +443,56 @@ namespace aspect template void - World::local_update_particles(const typename DoFHandler::active_cell_iterator &cell, + World::local_update_particles(Property::ParticleUpdateInputs &inputs, + small_vector> &positions, SolutionEvaluator &evaluator) { - const unsigned int n_particles_in_cell = particle_handler->n_particles_in_cell(cell); - typename ParticleHandler::particle_iterator_range particles = particle_handler->particles_in_cell(cell); + const unsigned int n_particles = particle_handler->n_particles_in_cell(inputs.current_cell); - std::vector> positions; - positions.reserve(n_particles_in_cell); + typename ParticleHandler::particle_iterator_range particles = particle_handler->particles_in_cell(inputs.current_cell); + positions.resize(n_particles); + unsigned int p = 0; for (const auto &particle : particles) - positions.push_back(particle.get_reference_location()); + { + positions[p] = particle.get_reference_location(); + ++p; + } const UpdateFlags update_flags = property_manager->get_needed_update_flags(); small_vector solution_values(this->get_fe().dofs_per_cell); - cell->get_dof_values(this->get_solution(), - solution_values.begin(), - solution_values.end()); + inputs.current_cell->get_dof_values(this->get_solution(), + solution_values.begin(), + solution_values.end()); if (update_flags & (update_values | update_gradients)) - evaluator.reinit(cell, positions, {solution_values.data(), solution_values.size()}, update_flags); + { + evaluator.reinit(inputs.current_cell, + {positions.data(), positions.size()}, + {solution_values.data(), solution_values.size()}, + update_flags); + } - std::vector> solution; if (update_flags & update_values) - solution.resize(n_particles_in_cell,Vector(this->introspection().n_components)); + inputs.solution.resize(n_particles,small_vector(evaluator.n_components())); - std::vector>> gradients; if (update_flags & update_gradients) - gradients.resize(n_particles_in_cell,std::vector>(this->introspection().n_components)); + inputs.gradients.resize(n_particles,small_vector>(evaluator.n_components())); - for (unsigned int i = 0; iupdate_particles(particles, - solution, - gradients); + property_manager->update_particles(inputs,particles); } @@ -631,6 +636,9 @@ namespace aspect std::unique_ptr> evaluator = construct_solution_evaluator(*this, update_flags); + Property::ParticleUpdateInputs inputs; + small_vector> positions; + // Loop over all cells and update the particles cell-wise for (const auto &cell : this->get_dof_handler().active_cell_iterators()) if (cell->is_locally_owned()) @@ -638,7 +646,9 @@ namespace aspect // Only update particles if there are any in this cell if (particle_handler->n_particles_in_cell(cell) > 0) { - local_update_particles(cell, + inputs.current_cell = cell; + local_update_particles(inputs, + positions, *evaluator); } diff --git a/source/solution_evaluator.cc b/source/solution_evaluator.cc index a63de8be6e1..50d869c44a9 100644 --- a/source/solution_evaluator.cc +++ b/source/solution_evaluator.cc @@ -383,7 +383,7 @@ namespace aspect template void SolutionEvaluator::get_solution(const unsigned int evaluation_point, - const ArrayView &solution) + const ArrayView &solution) const { Assert(solution.size() == simulator_access.introspection().n_components, ExcDimensionMismatch(solution.size(), simulator_access.introspection().n_components)); @@ -421,7 +421,7 @@ namespace aspect template void SolutionEvaluator::get_gradients(const unsigned int evaluation_point, - const ArrayView> &gradients) + const ArrayView> &gradients) const { Assert(gradients.size() == simulator_access.introspection().n_components, ExcDimensionMismatch(gradients.size(), simulator_access.introspection().n_components)); @@ -456,6 +456,7 @@ namespace aspect } + template FEPointEvaluation & SolutionEvaluator::get_velocity_or_fluid_velocity_evaluator(const bool use_fluid_velocity) @@ -467,6 +468,9 @@ namespace aspect return velocity; } + + + template NonMatching::MappingInfo & SolutionEvaluator::get_mapping_info() @@ -475,7 +479,16 @@ namespace aspect } - // A function to create a pointer to a SolutionEvaluator object. + + template + unsigned int + SolutionEvaluator::n_components() const + { + return simulator_access.introspection().n_components; + } + + + template std::unique_ptr> construct_solution_evaluator (const SimulatorAccess &simulator_access, From e7f4b2d54a0e78f29497cc80495e676b68698ee4 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Fri, 27 Sep 2024 11:35:12 +0200 Subject: [PATCH 2/2] Reduce statically allocated memory size --- include/aspect/particle/property/interface.h | 4 ++-- source/particle/world.cc | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/aspect/particle/property/interface.h b/include/aspect/particle/property/interface.h index 46e734c7998..ca37116a274 100644 --- a/include/aspect/particle/property/interface.h +++ b/include/aspect/particle/property/interface.h @@ -54,14 +54,14 @@ namespace aspect * The solution vector at each particle position. This vector is * only filled if the update function requires the solution values. */ - small_vector> solution; + std::vector> solution; /** * The solution gradients at each particle position. * This vector is only filled if the update function requires the * gradients of the solution values. */ - small_vector>> gradients; + std::vector,50>> gradients; /** * Cell iterator of the cell that is currently being updated. diff --git a/source/particle/world.cc b/source/particle/world.cc index 5126d4a31ec..01e079231b2 100644 --- a/source/particle/world.cc +++ b/source/particle/world.cc @@ -476,10 +476,10 @@ namespace aspect } if (update_flags & update_values) - inputs.solution.resize(n_particles,small_vector(evaluator.n_components())); + inputs.solution.resize(n_particles,small_vector(evaluator.n_components())); if (update_flags & update_gradients) - inputs.gradients.resize(n_particles,small_vector>(evaluator.n_components())); + inputs.gradients.resize(n_particles,small_vector,50>(evaluator.n_components())); for (unsigned int i = 0; i