Skip to content

Compositional viscosity prefactors interface #6849

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 2 commits into
base: main
Choose a base branch
from
+59 −62
Open

Compositional viscosity prefactors interface #6849

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
6c4363e
CMake: print Python info
tjhei Feb 2, 2026
b518d10
Rheology::CompositionalViscosityPrefactors: change the interface
lhy11009 Feb 3, 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
4 changes: 3 additions & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,9 @@ message(STATUS "Using ASPECT_WITH_PYTHON = '${ASPECT_WITH_PYTHON}'")
if(ASPECT_WITH_PYTHON)
find_package(Python3 COMPONENTS Interpreter Development NumPy)
if(Python3_FOUND)
message(STATUS " PYTHON_EXECUTABLE: ${Python3_EXECUTABLE}")
message(STATUS " Python3_EXECUTABLE: ${Python3_EXECUTABLE}")
message(STATUS " Python3_INCLUDE_DIRS: ${Python3_INCLUDE_DIRS}")
message(STATUS " Python3_NumPy_INCLUDE_DIRS: ${Python3_NumPy_INCLUDE_DIRS}")
list(APPEND ASPECT_INCLUDE_DIRS ${Python3_INCLUDE_DIRS} ${Python3_NumPy_INCLUDE_DIRS})
else()
message(FATAL_ERROR "Python not found. Consider disabling ASPECT_WITH_PYTHON.")
Expand Down
15 changes: 9 additions & 6 deletions include/aspect/material_model/rheology/compositional_viscosity_prefactors.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ namespace aspect
void
parse_parameters (ParameterHandler &prm);

