Skip to content

Commit ffb71e1

Browse files
lhy11009lhy11009
committed
Rheology::CompositionalViscosityPrefactors: change the interface
1 parent 7613d51 commit ffb71e1

File tree

3 files changed

+57
-60
lines changed

3 files changed

+57
-60
lines changed

‎include/aspect/material_model/rheology/compositional_viscosity_prefactors.h‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -66,15 +66,19 @@ namespace aspect
6666
dislocation
6767
} modified_flow_laws;
6868

69+
6970
/**
7071
* Compute the viscosity.
7172
*/
7273
double
73-
compute_viscosity (const MaterialModel::MaterialModelInputs<dim> &in,
74-
const double base_viscosity,
75-
const unsigned int composition_index,
76-
const unsigned int q,
77-
const ModifiedFlowLaws &modified_flow_laws) const;
74+
compute_viscosity(
75+
const double temperature,
76+
const double pressure,
77+
const double bound_fluid_fraction,
78+
const double base_viscosity,
79+
const unsigned int composition_index,
80+
const ModifiedFlowLaws &modified_flow_laws) const;
81+
7882

7983
private:
8084
/**

‎source/material_model/rheology/compositional_viscosity_prefactors.cc‎

Lines changed: 13 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -41,13 +41,17 @@ namespace aspect
4141
CompositionalViscosityPrefactors<dim>::CompositionalViscosityPrefactors ()
4242
= default;
4343

44+
4445
template <int dim>
4546
double
46-
CompositionalViscosityPrefactors<dim>::compute_viscosity (const MaterialModel::MaterialModelInputs<dim> &in,
47-
const double base_viscosity,
48-
const unsigned int composition_index,
49-
const unsigned int q,
50-
const ModifiedFlowLaws &modified_flow_laws) const
47+
CompositionalViscosityPrefactors<dim>::compute_viscosity(
48+
const double temperature,
49+
const double pressure,
50+
const double bound_fluid_fraction,
51+
const double base_viscosity,
52+
const unsigned int composition_index,
53+
const ModifiedFlowLaws &modified_flow_laws) const
54+
5155
{
5256
double factored_viscosities = base_viscosity;
5357
switch (viscosity_prefactor_scheme)
@@ -62,28 +66,13 @@ namespace aspect
6266
// We calculate the atomic H/Si ppm (C_OH) at each point to compute the water fugacity of
6367
// olivine assuming a composition of 90 mol% Forsterite and 10 mol% Fayalite from Hirth
6468
// and Kohlstaedt 2004 10.1029/138GM06.
65-
const double temperature_for_fugacity = (this->simulator_is_past_initialization())
66-
?
67-
in.temperature[q]
68-
:
69-
this->get_adiabatic_conditions().temperature(in.position[q]);
70-
AssertThrow(temperature_for_fugacity != 0, ExcMessage(
71-
"The temperature used in the calculation for determining the water fugacity is zero. "
72-
"This is not allowed, because this value is used to divide through. It is probably "
73-
"being caused by the temperature being zero somewhere in the model. The relevant "
74-
"values for debugging are: temperature (" + Utilities::to_string(in.temperature[q]) +
75-
"), and adiabatic temperature ("
76-
+ Utilities::to_string(this->get_adiabatic_conditions().temperature(in.position[q])) +
77-
"). If the adiabatic temperature is 0, double check that you are correctly defining an "
78-
" 'Adiabatic conditions model'."));
79-
80-
const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
81-
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
69+
70+
const double mass_fraction_H2O = std::max(minimum_mass_fraction_water_for_dry_creep[composition_index], bound_fluid_fraction); // mass fraction of bound water
8271
const double mass_fraction_olivine = 1 - mass_fraction_H2O; // mass fraction of olivine
8372
const double COH = (mass_fraction_H2O/molar_mass_H2O) / (mass_fraction_olivine/molar_mass_olivine) * 1e6; // COH in H / Si ppm
8473
const double point_water_fugacity = COH / A_H2O *
85-
std::exp((activation_energy_H2O + in.pressure[q]*activation_volume_H2O)/
86-
(constants::gas_constant * temperature_for_fugacity));
74+
std::exp((activation_energy_H2O + pressure*activation_volume_H2O)/
75+
(constants::gas_constant * temperature));
8776
const double r = modified_flow_laws == diffusion
8877
?
8978
-diffusion_water_fugacity_exponents[composition_index]

‎source/material_model/rheology/visco_plastic.cc‎

Lines changed: 35 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -190,42 +190,54 @@ namespace aspect
190190
if (use_adiabatic_pressure_in_creep)
191191
pressure_for_creep = this->get_adiabatic_conditions().pressure(in.position[i]);
192192

193-
const double viscosity_diffusion
194-
= (viscous_flow_law != dislocation
193+
const double bound_fluid_fraction
194+
= (this->introspection().compositional_name_exists("bound_fluid")
195195
?
196-
diffusion_creep.compute_viscosity(pressure_for_creep, temperature_for_viscosity, j,
197-
phase_function_values,
198-
n_phase_transitions_per_composition)
196+
in.composition[i][this->introspection().compositional_index_for_name("bound_fluid")]
199197
:
200198
numbers::signaling_nan<double>());
201199

200+
double viscosity_diffusion = numbers::signaling_nan<double>();
201+
if (viscous_flow_law != dislocation)
202+
{
203+
const double viscosity_diffusion_base = diffusion_creep.compute_viscosity(pressure_for_creep, temperature_for_viscosity, j,
204+
phase_function_values,
205+
n_phase_transitions_per_composition);
206+
viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(temperature_for_viscosity, pressure_for_creep, bound_fluid_fraction, viscosity_diffusion_base, j,
207+
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
208+
209+
}
210+
202211
// Step 1b: compute viscosity from dislocation creep law
203-
const double viscosity_dislocation
204-
= (viscous_flow_law != diffusion
205-
?
206-
dislocation_creep.compute_viscosity(edot_ii,
207-
pressure_for_creep,
208-
temperature_for_viscosity,
209-
j,
210-
phase_function_values,
211-
n_phase_transitions_per_composition)
212-
:
213-
numbers::signaling_nan<double>());
212+
double viscosity_dislocation = numbers::signaling_nan<double>();
213+
if (viscous_flow_law != diffusion)
214+
{
215+
const double viscosity_dislocation_base = dislocation_creep.compute_viscosity(edot_ii,
216+
pressure_for_creep,
217+
temperature_for_viscosity,
218+
j,
219+
phase_function_values,
220+
n_phase_transitions_per_composition);
221+
222+
viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(temperature_for_viscosity, pressure_for_creep, bound_fluid_fraction, viscosity_dislocation_base, j,
223+
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
224+
225+
}
226+
214227

215228
// Step 1c: select which form of viscosity to use (diffusion, dislocation, their minimum or composite, or fk), and apply
216229
// pre-exponential weakening, if required.
217230
switch (viscous_flow_law)
218231
{
219232
case diffusion:
220233
{
221-
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \
222-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
234+
235+
non_yielding_viscosity = viscosity_diffusion;
223236
break;
224237
}
225238
case dislocation:
226239
{
227-
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \
228-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
240+
non_yielding_viscosity = viscosity_dislocation;
229241
break;
230242
}
231243
case frank_kamenetskii:
@@ -238,21 +250,13 @@ namespace aspect
238250
}
239251
case composite:
240252
{
241-
const double scaled_viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \
242-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
243-
const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \
244-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
245-
non_yielding_viscosity = (scaled_viscosity_diffusion * scaled_viscosity_dislocation)/
246-
(scaled_viscosity_diffusion + scaled_viscosity_dislocation);
253+
non_yielding_viscosity = (viscosity_diffusion * viscosity_dislocation)/
254+
(viscosity_diffusion + viscosity_dislocation);
247255
break;
248256
}
249257
case minimum_diffusion_dislocation:
250258
{
251-
const double scaled_viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \
252-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
253-
const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \
254-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
255-
non_yielding_viscosity = std::min(scaled_viscosity_diffusion, scaled_viscosity_dislocation);
259+
non_yielding_viscosity = std::min(viscosity_diffusion, viscosity_dislocation);
256260
break;
257261
}
258262
default:

0 commit comments

Comments
 (0)