allow different base_elements for compositional fields #5802

tjhei · 2024-06-05T21:51:45Z

We now store different base_element indices in introspection (one index for each compositional field). This is no functional change, as we do not support different FiniteElements so far.

To be able to make this change, I had to fix various places that use base_elements or assume that all compositional fields have the same base element.

part of #5748

tjhei · 2024-06-05T22:23:11Z

This one if for you @bangerth . :-)

bangerth · 2024-06-05T23:17:44Z

There are numerous errors such as these:

In file included from /Users/jenkins/jenkins/workspace/aspect-osx_PR-5802/build/CMakeFiles/aspect-release.dir/Unity/unity_40_cxx.cxx:16:
�[1m/Users/jenkins/jenkins/workspace/aspect-osx_PR-5802/source/simulator/helper_functions.cc:1690:74: �[0m�[0;1;31merror: �[0m�[1mno viable conversion from 'const std::vector<unsigned int>' to 'unsigned int'�[0m
    const Quadrature<dim> quadrature_C(dof_handler.get_fe().base_element(introspection.base_elements.compositional_fields).get_unit_support_points());
�[0;1;32m                                                                         ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
�[0m�[1m/Users/jenkins/jenkins/workspace/aspect-osx_PR-5802/source/simulator/helper_functions.cc:2578:22: �[0m�[0;1;30mnote: �[0min instantiation of member function 'aspect::Simulator<2>::compute_reactions' requested here�[0m
  ASPECT_INSTANTIATE(INSTANTIATE)
�[0;1;32m                     ^
�[0m�[1m/Users/jenkins/candi-9.5.1-r1b/deal.II-v9.5.1/include/deal.II/fe/fe.h:1622:35: �[0m�[0;1;30mnote: �[0mpassing argument to parameter 'index' here�[0m
  base_element(const unsigned int index) const;
�[0;1;32m                                  ^
�[0mIn file included from /Users/jenkins/jenkins/workspace/aspect-osx_PR-5802/build/CMakeFiles/aspect-release.dir/Unity/unity_40_cxx.cxx:16:
�[1m/Users/jenkins/jenkins/workspace/aspect-osx_PR-5802/source/simulator/helper_functions.cc:1785:74: �[0m�[0;1;31merror: �[0m�[1mno viable conversion from 'const std::vector<unsigned int>' to 'unsigned int'�[0m
              for (unsigned int j=0; j<dof_handler.get_fe().base_element(introspection.base_elements.compositional_fields).dofs_per_cell; ++j)
�[0;1;32m                                                                         ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
�[0m�[1m/Users/jenkins/candi-9.5.1-r1b/deal.II-v9.5.1/include/deal.II/fe/fe.h:1622:35: �[0m�[0;1;30mnote: �[0mpassing argument to parameter 'index' here�[0m
  base_element(const unsigned int index) const;
�[0;1;32m                                  ^
�[0mIn file included from /Users/jenkins/jenkins/workspace/aspect-osx_PR-5802/build/CMakeFiles/aspect-release.dir/Unity/unity_40_cxx.cxx:16:
�[1m/Users/jenkins/jenkins/workspace/aspect-osx_PR-5802/source/simulator/helper_functions.cc:1826:70: �[0m�[0;1;31merror: �[0m�[1mno viable conversion from 'const std::vector<unsigned int>' to 'unsigned int'�[0m
          for (unsigned int j=0; j<dof_handler.get_fe().base_element(introspection.base_elements.compositional_fields).dofs_per_cell; ++j)
�[0;1;32m                                                                     ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

I would assume that you get those on your machine too?

bangerth

First part.

bangerth · 2024-06-06T04:32:44Z

include/aspect/introspection.h

+      /**
+       * Return a vector with all used indices of base elements of the deal.II FiniteElement
+       * that are used for compositional fields. If several fields are the same type, they
+       * share base elements.
+      */
+      std::vector<unsigned int>
+      get_compositional_field_base_element_indices() const;


Am I correct that the output array has length # of compositional fields? If so, would you mind adding that.

No, currently it has length 1, which I thought is obvious from the comment. I will make it more clear.

bangerth · 2024-06-06T04:34:04Z

include/aspect/introspection.h

+      /**
+       * Return a vector with all compositional field indices that belong to a given
+       * base element index.
+      */
+      std::vector<unsigned int>
+      get_compositional_field_indices_with_base_element(const unsigned int base_element_index) const;


Is this indexed starting at the zeroth compositional element? Or at the velocity?

Or perhaps you are saying that the argument must be an element of the array that get_compositional_field_base_element_indices() returns? If so, please add that. (And assert?)

bangerth · 2024-06-06T04:35:05Z

source/mesh_refinement/compaction_length.cc

-            for (unsigned int i=0; i<this->get_fe().base_element(this->introspection().base_elements.compositional_fields).dofs_per_cell; ++i)
+            for (unsigned int i=0; i<in.n_evaluation_points(); ++i)


I would accept this commit separately if it's cumbersome to carry around.

source/mesh_refinement/composition_gradient.cc

source/simulator/introspection.cc

tjhei · 2024-06-06T17:08:05Z

Updated and rebased.

bangerth

Close.

include/aspect/introspection.h

bangerth · 2024-06-06T21:23:52Z

source/simulator/helper_functions.cc

+
+          for (unsigned int i=0; i<points.size(); ++i)
+            {
+              const auto &it = std::find(unique_support_points.begin(), unique_support_points.end(), points[i]);


You don't want a reference to a temporary.

Suggested change

const auto &it = std::find(unique_support_points.begin(), unique_support_points.end(), points[i]);

const auto it = std::find(unique_support_points.begin(), unique_support_points.end(), points[i]);

I think lifetime extension makes this work, but I agree.

source/simulator/helper_functions.cc

jdannberg

I think this is way better than it was before, I really like it! In a few places, I found it hard to follow what all the different indices are, but it's a clear improvement on my original version.

jdannberg · 2024-06-06T20:49:24Z

source/mesh_refinement/composition_gradient.cc

+                    // refining, so anything >1 should work. (note that the gradient
+                    // itself scales like 1/h, so multiplying it with any factor h^s, s>1
+                    // will yield convergence of the error indicators to zero as h->0)
+                    this_indicator *= std::pow(cell->diameter(), power);


Is this still the same as before? You're not resetting your this_indicator variable in the loop. So it would add up from the different compositional fields with the same base index.

jdannberg · 2024-06-06T22:53:20Z

source/simulator/helper_functions.cc

-    HeatingModel::HeatingModelOutputs heating_model_outputs_C(quadrature_C.size(), introspection.n_compositional_fields);
+    // First compute the union of all support points (temperature and compositions):
+    std::vector<Point<dim>> unique_support_points;
+    // For each field (T or composition) store a list of indices of support points (unique_support_points):


So these are, for each support point of a given field, the index that this support point has in the unique_support_points vector? I think that did not quite become clear to me from reading just the comment.

jdannberg · 2024-06-06T22:56:08Z

source/simulator/helper_functions.cc

    // on every cell, compute the update, and then on every cell put the result into the
    // distributed_vector vector. Only after the loop over all cells do we copy distributed_vector
    // back onto the solution vector.
    // So even though we touch some DoF more than once, we always start from the same value, compute the
    // same value, and then overwrite the same value in distributed_vector.
-    // TODO: make this more efficient.
+    // TODO: make this even more efficient.


jdannberg · 2024-06-06T23:04:05Z

source/simulator/introspection.cc

+  Introspection<dim>::is_composition_component (const unsigned int component_index) const
+  {
+    // All compositions live at the end. Just to be sure, there are no other components
+    // in our system after compositional fields, right?


Does that work in models with melt transport? We do have the melt velocity, fluid pressure etc. before the temperature, right?

Yes, that is correct. But I added an Assert here just in case somebody messes with the order later on.

jdannberg · 2024-06-06T23:14:35Z

source/simulator/helper_functions.cc

+                  // the how-manyth dof is that for the current component and where do we find this in the list
+                  // of evaluation points?
+                  const unsigned int index_within = comp_pair.second;
+                  const unsigned int field_index = component_idx-component_idx_T;
+                  const unsigned int idx = support_point_index_by_field[field_index][index_within];


I find this hard to follow. I understand the general concept, but if you could add another sentence of explanation (either here, or where you fill the support_point_index_by_field vector) I think that would be helpful.

tjhei · 2024-06-06T23:53:22Z

Updated. Thank you for the comments! 👍

tjhei · 2024-06-11T19:39:36Z

Fingers crossed, this might be good to go now.

tjhei · 2024-06-11T20:19:53Z

The two tests are changing very slightly due to the reordering of the loops in the composition mesh refinement plugin.

gassmoeller

Looks all correct, but could you look into my comments about the documentation? I think the current text is slightly confusing. Also some small suggestions for performance improvements.

gassmoeller · 2024-06-11T21:32:47Z

include/aspect/introspection.h

+       * Return a vector with all used indices of base elements of the deal.II FiniteElement
+       * that are used for compositional fields. If several fields are the same type, they


This documentation is not clear to me. Do you mean: Return a vector with as many entries as compositional fields. Each entry of the vector represents the index of the base element that is used for this compositional field.?

EDIT: After reading the function I understood. What about rewording to: Return a vector that contains the base element indices of the deal.II FiniteElement that are used for compositional fields. Note that compositional fields can share the same base element, so this vector can (and usually will) be smaller than the number of compositional fields. The function get_compositional_field_indices_with_base_element() can be used to translate from base element index to compositional field indices.

Also in the second sentence the phrase are the same type suggests a relation with the compositional_field_types feature, but you probably mean If several fields use the same finite element and degree, they share base elements.

gassmoeller · 2024-06-12T14:01:44Z

source/mesh_refinement/composition_gradient.cc

@@ -19,6 +19,7 @@
 */


+#include <aspect/simulator.h>


Ah, we really need to find a solution to put AdvectionField out of the simulator.h header. Including simulator.h in a lot of places is not a clean design. Not for this PR though.

agreed, see #5887

gassmoeller · 2024-06-12T14:14:39Z

source/simulator/helper_functions.cc

-              if (dof_handler.locally_owned_dofs().is_element(local_dof_indices[temperature_idx]))
+              const auto comp_pair = dof_handler.get_fe().system_to_component_index(dof_idx);
+              const unsigned int component_idx = comp_pair.first;
+              if (component_idx>=component_idx_T) // ignore velocity, pressure, etc.


This relies on implicit knowledge about the component ordering. What about:

Suggested change

if (component_idx>=component_idx_T) // ignore velocity, pressure, etc.

if (introspection.is_stokes_component(component_idx) == false) // ignore velocity, pressure, etc.

Neat idea, but this won't work with melt or other FE modifications, I think.

source/simulator/helper_functions.cc

source/simulator/introspection.cc

gassmoeller · 2024-06-12T14:30:38Z

source/simulator/introspection.cc

+    std::vector<unsigned int> result;
+    const unsigned int first = this->base_elements.compositional_fields[0];
+    const unsigned int last = *std::max_element(this->base_elements.compositional_fields.begin(), this->base_elements.compositional_fields.end());
+    for (unsigned int i = first; i<=last; ++i)
+      result.emplace_back(i);
+    return result;


Save a few allocations:

Suggested change

std::vector<unsigned int> result;

const unsigned int first = this->base_elements.compositional_fields[0];

const unsigned int last = *std::max_element(this->base_elements.compositional_fields.begin(), this->base_elements.compositional_fields.end());

for (unsigned int i = first; i<=last; ++i)

result.emplace_back(i);

return result;

const unsigned int first = this->base_elements.compositional_fields[0];

const unsigned int last = *std::max_element(this->base_elements.compositional_fields.begin(), this->base_elements.compositional_fields.end());

std::vector<unsigned int> result(last-first);

std::iota(result.begin(), result.end(), first);

return result;

Even better: I reworked some parts of introspection to cache this information (see get_names_for_fields_of_type()). This way we can avoid allocating a vector in tight loops (the caller of this function might not know it is expensive).

gassmoeller · 2024-06-12T14:52:33Z

source/simulator/introspection.cc

+    std::vector<unsigned int> result;
+    Assert(base_element_index <= get_compositional_field_base_element_indices().back(),
+           ExcMessage("Invalid base_element_index specified."));
+
+    unsigned int idx = 0;
+    for (const auto base_idx : this->base_elements.compositional_fields)
+      {
+        if (base_idx == base_element_index)
+          result.emplace_back(idx);
+
+        ++idx;
+      }


Also here, I think it would be better if this could be computed in the constructor once and then just returned by reference (it is called from a loop over all base elements x over all cells). I think of introspection as a sort of cache for internal information so it would be ok to precompute things.

are you ok with a std::map? See the last commit.

I am ok with it until it shows up in a profiling run ;-) ... we can keep it for now. But can you change both functions to return references to vectors? Also I think you need to indent.

tjhei · 2024-06-12T19:51:27Z

updated.

tjhei added the work in progress label Jun 5, 2024

tjhei force-pushed the base_elements branch 3 times, most recently from fb601c0 to 7fdefbc Compare June 6, 2024 03:57

bangerth reviewed Jun 6, 2024

View reviewed changes

source/simulator/introspection.cc Outdated Show resolved Hide resolved

source/simulator/introspection.cc Outdated Show resolved Hide resolved

source/simulator/introspection.cc Outdated Show resolved Hide resolved

tjhei force-pushed the base_elements branch from e99756a to e5e9ae7 Compare June 6, 2024 13:50

tjhei mentioned this pull request Jun 6, 2024

simplify compaction length refinement #5816

Merged

tjhei force-pushed the base_elements branch 2 times, most recently from 7774d04 to 86d521f Compare June 6, 2024 17:06

tjhei force-pushed the base_elements branch from 86d521f to 3642cf9 Compare June 6, 2024 17:37

tjhei mentioned this pull request Jun 6, 2024

simplify code #5832

Merged

bangerth reviewed Jun 6, 2024

View reviewed changes

jdannberg approved these changes Jun 6, 2024

View reviewed changes

tjhei force-pushed the base_elements branch 3 times, most recently from 12f4781 to 7831a79 Compare June 11, 2024 18:04

allow different base_elements for compositional fields

e32ba66

tjhei force-pushed the base_elements branch from 7831a79 to e32ba66 Compare June 11, 2024 18:41

update test output

a0ec101

tjhei removed the work in progress label Jun 11, 2024

gassmoeller reviewed Jun 12, 2024

View reviewed changes

tjhei mentioned this pull request Jun 12, 2024

move AdvectionField out of Simulator #5887

Open

address comments

e1d8b0f

cache in introspection

f1551e9

tjhei force-pushed the base_elements branch from ae7b9fe to f1551e9 Compare June 12, 2024 18:55

gassmoeller approved these changes Jun 13, 2024

View reviewed changes

gassmoeller merged commit c6f4ed7 into geodynamics:main Jun 13, 2024
8 checks passed

gassmoeller mentioned this pull request Jun 13, 2024

Return const reference #5893

Merged

		for (unsigned int i=0; i<this->get_fe().base_element(this->introspection().base_elements.compositional_fields).dofs_per_cell; ++i)
		for (unsigned int i=0; i<in.n_evaluation_points(); ++i)

	const auto &it = std::find(unique_support_points.begin(), unique_support_points.end(), points[i]);
	const auto it = std::find(unique_support_points.begin(), unique_support_points.end(), points[i]);

		* Return a vector with all used indices of base elements of the deal.II FiniteElement
		* that are used for compositional fields. If several fields are the same type, they

	if (component_idx>=component_idx_T) // ignore velocity, pressure, etc.
	if (introspection.is_stokes_component(component_idx) == false) // ignore velocity, pressure, etc.

allow different base_elements for compositional fields #5802

allow different base_elements for compositional fields #5802

Conversation

tjhei commented Jun 5, 2024 • edited Loading

tjhei commented Jun 5, 2024

bangerth commented Jun 5, 2024

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tjhei Jun 6, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tjhei commented Jun 6, 2024

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jdannberg left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tjhei Jun 6, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tjhei commented Jun 6, 2024 • edited Loading

tjhei commented Jun 11, 2024

tjhei commented Jun 11, 2024

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tjhei commented Jun 12, 2024

tjhei commented Jun 5, 2024 •

edited

Loading

tjhei Jun 6, 2024 •

edited

Loading

tjhei Jun 6, 2024 •

edited

Loading

tjhei commented Jun 6, 2024 •

edited

Loading