Skip to content

Commit f3dc38f

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

File tree

3 files changed

+68
-46
lines changed

3 files changed

+68
-46
lines changed

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

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

69+
70+
/**
71+
* Return whether the current prefactor scheme requires viscosity scaling.
72+
* If this returns false, callers can skip invoking compute_viscosity.
73+
*/
74+
bool needs_scaling() const;
75+
76+
6977
/**
7078
* Compute the viscosity.
7179
*/
80+
// todo_new
7281
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;
82+
compute_viscosity(
83+
const double temperature,
84+
const double pressure,
85+
const double bound_fluid_fraction,
86+
const double base_viscosity,
87+
const unsigned int composition_index,
88+
const ModifiedFlowLaws &modified_flow_laws) const;
89+
7890

7991
private:
8092
/**

‎source/material_model/rheology/compositional_viscosity_prefactors.cc‎

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

44+
template <int dim>
45+
bool
46+
CompositionalViscosityPrefactors<dim>::needs_scaling() const
47+
{
48+
return (viscosity_prefactor_scheme != none);
49+
}
50+
4451
template <int dim>
4552
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
53+
CompositionalViscosityPrefactors<dim>::compute_viscosity(
54+
const double temperature,
55+
const double pressure,
56+
const double bound_fluid_fraction,
57+
const double base_viscosity,
58+
const unsigned int composition_index,
59+
const ModifiedFlowLaws &modified_flow_laws) const
60+
5161
{
5262
double factored_viscosities = base_viscosity;
5363
switch (viscosity_prefactor_scheme)
@@ -62,28 +72,13 @@ namespace aspect
6272
// We calculate the atomic H/Si ppm (C_OH) at each point to compute the water fugacity of
6373
// olivine assuming a composition of 90 mol% Forsterite and 10 mol% Fayalite from Hirth
6474
// 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
75+
76+
const double mass_fraction_H2O = std::max(minimum_mass_fraction_water_for_dry_creep[composition_index], bound_fluid_fraction); // mass fraction of bound water
8277
const double mass_fraction_olivine = 1 - mass_fraction_H2O; // mass fraction of olivine
8378
const double COH = (mass_fraction_H2O/molar_mass_H2O) / (mass_fraction_olivine/molar_mass_olivine) * 1e6; // COH in H / Si ppm
8479
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));
80+
std::exp((activation_energy_H2O + pressure*activation_volume_H2O)/
81+
(constants::gas_constant * temperature));
8782
const double r = modified_flow_laws == diffusion
8883
?
8984
-diffusion_water_fugacity_exponents[composition_index]

‎source/material_model/rheology/visco_plastic.cc‎

Lines changed: 32 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -190,7 +190,15 @@ 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
193+
const double bound_fluid_fraction
194+
= (compositional_viscosity_prefactors.needs_scaling()
195+
?
196+
in.composition[i][this->introspection().compositional_index_for_name("bound_fluid")]
197+
:
198+
numbers::signaling_nan<double>());
199+
200+
201+
const double viscosity_diffusion_base
194202
= (viscous_flow_law != dislocation
195203
?
196204
diffusion_creep.compute_viscosity(pressure_for_creep, temperature_for_viscosity, j,
@@ -199,8 +207,16 @@ namespace aspect
199207
:
200208
numbers::signaling_nan<double>());
201209

210+
const double viscosity_diffusion
211+
= (viscous_flow_law != dislocation
212+
?
213+
compositional_viscosity_prefactors.compute_viscosity(temperature_for_viscosity, pressure_for_creep, bound_fluid_fraction, viscosity_diffusion_base, j,
214+
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion)
215+
:
216+
numbers::signaling_nan<double>());
217+
202218
// Step 1b: compute viscosity from dislocation creep law
203-
const double viscosity_dislocation
219+
const double viscosity_dislocation_base
204220
= (viscous_flow_law != diffusion
205221
?
206222
dislocation_creep.compute_viscosity(edot_ii,
@@ -212,20 +228,27 @@ namespace aspect
212228
:
213229
numbers::signaling_nan<double>());
214230

231+
const double viscosity_dislocation
232+
= (viscous_flow_law != diffusion
233+
?
234+
compositional_viscosity_prefactors.compute_viscosity(temperature_for_viscosity, pressure_for_creep, bound_fluid_fraction, viscosity_dislocation_base, j,
235+
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation)
236+
:
237+
numbers::signaling_nan<double>());
238+
215239
// Step 1c: select which form of viscosity to use (diffusion, dislocation, their minimum or composite, or fk), and apply
216240
// pre-exponential weakening, if required.
217241
switch (viscous_flow_law)
218242
{
219243
case diffusion:
220244
{
221-
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \
222-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
245+
246+
non_yielding_viscosity = viscosity_diffusion;
223247
break;
224248
}
225249
case dislocation:
226250
{
227-
non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \
228-
CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
251+
non_yielding_viscosity = viscosity_dislocation;
229252
break;
230253
}
231254
case frank_kamenetskii:
@@ -238,21 +261,13 @@ namespace aspect
238261
}
239262
case composite:
240263
{
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);
264+
non_yielding_viscosity = (viscosity_diffusion * viscosity_dislocation)/
265+
(viscosity_diffusion + viscosity_dislocation);
247266
break;
248267
}
249268
case minimum_diffusion_dislocation:
250269
{
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);
270+
non_yielding_viscosity = std::min(viscosity_diffusion, viscosity_dislocation);
256271
break;
257272
}
258273
default:

0 commit comments

Comments
 (0)