Skip to content

Use current-surface depth for adiabatic pressure in visco-plastic rheology #2

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
Minerallo wants to merge 2 commits into
base: current-surface-depth-infrastructure
Choose a base branch
from
Open

Use current-surface depth for adiabatic pressure in visco-plastic rheology #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
14481a1
Use current-surface depth for adiabatic visco-plastic rheology
Minerallo Jul 26, 2026
99489b1
Add changelog for current-surface adiabatic pressure
Minerallo Jul 28, 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
6 changes: 6 additions & 0 deletions doc/modules/changes/20260728_Minerallo
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
Fixed: When mesh deformation is enabled, the adiabatic pressure used by the
visco-plastic material model for plastic yielding and viscous creep is now
evaluated using depth below the current deformed surface. This prevents surface
depressions such as trenches from becoming artificially strong.
<br>
(Michaël Pons, 2026/07/28)
4 changes: 2 additions & 2 deletions doc/sphinx/parameters/Material_20model.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5547,7 +5547,7 @@ If a compositional field named &rsquo;noninitial\_plastic\_strain&rsquo; is incl

**Pattern:** [Bool]

**Documentation:** Whether to use the adiabatic pressure instead of the full pressure (default) when calculating viscous creep. This may be helpful in models where the full pressure has an unusually large negative value arising from large negative dynamic pressure, resulting in solver convergence issue and in some cases a viscosity of zero.
**Documentation:** Whether to use the adiabatic pressure instead of the full pressure (default) when calculating viscous creep. This may be helpful in models where the full pressure has an unusually large negative value arising from large negative dynamic pressure, resulting in solver convergence issue and in some cases a viscosity of zero. With mesh deformation, the adiabatic profile is evaluated using depth below the current deformed surface.
::::

::::{dropdown} __Parameter:__ {ref}`Use adiabatic pressure in plasticity<parameters:Material_20model/Visco_20Plastic/Use_20adiabatic_20pressure_20in_20plasticity>`
Expand All @@ -5556,7 +5556,7 @@ If a compositional field named &rsquo;noninitial\_plastic\_strain&rsquo; is incl

**Pattern:** [Bool]

**Documentation:** Whether to use the adiabatic pressure instead of the full pressure when calculating plastic yield stress. This may be helpful in models where the full pressure has unusually large variations, resulting in solver convergence issues. Be aware that this setting will change the plastic shear band angle.
**Documentation:** Whether to use the adiabatic pressure instead of the full pressure when calculating plastic yield stress. This may be helpful in models where the full pressure has unusually large variations, resulting in solver convergence issues. With mesh deformation, the adiabatic profile is evaluated using depth below the current deformed surface, so surface depressions such as oceanic trenches remain at the adiabatic surface pressure. Be aware that this setting will also change the plastic shear band angle.
::::

::::{dropdown} __Parameter:__ {ref}`Use dominant phase for viscosity<parameters:Material_20model/Visco_20Plastic/Use_20dominant_20phase_20for_20viscosity>`
Expand Down
8 changes: 8 additions & 0 deletions include/aspect/adiabatic_conditions/interface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,14 @@ namespace aspect
virtual
double pressure (const Point<dim> &p) const = 0;

/**
* Return the adiabatic pressure at the representative point
* corresponding to @p depth. The depth must lie between zero and the
* maximal depth of the geometry model.
*/
double
pressure_at_depth (const double depth) const;