// The flow laws that can be
// The flow law that can be
// currently modified.
enum ModifiedFlowLaws
{
Expand All @@ -70,11 +70,14 @@ namespace aspect
* Compute the viscosity.
*/
double
compute_viscosity (const MaterialModel::MaterialModelInputs<dim> &in,
const double base_viscosity,
const unsigned int composition_index,
const unsigned int q,
const ModifiedFlowLaws &modified_flow_laws) const;
compute_viscosity(
const double temperature,
const double pressure,
const double bound_fluid_fraction,
const double base_viscosity,
const unsigned int composition_index,
const ModifiedFlowLaws &modified_flow_laws) const;


private:
/**
Expand Down
36 changes: 12 additions & 24 deletions source/material_model/rheology/compositional_viscosity_prefactors.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,11 +43,14 @@ namespace aspect

template <int dim>
double
CompositionalViscosityPrefactors<dim>::compute_viscosity (const MaterialModel::MaterialModelInputs<dim> &in,
const double base_viscosity,
const unsigned int composition_index,
const unsigned int q,
const ModifiedFlowLaws &modified_flow_laws) const
CompositionalViscosityPrefactors<dim>::compute_viscosity(
const double temperature,
const double pressure,
const double bound_fluid_fraction,
const double base_viscosity,
const unsigned int composition_index,
const ModifiedFlowLaws &modified_flow_laws) const

{
double factored_viscosities = base_viscosity;
switch (viscosity_prefactor_scheme)
Expand All @@ -62,28 +65,13 @@ namespace aspect
// We calculate the atomic H/Si ppm (C_OH) at each point to compute the water fugacity of
// olivine assuming a composition of 90 mol% Forsterite and 10 mol% Fayalite from Hirth
// and Kohlstaedt 2004 10.1029/138GM06.
const double temperature_for_fugacity = (this->simulator_is_past_initialization())
?
in.temperature[q]
:
this->get_adiabatic_conditions().temperature(in.position[q]);
AssertThrow(temperature_for_fugacity != 0, ExcMessage(
"The temperature used in the calculation for determining the water fugacity is zero. "
"This is not allowed, because this value is used to divide through. It is probably "
"being caused by the temperature being zero somewhere in the model. The relevant "
"values for debugging are: temperature (" + Utilities::to_string(in.temperature[q]) +
"), and adiabatic temperature ("
+ Utilities::to_string(this->get_adiabatic_conditions().temperature(in.position[q])) +
"). If the adiabatic temperature is 0, double check that you are correctly defining an "
" 'Adiabatic conditions model'."));

const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
const double mass_fraction_H2O = std::max(minimum_mass_fraction_water_for_dry_creep[composition_index], in.composition[q][bound_fluid_idx]); // mass fraction of bound water

const double mass_fraction_H2O = std::max(minimum_mass_fraction_water_for_dry_creep[composition_index], bound_fluid_fraction); // mass fraction of bound water
const double mass_fraction_olivine = 1 - mass_fraction_H2O; // mass fraction of olivine
const double COH = (mass_fraction_H2O/molar_mass_H2O) / (mass_fraction_olivine/molar_mass_olivine) * 1e6; // COH in H / Si ppm
const double point_water_fugacity = COH / A_H2O *
std::exp((activation_energy_H2O + in.pressure[q]*activation_volume_H2O)/
(constants::gas_constant * temperature_for_fugacity));
std::exp((activation_energy_H2O + pressure*activation_volume_H2O)/
(constants::gas_constant * temperature));
const double r = modified_flow_laws == diffusion
?
-diffusion_water_fugacity_exponents[composition_index]
Expand Down
66 changes: 35 additions & 31 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 @@ -190,42 +190,54 @@ namespace aspect
if (use_adiabatic_pressure_in_creep)
pressure_for_creep = this->get_adiabatic_conditions().pressure(in.position[i]);

const double viscosity_diffusion
= (viscous_flow_law != dislocation
const double bound_fluid_fraction
= (this->introspection().compositional_name_exists("bound_fluid")
?
diffusion_creep.compute_viscosity(pressure_for_creep, temperature_for_viscosity, j,
phase_function_values,
n_phase_transitions_per_composition)
in.composition[i][this->introspection().compositional_index_for_name("bound_fluid")]
:
numbers::signaling_nan<double>());

double viscosity_diffusion = numbers::signaling_nan<double>();
if (viscous_flow_law != dislocation)
{
const double viscosity_diffusion_base = diffusion_creep.compute_viscosity(pressure_for_creep, temperature_for_viscosity, j,
phase_function_values,
n_phase_transitions_per_composition);
viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(temperature_for_viscosity, pressure_for_creep, bound_fluid_fraction, viscosity_diffusion_base, j,
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);

}

// Step 1b: compute viscosity from dislocation creep law
const double viscosity_dislocation
= (viscous_flow_law != diffusion
?
dislocation_creep.compute_viscosity(edot_ii,
pressure_for_creep,
temperature_for_viscosity,
j,
phase_function_values,
n_phase_transitions_per_composition)
:
numbers::signaling_nan<double>());
double viscosity_dislocation = numbers::signaling_nan<double>();
if (viscous_flow_law != diffusion)
{
const double viscosity_dislocation_base = dislocation_creep.compute_viscosity(edot_ii,
pressure_for_creep,
temperature_for_viscosity,
j,
phase_function_values,
n_phase_transitions_per_composition);

viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(temperature_for_viscosity, pressure_for_creep, bound_fluid_fraction, viscosity_dislocation_base, j,
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);

}


// Step 1c: select which form of viscosity to use (diffusion, dislocation, their minimum or composite, or fk), and apply
// pre-exponential weakening, if required.
switch (viscous_flow_law)
{
case diffusion:
{
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);

non_yielding_viscosity = viscosity_diffusion;
break;
}
case dislocation:
{
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
non_yielding_viscosity = viscosity_dislocation;
break;
}
case frank_kamenetskii:
Expand All @@ -238,21 +250,13 @@ namespace aspect
}
case composite:
{
const double scaled_viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
non_yielding_viscosity = (scaled_viscosity_diffusion * scaled_viscosity_dislocation)/
(scaled_viscosity_diffusion + scaled_viscosity_dislocation);
non_yielding_viscosity = (viscosity_diffusion * viscosity_dislocation)/
(viscosity_diffusion + viscosity_dislocation);
break;
}
case minimum_diffusion_dislocation:
{
const double scaled_viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
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);
non_yielding_viscosity = std::min(viscosity_diffusion, viscosity_dislocation);
break;
}
default:
Expand Down