Only require old stresses for elasticity if necessary#7007
Conversation
|
Hi @gassmoeller , Thank you for working on this. I think this is a very useful improvement. I have run into a similar issue in my recent 3-D viscoelastic models. In my tests, requiring the old stress fields did not seem to provide a noticeable improvement in accuracy for the cases I examined, but it did have a large impact on computational cost. The additional stress fields significantly affect memory use, advection cost, matrix assembly, and overall solve time in large 3-D calculations. In my local setup, I have been using an option similar to use old stress = false. I have tested this with several benchmarks and model comparisons, and so far I have not found a significant difference in the results for the cases where the old unrotated stress is not required by the time-stepping scheme. I had been planning to discuss this after finishing my current project, so I am very glad to see this issue being addressed in this PR. This looks like a very helpful change for large 3-D viscoelastic models. |
anne-glerum
left a comment
anne-glerum
left a comment
There was a problem hiding this comment.
Hi @gassmoeller, this is going to be a great speed-up, thank you! I have a few tiny comments, and think that the logic in solver_schemes.cc needs to be changed.
|
|
||
| // Compute the total stress at time t. | ||
| const SymmetricTensor<2, dim> | ||
| // Compute the total stress at time $t+\Delta t_c$. |
There was a problem hiding this comment.
Since this PR doesn't change when the reactions are computed/applied, this is still computing the total stress at time t, at the beginning of the next timestep
| // Compute the total stress at time $t+\Delta t_c$. | |
| // Compute the total stress at time $t$. |
| const std::shared_ptr<MaterialModel::ReactionRateOutputs<dim>> reaction_rate_outputs | ||
| = material_outputs.template get_additional_output_object<MaterialModel::ReactionRateOutputs<dim>>(); | ||
|
|
There was a problem hiding this comment.
These lines are not needed in this PR, probably a leftover from anne-glerum#19.
| for (auto &c : stress_field_indices) | ||
| stress_initial_residual += system_rhs.block(introspection.block_indices.compositional_fields[c]).l2_norm() / n_stress_fields; |
There was a problem hiding this comment.
I think this should be restricted to the current fields ve_stress*, not the current and old fields. n_stress_fields could be 2*n_independent_components here, and we don't want to average the residual over both tensors.
|
|
||
| // Overwrite the initial residual | ||
| for (auto &c : stress_indices) | ||
| if (c<SymmetricTensor<2,dim>::n_independent_components) |
There was a problem hiding this comment.
The stresses don't need to be listed as the first fields, so their compositional field indexes aren't necessarily <n_independent_components. The index of each stress field within stress_indices would be for the first/current set of stresses.
At the moment, in some places decisions are based on the number of fields of type |
gassmoeller
force-pushed
the
only_use_as_many_fields_as_necessary
branch
from
a390921 to
85e1b98
Compare
anne-glerum
left a comment
anne-glerum
left a comment
There was a problem hiding this comment.
Hi Rene, thanks for going through all the prm files! The ones that have changed so far look good. I think some comments in the material model and one more thing in the solver_scheme still needs to be changed.
| // second set of n_independent_components fields, e.g. in 2D on field 3, 4 and 5. | ||
| // The old stress was advected into the current timestep, but not rotated. | ||
| // Below we update it to full stress of the current timestep after using it to | ||
| // compute the full stress. |
There was a problem hiding this comment.
These 3 lines were changed, but the old version was correct. In this PR the reaction rates are still used at the beginning of the current timestep to compute the full stress of the previous timestep.
| // Also update the second set of stresses, stress_old, with the newly computed stress, | ||
| // which in the next timestep will serve as the old stress advected but not rotated | ||
| // into the current timestep. For fields, this function fill_reaction_rates is only | ||
| // called at the end of the timestep, and so this update only happens once. |
There was a problem hiding this comment.
This comment has also been changed, but for the current PR the reaction rates are still used at the beginning of the timpestep to update the stress of the previous timestep.
| old_stress_indices.push_back(introspection.compositional_index_for_name("ve_stress_yz_old")); | ||
| } | ||
|
|
||
| for (unsigned int c=0; c<SymmetricTensor<2,dim>::n_independent_components; ++c) |
There was a problem hiding this comment.
As the stress fields do not have to be the first fields, looping over field 0 to n_independent_components is not necessarily correct. The loop will have to be over the first n_independent_components of stress_field_indices.
3 participants
@anne-glerum as discussed this is a rebased version of anne-glerum#20. Can you give this a proper review? I am missing a changelog entry, I will need to write that as well.
Summary of PR: Our elasticity implementation under certain conditions needs two sets of stress tensors, one for the current time step, and one for the old time step partially modified into the current time (advected with the flow, but not rotated according to it). However, this is only necessary under certain conditions (see: https://github.com/geodynamics/aspect/compare/main...gassmoeller:aspect:only_use_as_many_fields_as_necessary?expand=1#diff-06e3c0241f4cd5f26a2dd149ee3e25bf3217d99aafcff84d02a98fe32629db6bR806). If these conditions are not met, we can get away with one advected stress tensor, which significantly reduces the compute time and memory requirements.
We will also have to look into all the elasticity benchmarks and tests to see which ones need to be modified. And we should decide if it should be allowed to specify two stress tensors (current and old) even if the old one is not needed. We could allow it for backward compatibility (and to make it easier to switch between the two approaches). I havent actually thought about if my current implementation allows both or not.