added isostress averaging to CompositeViscoPlastic #5990
Conversation
bobmyhill
changed the title
|
/rebuild |
| calculate_log_strain_rate_and_derivative(const std::array<std::pair<double, double>, 4> &logarithmic_strain_rates_and_stress_derivatives, | ||
| const double viscoplastic_stress, | ||
| std::vector<double> &partial_strain_rates) const; | ||
| calculate_composition_log_strain_rate_and_derivative(const std::array<std::pair<double, double>, 4> &logarithmic_strain_rates_and_stress_derivatives, |
There was a problem hiding this comment.
actually, why is this function and one above (compute_composition_viscosity) public? they are also just called from inside the rheology model?
And can you explain the name change? Is it to align this function with the compute_composition_viscosity? From the description it is not obvious what connects this function to a composition. For symmetry reasons with calculate_isostress_log_strain_rate_and_derivative would it be correct to call this function calculate_isostrain_log_strain_rate_and_derivative (and similar for the compute_composition_viscosity function)?
From a structure perspective these are the function names I would find cleanest (but I will let you comment if this makes sense):
public:
double compute_viscosity (...);
private:
... compute_isostrain_viscosity (...);
... compute_isostress_viscosity (...);
... calculate_isostrain_log_strain_rate_and_derivative (...);
... calculate_isostress_log_strain_rate_and_derivative (...);
There was a problem hiding this comment.
I'd also like this symmetry, but I don't see a way to do it neatly.
The problem: while the isostress viscosity solves for a single stress, the isostrain rheology solves for one stress per compositional field. It felt natural to create inner functions that only do one solve, leading to the two functions for "composition" (used several times by the isostrain branch) and the two functions for "isostress" (only used once by the isostress branch).
I'm happy to group the loop over compositions into a compute_isostrain_viscosity if you think that would be neater, but that would still leave a calculate_composition_log_strain_rate_and_derivative.
What would you prefer?
There was a problem hiding this comment.
Thanks for the explanation, I understand now. Yes outsourcing the loop over compositions into a compute_isostrain_viscosity and then mentioning in the calculate_composition... functions that they are only computing things for one composition would be most intuitive I think. What confused me most about the current situation is that there was compute_isostress but no compute_isostrain, but instead compute_composition_....
There was a problem hiding this comment.
Great, I've now made a new compute_isostrain_viscosity function. I agree that this makes more logical sense.
| /** | ||
| * Enumeration for selecting which type of viscosity averaging to use. | ||
| */ | ||
| enum ViscosityAveragingScheme | ||
| { | ||
| isostrain, | ||
| isostress | ||
| } viscosity_averaging_scheme; |
There was a problem hiding this comment.
I think this is something that may be useful for other rheology models as well. Can you move it outside the class into its own namespace? Something like:
namespace ViscosityAveraging
{
enum Kind
{
isostrain,
isostress
};
}
Then you can use that as type for the private variable viscosity_averaging_scheme inside the rheology.
| // a certain (small) fraction. | ||
| if (volume_fractions[composition] > 2.*std::numeric_limits<double>::epsilon()) | ||
| { | ||
| std::vector<double> partial_strain_rates_composition(5, 0.); |
There was a problem hiding this comment.
Since you know at compile time the size of this vector it may as well be an array with 5 entries (composition viscosity will need a different interface for this).
And since you now use the magic number 5 in several places, it would be worth adding it as static constexpr n_deformation_mechanisms = 5; to the private part of the header file and use it instead of 5 in all places.
| if (use_diffusion_creep) | ||
| { | ||
| log_edot_and_deriv[i][0] = diffusion_creep->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, grain_size, diffusion_creep_parameters[i]); | ||
| log_edot_and_deriv[i][0].first += log_volume_fraction; |
There was a problem hiding this comment.
why do you add the log_volume_fraction to the log strain rate? is it because the addition in log-space is the same as the multiplication in without the log? please add a comment to explain.
| // between these components. | ||
|
|
||
| // The following while loop contains a Newton iteration to obtain the | ||
| // viscoplastic stress that is consistent with the total strain rate. |
There was a problem hiding this comment.
please add a comment saying that this also results in a correct partitioning between the flow mechanisms (the introductory sentences of the comment refer to the partitioning, but the last sentence doesnt explain that the Newton iteration will result in a correct partitioning).
| template <int dim> | ||
| void | ||
| CompositeViscoPlastic<dim>::declare_parameters (ParameterHandler &prm) | ||
| { | ||
| prm.declare_entry ("Viscosity averaging scheme", "isostress", | ||
| Patterns::Selection("isostress|isostrain|unspecified"), |
There was a problem hiding this comment.
no need to have an unspecified option if it is not working and not the default
|
@gassmoeller , brilliant, thanks for the comments. I've updated the code as suggested and modified the tests to use the new structure (the functions that they tested are now hidden). I haven't yet tested the use of phases, so I've removed them from the test cc files. I'll add more tests for phases in a later PR. This should be reasy for another review now. |
This PR implements isostress (Reuss) component averaging to the CompositeViscoPlastic rheology.
Aside from giving users the option of choosing from two averaging schemes, isostress averaging is the easiest way to implement a logically consistent viscoelastic model. This is because multiple elastic elements can be combined into a single element, and the "pseudoviscous" part of the elastic element (see Moresi et al., 2003) can be combined with the other viscous elements. Therefore, #5984 should build on this PR.