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.

Sign up for GitHub

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

Jump to bottom

added isostress averaging to CompositeViscoPlastic #5990

Merged
merged 8 commits into from
Jul 31, 2024
Merged

added isostress averaging to CompositeViscoPlastic #5990

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
b1e55a0
added isostress averaging
bobmyhill Jul 25, 2024
4cc572a
added test
bobmyhill Jul 25, 2024
82c7607
added changelog
bobmyhill Jul 25, 2024
f41a86e
set averaging schemes in tests
bobmyhill Jul 25, 2024
791ea55
sanitize function order
bobmyhill Jul 26, 2024
84ab5db
improve starting guess
bobmyhill Jul 26, 2024
56a8c6a
fix comment
bobmyhill Jul 26, 2024
2f09e15
addressed comments
bobmyhill Jul 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions doc/modules/changes/20240725_myhill
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
New: The CompositeViscoPlastic rheology in ASPECT now
allows the user to select isostress (Reuss) or
isostrain (Voigt) averaging of viscosities for
multicomponent materials via the parameter
"Viscosity averaging scheme".
<br>
(Bob Myhill, 2024/07/25)
107 changes: 96 additions & 11 deletions include/aspect/material_model/rheology/composite_visco_plastic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,31 @@ namespace aspect

namespace Rheology
{
/**
* A namespace for selecting how to average the viscosity
* over multiple chemical compositions.
*
* Current options are:
* 'isostress': Reuss (harmonic) averaging of viscosities.
* 'isostrain': Voigt (arithmetic) averaging of viscosities.
*/
namespace ViscosityAveraging
{
enum Kind
{
isostress,
isostrain
};

/**
* Read the viscosity averaging scheme from the parameter file
* using the parameter name given in @p parameter_name and return
* the enum that corresponds to this operation.
*/
ViscosityAveraging::Kind
parse (const std::string &parameter_name,
const ParameterHandler &prm);
}

template <int dim>
class CompositeViscoPlastic : public ::aspect::SimulatorAccess<dim>
Expand Down Expand Up @@ -83,8 +108,56 @@ namespace aspect
const std::vector<double> &phase_function_values = std::vector<double>(),
const std::vector<unsigned int> &n_phase_transitions_per_composition = std::vector<unsigned int>()) const;

private:
/**
* Compute the isostress viscosity over all compositional fields
* based on the composite viscous creep law.
* If @p expected_n_phases_per_composition points to a vector of
* unsigned integers this is considered the number of phase transitions
* for each compositional field and viscosity will be first computed on
* each phase and then averaged for each compositional field.
*/
double
compute_isostress_viscosity (const double pressure,
bobmyhill marked this conversation as resolved.
Show resolved Hide resolved
const double temperature,
const double grain_size,
const std::vector<double> &volume_fractions,
const SymmetricTensor<2,dim> &strain_rate,
std::vector<double> &partial_strain_rates,
const std::vector<double> &phase_function_values = std::vector<double>(),
const std::vector<unsigned int> &n_phase_transitions_per_composition = std::vector<unsigned int>()) const;

/**
* Compute the total strain rate and the first derivative of log strain rate
* with respect to log viscoplastic stress from individual log strain rate components
* over all compositional fields. Also updates the partial_strain_rates vector.
*/
std::pair<double, double>
calculate_isostress_log_strain_rate_and_derivative(const std::vector<std::array<std::pair<double, double>, 4>> &logarithmic_strain_rates_and_stress_derivatives,
bobmyhill marked this conversation as resolved.
Show resolved Hide resolved
const double viscoplastic_stress,
std::vector<double> &partial_strain_rates) const;

/**
* Compute the compositional field viscosity
* Compute the isostrain viscosity over all compositional fields
* based on the composite viscous creep law.
* If @p expected_n_phases_per_composition points to a vector of
* unsigned integers this is considered the number of phase transitions
* for each compositional field and viscosity will be first computed on
* each phase and then averaged for each compositional field.
*/
double
compute_isostrain_viscosity (const double pressure,
const double temperature,
const double grain_size,
const std::vector<double> &volume_fractions,
const SymmetricTensor<2,dim> &strain_rate,
std::vector<double> &partial_strain_rates,
const std::vector<double> &phase_function_values = std::vector<double>(),
const std::vector<unsigned int> &n_phase_transitions_per_composition = std::vector<unsigned int>()) const;


/**
* Compute the viscosity for a single compositional field
* based on the composite viscous creep law.
* If @p n_phase_transitions_per_composition points to a vector of
* unsigned integers this is considered the number of phase transitions
Expand All @@ -103,29 +176,41 @@ namespace aspect

/**
* Compute the total strain rate and the first derivative of log strain rate
* with respect to log viscoplastic stress from individual log strain rate components.
* Also updates the partial_strain_rates vector
* with respect to log viscoplastic stress from individual log strain rate components
* for a single compositional field.
* Also updates the partial_strain_rates vector.
*/
std::pair<double, double>
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,
Copy link
Member

Choose a reason for hiding this comment

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

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 (...);

Copy link
Member Author

Choose a reason for hiding this comment

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

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?

Copy link
Member

Choose a reason for hiding this comment

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

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_....

Copy link
Member Author

Choose a reason for hiding this comment

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

Great, I've now made a new compute_isostrain_viscosity function. I agree that this makes more logical sense.

const double viscoplastic_stress,
std::vector<double> &partial_strain_rates) const;

private:
/**
* Enumeration for selecting which type of viscosity averaging to use.
*/
ViscosityAveraging::Kind viscosity_averaging_scheme;

/**
* Whether to use different deformation mechanisms
* Whether to use different deformation mechanisms.
*/
bool use_diffusion_creep;
bool use_dislocation_creep;
bool use_peierls_creep;
bool use_drucker_prager;

/**
* Vector storing which flow mechanisms are active
* Vector storing which flow mechanisms are active.
*/
std::vector<unsigned int> active_flow_mechanisms;

/**
* Total number of decomposed strain rates. This includes the maximum
* viscosity limiter but not the minimum viscosity limiter,
* which is arranged in parallel with the viscoplastic elements and
* therefore does not contribute to the total strain rate.
*/
static constexpr unsigned int n_decomposed_strain_rates = 5;

/**
* Pointers to objects for computing deformation mechanism
* strain rates and effective viscosities.
Expand All @@ -142,7 +227,7 @@ namespace aspect

/**
* The maximum viscosity, imposed via an isoviscous damper
* in series with the composite viscoplastic element
* in series with the composite viscoplastic element.
*/
double maximum_viscosity;

Expand Down Expand Up @@ -185,7 +270,7 @@ namespace aspect
unsigned int stress_max_iteration_number;

/**
* Useful number to aid in adding together exponentials
* Useful number to aid in adding together exponentials.
*/
const double logmin = std::log(std::numeric_limits<double>::min());
};
Expand Down
Loading