Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Combine interpolating fields in one operation #6202

Merged
merged 1 commit into from
Mar 2, 2025

Conversation

gassmoeller
Copy link
Member

@Wang-yijun and @KiralyAgi showed me a model with 36 prescribed compositional fields (a lot of components of a larger tensor), in which the interpolation from material model output into compositional fields took >50% of the total model runtime. This made me look into the algorithm and search for optimization potential. When we prescribe compositional fields from the material model output we currently do the following steps consecutively for each compositional field:

  • evaluate the solution at the support points of the current compositional field
  • evaluate the material model at the support points incl. prescribed outputs for all fields
  • only write the material model output for the current field into the solution vector

All steps of this operation could be combined for all compositional fields that need to be interpolated, which would save significant time for models with many prescribed fields and expensive material models. This PR implements this change by copying the approach we already take for particles. We essentially collect all fields that are prescribed while looping through the compositional fields and then interpolate all of them together after all other fields have been solved.

Some intricacies:

  • the order in which we solve/interpolate compositional fields now changes if prescribed fields are mixed in between "normal" fields, but this doesnt matter much, because the current linearization point for all fields is only updated once all fields are solved (to avoid inconsistent time stepping / reactions between them). So the reordering of fields will only affect the screen and statistics output
  • the fields may have different finite elements. I solved this the same way we do in the function compute_reactions, by collecting all support points of all fields and creating a combined quadrature for all of them. I refactored the code from compute_reactions into its own function to avoid duplicating it.
  • if we interpolate several compositional fields at the same time we need to change the screen output, it now no longer prints the field name

First tests by @Wang-yijun (and I could reproduce this) showed that this change reduces the time spent in this algorithm from >50% to <5% of the total model time for her model.

Copy link
Member

@bobmyhill bobmyhill left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice! Just a couple of comments.

for (const auto &support_point: support_points)
{
// Naive n^2 algorithm to merge all points into the data structure, speed is likely irrelevant, because
// n is small and this is only done rarely.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
// n is small and this is only done rarely.
// n is small and this function is rarely called.

@@ -8,8 +8,8 @@
# 8: Number of nonlinear iterations
# 9: Iterations for temperature solver
# 10: Iterations for composition solver 1
# 11: Iterations for composition solver 2
# 12: Iterations for composition solver 3
# 11: Iterations for composition solver 3
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the reason for this change in ordering?

@@ -23,5 +23,5 @@
# 23: Outward heat flux through boundary with indicator 1 ("top") (W)
# 24: Outward heat flux through boundary with indicator 2 ("left") (W)
# 25: Outward heat flux through boundary with indicator 3 ("right") (W)
0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499 0 0 4294967295 0 14 16 16 output-steinberger_projected_density/solution/solution-00000 5.81342027e-01 1.11483000e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01 4.85706828e+05 -4.25073533e+06 0.00000000e+00 0.00000000e+00
1 1.000000000000e+05 1.000000000000e+05 192 1891 833 2499 13 11 4294967295 11 12 14 14 output-steinberger_projected_density/solution/solution-00001 5.78335271e-01 1.11468431e+00 2.73000000e+02 2.43540495e+03 4.25000000e+03 5.43727671e-01 -2.81356102e+06 -1.21533631e+06 0.00000000e+00 0.00000000e+00
0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499 0 0 0 4294967295 14 16 16 output-steinberger_projected_density/solution/solution-00000 5.81342027e-01 1.11483000e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01 4.85706828e+05 -4.25073533e+06 0.00000000e+00 0.00000000e+00
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know this existed in the code prior to this PR, but why is the code returning the max 32-bit unsigned integer for "composition solver 2"?

@@ -24,5 +24,5 @@
# 24: Outward heat flux through boundary with indicator 1 ("top") (W)
# 25: Outward heat flux through boundary with indicator 2 ("left") (W)
# 26: Outward heat flux through boundary with indicator 3 ("right") (W)
0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499 5 0 0 4294967291 0 54 64 64 output-steinberger_projected_density_dcPicard/solution/solution-00000 5.78314979e-01 1.09945267e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01 4.99327014e+05 -4.27418614e+06 0.00000000e+00 0.00000000e+00
1 1.000000000000e+05 1.000000000000e+05 192 1891 833 2499 5 49 41 4294967291 40 59 69 69 output-steinberger_projected_density_dcPicard/solution/solution-00001 5.77687184e-01 1.11532571e+00 2.73000000e+02 2.43583623e+03 4.25000000e+03 5.43836114e-01 -2.71821909e+06 -1.13647007e+06 0.00000000e+00 0.00000000e+00
0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499 5 0 0 0 4294967291 54 64 64 output-steinberger_projected_density_dcPicard/solution/solution-00000 5.78314979e-01 1.09945267e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01 4.99327014e+05 -4.27418614e+06 0.00000000e+00 0.00000000e+00
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not quite the max 32-bit unsigned integer

if (adv_field.is_temperature() ||
adv_field.compositional_variable != introspection.find_composition_type(CompositionalFieldDescription::density))
current_constraints.distribute (distributed_vector);
if (adv_field.is_temperature() ||
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add a comment to explain the logic behind this if statement?

@tjhei
Copy link
Member

tjhei commented Feb 14, 2025

Thank you for cleaning up the code I wrote. 😉

@gassmoeller
Copy link
Member Author

Thanks for the reviews! I have addressed the comments as far as they are connected to this PR. I will open a separate PR to adress the two issues with global statistics.

I know this existed in the code prior to this PR, but why is the code returning the max 32-bit unsigned integer for "composition solver 2"?

This is caused by behavior in the global_statistics postprocessor. If we do not solve an advection field (e.g. because we prescribe it) it reports the number of iterations as numbers::invalid_unsigned_int which is std::numeric_limits<unsigned int>::max() - 1 (=4294967295). I am also not sure why the one test has a slightly different number (4294967291). I will address this in a separate PR.

What is the reason for this change in ordering?

This is also caused by behavior in the global_statistics postprocessor. It reports the statistics columns in the order the advection fields have been solved. Because we now group all prescribed fields and interpolate them in a combined operation after all other fields this order has changed. I am actually not convinced that reporting the columns in the order of solving the fields makes much sense, we should just report temperature first, then all other fields in order. So far the two options were the same, but now they arent any more.

If you are ok with me addressing the issue with global_statistics in a separate PR this should be good to go.

@gassmoeller gassmoeller force-pushed the unify_prescribing_fields branch from 1e274eb to 813e955 Compare February 28, 2025 16:13
@bobmyhill bobmyhill merged commit be825e7 into geodynamics:main Mar 2, 2025
8 checks passed
@gassmoeller gassmoeller deleted the unify_prescribing_fields branch March 3, 2025 10:19
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants