Skip to content

Add viscosity additional outputs #6974

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

Open
lhy11009 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Add viscosity additional outputs #6974

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e66d521
Add viscosity additional outputs
lhy11009 May 8, 2026
9795977
Update include/aspect/material_model/rheology/visco_plastic.h
lhy11009 May 22, 2026
2668740
Update include/aspect/material_model/rheology/visco_plastic.h
lhy11009 May 22, 2026
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
3 changes: 3 additions & 0 deletions doc/modules/changes/20260512_lhy11009
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
Added: additional named outputs for diffusion and dislocation viscosities in the visco-plastic material model, including rheology-specific handling of inactive flow laws and visualization support through the additional material output system.
<br>
(Haoyuan Li, 2026/05/12)
51 changes: 51 additions & 0 deletions include/aspect/material_model/rheology/visco_plastic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,32 @@ namespace aspect

};

/**
* Additional output fields for diffusion and dislocation viscosities.
*/
template <int dim>
class ViscosityAdditionalOutputs : public NamedAdditionalMaterialOutputs<dim>
{
public:
/**
* Constructor.
*/
ViscosityAdditionalOutputs(const unsigned int n_points);

std::vector<double>
get_nth_output(const unsigned int idx) const override;

/**
* Diffusion viscosities.
*/

@alarshi alarshi May 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Here, it would be helpful to add the details. As I understand, these are the diffusion viscosities that are computed in Rheology::DiffusionCreep::compute_viscosity() function before yielding. The values are only relevant when diffusion creep is present, i.e., viscous flow law is either diffusion or composite.
Or something along these lines.

std::vector<double> diffusion_viscosities;

/**
* Dislocation viscosities.
*/
Comment on lines +107 to +109

@alarshi alarshi May 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Same as above.

std::vector<double> dislocation_viscosities;
};

/**
* A data structure with the output of calculate_isostrain_viscosities.
*/
Expand Down Expand Up @@ -117,6 +143,16 @@ namespace aspect
* MaterialModel::PrescribedDilation::dilation_rhs_term.
*/
std::vector<double> dilation_rhs_terms;

/**
* Diffusion viscosities for each composition.
*/
std::vector<double> diffusion_viscosities;

/**
* Dislocation viscosities for each composition.
*/
std::vector<double> dislocation_viscosities;
};

namespace Rheology
Expand Down Expand Up @@ -209,6 +245,21 @@ namespace aspect
const MaterialModel::MaterialModelInputs<dim> &in,
MaterialModel::MaterialModelOutputs<dim> &out,
const IsostrainViscosities &isostrain_viscosities) const;
/**
* Create additional outputs for diffusion and dislocation viscosities.
*/
void
create_viscosity_outputs(MaterialModel::MaterialModelOutputs<dim> &out) const;


/**
* Fill additional outputs for diffusion and dislocation viscosities,
* if viscosity additional output object is created.
*/
Comment thread
lhy11009 marked this conversation as resolved.
void fill_viscosity_outputs(const unsigned int point_index,
const std::vector<double> &volume_fractions,
MaterialModel::MaterialModelOutputs<dim> &out,
const IsostrainViscosities &isostrain_viscosities) const;

/**
* Minimum strain rate used to stabilize the strain rate dependent rheology.
Expand Down
136 changes: 136 additions & 0 deletions source/material_model/rheology/visco_plastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,42 @@ namespace aspect



namespace
{
std::vector<std::string> make_viscosity_additional_outputs_names()
{
std::vector<std::string> names;
names.emplace_back("diffusion_viscosity");
names.emplace_back("dislocation_viscosity");
return names;
}
}



template <int dim>
ViscosityAdditionalOutputs<dim>::ViscosityAdditionalOutputs(const unsigned int n_points)
: NamedAdditionalMaterialOutputs<dim>(make_viscosity_additional_outputs_names()),
diffusion_viscosities(n_points, numbers::signaling_nan<double>()),
dislocation_viscosities(n_points, numbers::signaling_nan<double>())
{}



template <int dim>
std::vector<double>
ViscosityAdditionalOutputs<dim>::get_nth_output(const unsigned int idx) const
{
AssertIndexRange(idx, 2);

if (idx == 0)
return diffusion_viscosities;
else
return dislocation_viscosities;
}



namespace Rheology
{

Expand All @@ -110,6 +146,8 @@ namespace aspect
output_parameters.drucker_prager_parameters.resize(volume_fractions.size());
output_parameters.dilation_lhs_terms.resize(volume_fractions.size(), numbers::signaling_nan<double>());
output_parameters.dilation_rhs_terms.resize(volume_fractions.size(), numbers::signaling_nan<double>());
output_parameters.diffusion_viscosities.resize(volume_fractions.size(), numbers::signaling_nan<double>());

@lhy11009 lhy11009 May 12, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

However, in the calculation step, it seems the nan value won't trigger issue, as indicated by the examples of other variables.

output_parameters.dislocation_viscosities.resize(volume_fractions.size(), numbers::signaling_nan<double>());

// Assemble current and old stress tensor if elastic behavior is enabled
SymmetricTensor<2, dim> stress_0_advected = numbers::signaling_nan<SymmetricTensor<2, dim>>();
Expand Down Expand Up @@ -220,12 +258,14 @@ namespace aspect
{
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i,
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
output_parameters.diffusion_viscosities[j] = non_yielding_viscosity;
break;
}
case dislocation:
{
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i,
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
output_parameters.dislocation_viscosities[j] = non_yielding_viscosity;
break;
}
case frank_kamenetskii:
Expand All @@ -244,6 +284,9 @@ namespace aspect
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
non_yielding_viscosity = (scaled_viscosity_diffusion * scaled_viscosity_dislocation)/
(scaled_viscosity_diffusion + scaled_viscosity_dislocation);

output_parameters.diffusion_viscosities[j] = scaled_viscosity_diffusion;
output_parameters.dislocation_viscosities[j] = scaled_viscosity_dislocation;
break;
}
case minimum_diffusion_dislocation:
Expand All @@ -253,6 +296,9 @@ namespace aspect
const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i,
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
non_yielding_viscosity = std::min(scaled_viscosity_diffusion, scaled_viscosity_dislocation);

output_parameters.diffusion_viscosities[j] = scaled_viscosity_diffusion;
output_parameters.dislocation_viscosities[j] = scaled_viscosity_dislocation;
break;
}
default:
Expand Down Expand Up @@ -984,6 +1030,95 @@ namespace aspect
}
}

template <int dim>
void
Rheology::ViscoPlastic<dim>::create_viscosity_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
{
if (out.template has_additional_output_object<ViscosityAdditionalOutputs<dim>>() == false)
{
const unsigned int n_points = out.n_evaluation_points();

out.additional_outputs.push_back(
std::make_unique<ViscosityAdditionalOutputs<dim>>(n_points));
}
}

template <int dim>
void
Rheology::ViscoPlastic<dim>::fill_viscosity_outputs(
const unsigned int i,
const std::vector<double> &volume_fractions,
MaterialModel::MaterialModelOutputs<dim> &out,
const IsostrainViscosities &isostrain_viscosities) const
{
if (const std::shared_ptr<ViscosityAdditionalOutputs<dim>>
viscosity_out = out.template get_additional_output_object<ViscosityAdditionalOutputs<dim>>())
{
switch (viscous_flow_law)

@alarshi alarshi May 15, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Since you are doing almost the same operation in the each case, you could replace ln 1057-1114 by initializing the viscosities using a default value (std::max in your case) and a flag to decide whether to output diffusion viscosities or/and dislocation viscosities.
Something like
bool compute_dislocation_viscosity = (viscous_flow_law != dislocation && viscous_flow_law != frank_kamenetskii); and

if (compute_dislocation_viscosity)
     viscosity_out->dislocation_viscosities[i]
                    = MaterialUtilities::average_value(
                        volume_fractions,
                        isostrain_viscosities.dislocation_viscosities,
                        viscosity_averaging);

and similarly for the diffusion_viscosities output.

{
case diffusion:

@lhy11009 lhy11009 May 12, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Currently, I assign the value to max limit, if the viscosity is not computed by the viscous flow law (e.g. the dislocation creep would not be calculated if the flow law is diffusion). I'm not entirely sure this is the best approach. I have tested that assigning values to nan doesn't work and will trigger error in the post-processor.

The values will be written as "inf" in outputs.

@alarshi alarshi May 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Hmm, I see. Does the issue persist after #6973 was merged?
But in either case, I am not sure if having both dislocation and diffusion viscosity outputs for each viscous flow law is the right approach. See my comment above.

@lhy11009 lhy11009 May 29, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Oh, this and 6973 are separated.

{
viscosity_out->diffusion_viscosities[i]
= MaterialUtilities::average_value(
volume_fractions,
isostrain_viscosities.diffusion_viscosities,
viscosity_averaging);

viscosity_out->dislocation_viscosities[i]
= std::numeric_limits<double>::max();

break;
}

case dislocation:
{
viscosity_out->diffusion_viscosities[i]
= std::numeric_limits<double>::max();

viscosity_out->dislocation_viscosities[i]
= MaterialUtilities::average_value(
volume_fractions,
isostrain_viscosities.dislocation_viscosities,
viscosity_averaging);

break;
}

case composite:
case minimum_diffusion_dislocation:
{
viscosity_out->diffusion_viscosities[i]
= MaterialUtilities::average_value(
volume_fractions,
isostrain_viscosities.diffusion_viscosities,
viscosity_averaging);

viscosity_out->dislocation_viscosities[i]
= MaterialUtilities::average_value(
volume_fractions,
isostrain_viscosities.dislocation_viscosities,
viscosity_averaging);

break;
}

case frank_kamenetskii:
{
viscosity_out->diffusion_viscosities[i]
= std::numeric_limits<double>::max();

viscosity_out->dislocation_viscosities[i]
= std::numeric_limits<double>::max();

break;
}

default:
AssertThrow(false, ExcInternalError());
}
}
}

}
}
}
Expand All @@ -995,6 +1130,7 @@ namespace aspect
{
#define INSTANTIATE(dim) \
template class PlasticAdditionalOutputs<dim>; \
template class ViscosityAdditionalOutputs<dim>; \
\
namespace Rheology \
{ \
Expand Down
10 changes: 8 additions & 2 deletions source/material_model/visco_plastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -251,6 +251,8 @@ namespace aspect
isostrain_viscosities.composition_yielding.clear();
isostrain_viscosities.composition_viscosities.clear();
isostrain_viscosities.drucker_prager_parameters.clear();
isostrain_viscosities.diffusion_viscosities.clear();
isostrain_viscosities.dislocation_viscosities.clear();

out.viscosities[i] = numbers::signaling_nan<double>();

Expand All @@ -270,11 +272,14 @@ namespace aspect
// TODO only when requests_property is set to reaction_terms
rheology->strain_rheology.fill_reaction_outputs(in, i, rheology->min_strain_rate, plastic_yielding, out);

// Fill plastic outputs if they exist.
// Fill plastic outputs and additional viscosity outputs if they exist.
// The values in isostrain_viscosities only make sense when the calculate_isostrain_viscosities function
// has been called.
if (in.requests_property(MaterialProperties::additional_outputs))
rheology->fill_plastic_outputs(i, volume_fractions, plastic_yielding, in, out, isostrain_viscosities);
{
rheology->fill_plastic_outputs(i, volume_fractions, plastic_yielding, in, out, isostrain_viscosities);
rheology->fill_viscosity_outputs(i, volume_fractions, out, isostrain_viscosities);
}

if (this->get_parameters().enable_elasticity)
{
Expand Down Expand Up @@ -476,6 +481,7 @@ namespace aspect
ViscoPlastic<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
{
rheology->create_plastic_outputs(out);
rheology->create_viscosity_outputs(out);

if (this->get_parameters().enable_elasticity)
rheology->elastic_rheology.create_elastic_additional_outputs(out);
Expand Down
Loading
Loading