/**
* Return the reference_density at a given point of the domain.
*/
Expand Down
14 changes: 14 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 @@ -204,11 +204,23 @@ namespace aspect
IsostrainViscosities
calculate_isostrain_viscosities ( const MaterialModel::MaterialModelInputs<dim> &in,
const unsigned int i,
const double current_surface_adiabatic_pressure,
const std::vector<double> &volume_fractions,
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;

/**
* Return adiabatic pressures evaluated using depth below the current
* surface at all evaluation points. These values are used by the
* pressure-dependent parts of the rheology for which adiabatic
* pressure is enabled. A NaN entry indicates that the existing
* pressure selection in the rheology should be used.
*/
std::vector<double>
compute_current_surface_adiabatic_pressures(
const MaterialModel::MaterialModelInputs<dim> &in) const;

/**
* A function that fills the viscosity derivatives in the
* MaterialModelOutputs object that is handed over, if they exist.
Expand All @@ -221,6 +233,7 @@ namespace aspect
void compute_viscosity_derivatives(const unsigned int point_index,
const std::vector<double> &volume_fractions,
const IsostrainViscosities &isostrain_values,
const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs<dim> &in,
MaterialModel::MaterialModelOutputs<dim> &out,
const std::vector<double> &phase_function_values = std::vector<double>(),
Expand Down Expand Up @@ -268,6 +281,7 @@ namespace aspect
void fill_plastic_outputs(const unsigned int point_index,
const std::vector<double> &volume_fractions,
const bool plastic_yielding,
const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs<dim> &in,
MaterialModel::MaterialModelOutputs<dim> &out,
const IsostrainViscosities &isostrain_viscosities) const;
Expand Down
9 changes: 9 additions & 0 deletions source/adiabatic_conditions/interface.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,15 @@ namespace aspect
{
namespace AdiabaticConditions
{
template <int dim>
double
Interface<dim>::pressure_at_depth(const double depth) const
{
return pressure(this->get_geometry_model().representative_point(depth));
}



template <int dim>
void Interface<dim>::get_adiabatic_temperature_profile(std::vector<double> &values) const
{
Expand Down
82 changes: 72 additions & 10 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 @@ -19,6 +19,7 @@
*/

#include <aspect/material_model/rheology/visco_plastic.h>
#include <aspect/mesh_deformation/interface.h>
#include <aspect/material_model/utilities.h>
#include <aspect/utilities.h>
#include <aspect/newton.h>
Expand Down Expand Up @@ -158,6 +159,7 @@ namespace aspect
ViscoPlastic<dim>::
calculate_isostrain_viscosities (const MaterialModel::MaterialModelInputs<dim> &in,
const unsigned int i,
const double current_surface_adiabatic_pressure,
const std::vector<double> &volume_fractions,
const std::vector<double> &phase_function_values,
const std::vector<unsigned int> &n_phase_transitions_per_composition) const
Expand Down Expand Up @@ -250,7 +252,12 @@ namespace aspect
double pressure_for_creep = in.pressure[i];

if (use_adiabatic_pressure_in_creep)
pressure_for_creep = this->get_adiabatic_conditions().pressure(in.position[i]);
pressure_for_creep =
(std::isnan(current_surface_adiabatic_pressure)
?
this->get_adiabatic_conditions().pressure(in.position[i])
:
current_surface_adiabatic_pressure);

const double viscosity_diffusion
= (viscous_flow_law != dislocation
Expand Down Expand Up @@ -429,9 +436,13 @@ namespace aspect
// than the lithostatic pressure.

double pressure_for_plasticity = in.pressure[i];

if (use_adiabatic_pressure_in_plasticity)
pressure_for_plasticity = this->get_adiabatic_conditions().pressure(in.position[i]);
pressure_for_plasticity =
(std::isnan(current_surface_adiabatic_pressure)
?
this->get_adiabatic_conditions().pressure(in.position[i])
:
current_surface_adiabatic_pressure);

if (allow_negative_pressures_in_plasticity == false)
pressure_for_plasticity = std::max(pressure_for_plasticity,0.0);
Expand Down Expand Up @@ -514,12 +525,48 @@ namespace aspect



template <int dim>
std::vector<double>
ViscoPlastic<dim>::compute_current_surface_adiabatic_pressures(
const MaterialModel::MaterialModelInputs<dim> &in) const
{
std::vector<double> pressures(
in.n_evaluation_points(),
std::numeric_limits<double>::quiet_NaN());

if ((use_adiabatic_pressure_in_plasticity == false
&& use_adiabatic_pressure_in_creep == false)
|| this->get_parameters().mesh_deformation_enabled == false
|| in.current_cell.state() != IteratorState::valid)
return pressures;

const std::vector<double> depths =
this->get_mesh_deformation_handler().depth_below_current_surface(
in.current_cell, in.position);
const double maximal_depth = this->get_geometry_model().maximal_depth();

for (unsigned int q = 0; q < in.n_evaluation_points(); ++q)
{
// A deformed domain can extend beyond the reference geometry.
// representative_point() is only defined inside its depth range.
const double bounded_depth =
std::min(std::max(depths[q], 0.0), maximal_depth);
pressures[q] =
this->get_adiabatic_conditions().pressure_at_depth(bounded_depth);
}

return pressures;
}



template <int dim>
void
ViscoPlastic<dim>::
compute_viscosity_derivatives(const unsigned int i,
const std::vector<double> &volume_fractions,
const IsostrainViscosities &current_isostrain_values,
const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs<dim> &in,
MaterialModel::MaterialModelOutputs<dim> &out,
const std::vector<double> &phase_function_values,
Expand Down Expand Up @@ -565,7 +612,9 @@ namespace aspect
in_derivatives.strain_rate[i] = forward_strain_rate;

const IsostrainViscosities forward_isostrain_values =
calculate_isostrain_viscosities(in_derivatives, i, volume_fractions,
calculate_isostrain_viscosities(in_derivatives, i,
current_surface_adiabatic_pressure,
volume_fractions,
phase_function_values, n_phase_transitions_per_composition);

// For each composition of the independent component, compute the derivative.
Expand Down Expand Up @@ -603,7 +652,9 @@ namespace aspect
in_derivatives.strain_rate[i] = current_strain_rate;

const IsostrainViscosities forward_isostrain_values =
calculate_isostrain_viscosities(in_derivatives, i, volume_fractions,
calculate_isostrain_viscosities(in_derivatives, i,
current_surface_adiabatic_pressure,
volume_fractions,
phase_function_values, n_phase_transitions_per_composition);

for (unsigned int composition_index = 0; composition_index < n_compositions; ++composition_index)
Expand Down Expand Up @@ -769,15 +820,21 @@ namespace aspect
"This may be helpful in models where the "
"full pressure has an unusually large negative value arising from "
"large negative dynamic pressure, resulting in solver convergence "
"issue and in some cases a viscosity of zero.");
"issue and in some cases a viscosity of zero. With mesh deformation, "
"the adiabatic profile is evaluated using depth below the current "
"deformed surface.");
prm.declare_entry ("Use adiabatic pressure in plasticity", "false",
Patterns::Bool (),
"Whether to use the adiabatic pressure instead of the full "
"pressure when calculating plastic yield stress. "
"This may be helpful in models where the "
"full pressure has unusually large variations, resulting "
"in solver convergence issues. Be aware that this setting "
"will change the plastic shear band angle.");
"in solver convergence issues. With mesh deformation, the "
"adiabatic profile is evaluated using depth below the current "
"deformed surface, so surface depressions such as oceanic "
"trenches remain at the adiabatic surface pressure. Be aware "
"that this setting will also change the "
"plastic shear band angle.");

// Diffusion creep parameters
Rheology::DiffusionCreep<dim>::declare_parameters(prm);
Expand Down Expand Up @@ -1014,6 +1071,7 @@ namespace aspect
fill_plastic_outputs(const unsigned int i,
const std::vector<double> &volume_fractions,
const bool plastic_yielding,
const double current_surface_adiabatic_pressure,
const MaterialModel::MaterialModelInputs<dim> &in,
MaterialModel::MaterialModelOutputs<dim> &out,
const IsostrainViscosities &isostrain_viscosities) const
Expand All @@ -1032,9 +1090,13 @@ namespace aspect
plastic_out->yielding[i] = plastic_yielding ? 1 : 0;

double pressure_for_plasticity = in.pressure[i];

if (use_adiabatic_pressure_in_plasticity)
pressure_for_plasticity = this->get_adiabatic_conditions().pressure(in.position[i]);
pressure_for_plasticity =
(std::isnan(current_surface_adiabatic_pressure)
?
this->get_adiabatic_conditions().pressure(in.position[i])
:
current_surface_adiabatic_pressure);

if (allow_negative_pressures_in_plasticity == false)
pressure_for_plasticity = std::max(pressure_for_plasticity, 0.0);
Expand Down
28 changes: 24 additions & 4 deletions source/material_model/visco_plastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,14 @@ namespace aspect
/* The following returns whether or not the material is plastically yielding
* as documented in evaluate.
*/
const IsostrainViscosities isostrain_viscosities = rheology->calculate_isostrain_viscosities(in, 0, volume_fractions, phase_function_values, phase_function.n_phase_transitions_for_each_composition());
const std::vector<double> current_surface_adiabatic_pressures =
rheology->compute_current_surface_adiabatic_pressures(in);
const IsostrainViscosities isostrain_viscosities =
rheology->calculate_isostrain_viscosities(in, 0,
current_surface_adiabatic_pressures[0],
volume_fractions,
phase_function_values,
phase_function.n_phase_transitions_for_each_composition());

std::vector<double>::const_iterator max_composition = std::max_element(volume_fractions.begin(), volume_fractions.end());
const bool plastic_yielding = isostrain_viscosities.composition_yielding[std::distance(volume_fractions.begin(), max_composition)];
Expand Down Expand Up @@ -116,6 +123,8 @@ namespace aspect
std::vector<double> phase_function_discrete_values = (use_dominant_phase_for_viscosity?
std::vector<double>(phase_function_discrete->n_phase_transitions(), 0.0): std::vector<double>());

const std::vector<double> current_surface_adiabatic_pressures =
rheology->compute_current_surface_adiabatic_pressures(in);

// Loop through all requested points
for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
Expand Down Expand Up @@ -214,12 +223,20 @@ namespace aspect
phase_function_discrete_values[j] = phase_function_discrete->compute_value(phase_inputs);
}
isostrain_viscosities =
rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_discrete_values, phase_function_discrete->n_phase_transitions_for_each_chemical_composition());
rheology->calculate_isostrain_viscosities(in, i,
current_surface_adiabatic_pressures[i],
volume_fractions,
phase_function_discrete_values,
phase_function_discrete->n_phase_transitions_for_each_chemical_composition());
}
else
{
isostrain_viscosities =
rheology->calculate_isostrain_viscosities(in, i, volume_fractions, phase_function_values, n_phase_transitions_for_each_chemical_composition);
rheology->calculate_isostrain_viscosities(in, i,
current_surface_adiabatic_pressures[i],
volume_fractions,
phase_function_values,
n_phase_transitions_for_each_chemical_composition);
}

// The isostrain condition implies that the viscosity averaging should be arithmetic (see above).
Expand All @@ -241,6 +258,7 @@ namespace aspect

rheology->compute_viscosity_derivatives(i, volume_fractions,
isostrain_viscosities,
current_surface_adiabatic_pressures[i],
in, out, phase_function_values,
n_phase_transitions_for_each_chemical_composition);
}
Expand Down Expand Up @@ -277,7 +295,9 @@ namespace aspect
// 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,
current_surface_adiabatic_pressures[i],
in, out, isostrain_viscosities);
rheology->fill_viscosity_outputs(i, volume_fractions, out, isostrain_viscosities);
}

Expand Down
Loading
Loading