diff --git a/source/material_model/latent_heat.cc b/source/material_model/latent_heat.cc index b6f1f340d11..69c420c5570 100644 --- a/source/material_model/latent_heat.cc +++ b/source/material_model/latent_heat.cc @@ -73,8 +73,8 @@ namespace aspect if ((composition_viscosity_prefactor != 1.0) && (composition.size() > 0)) { // geometric interpolation - out.viscosities[i] = (std::pow(10, ((1-composition[0]) * std::log10(eta*visc_temperature_dependence) - + composition[0] * std::log10(eta*composition_viscosity_prefactor*visc_temperature_dependence)))); + out.viscosities[i] = (std::pow(10., ((1-composition[0]) * std::log10(eta*visc_temperature_dependence) + + composition[0] * std::log10(eta*composition_viscosity_prefactor*visc_temperature_dependence)))); } else out.viscosities[i] = visc_composition_dependence * visc_temperature_dependence * eta; @@ -130,10 +130,8 @@ namespace aspect const double phaseFunction = phase_function.compute_value(phase_in); - // Note that for the densities we have a list of jumps, so the index used - // in the loop corresponds to the index of the phase transition. For the - // viscosities we have a list of prefactors, which has one more entry - // for the first layer, so we have to use phase+1 as the index. + // For the densities we have a list of jumps, so the index used + // in the loop corresponds to the index of the phase transition. if (composition.size()==0) { phase_dependence += phaseFunction * density_jumps[transition_index]; @@ -146,14 +144,16 @@ namespace aspect phase_dependence += phaseFunction * density_jumps[transition_index] * composition[0]; } - // Viscosities + // For the viscosities we have a list of prefactors, which has one more entry + // for the first layer for each composition. Consequently, we have to use + // transition_index + 1 (+1 for additional compositions) as the index. if (transition_index+1 < phase_function.n_phases_for_each_composition()[0]) // 1st compositional field - viscosity_phase_prefactors[0] = std::pow(10, + viscosity_phase_prefactors[0] = std::pow(10., phaseFunction * std::log10(viscosity_phase_prefactors[0] * phase_prefactors[transition_index+1]) + (1. - phaseFunction) * std::log10(viscosity_phase_prefactors[0])); else if (transition_index+2 < phase_function.n_phases_for_each_composition()[0] + phase_function.n_phases_for_each_composition()[1]) // 2nd compositional field - viscosity_phase_prefactors[1] = std::pow(10, + viscosity_phase_prefactors[1] = std::pow(10., phaseFunction * composition[0] * std::log10(viscosity_phase_prefactors[1] * phase_prefactors[transition_index+2]) + (1. - phaseFunction) * std::log10(viscosity_phase_prefactors[1])); @@ -181,8 +181,8 @@ namespace aspect if (composition.size()==0) out.viscosities[i] = out.viscosities[i]*viscosity_phase_prefactors[0]; else - out.viscosities[i] = (std::pow(10, ((1-composition[0]) * std::log10(out.viscosities[i]*viscosity_phase_prefactors[0]) - + composition[0] * std::log10(out.viscosities[i]*viscosity_phase_prefactors[1])))); + out.viscosities[i] = (std::pow(10., ((1-composition[0]) * std::log10(out.viscosities[i]*viscosity_phase_prefactors[0]) + + composition[0] * std::log10(out.viscosities[i]*viscosity_phase_prefactors[1])))); out.viscosities[i] = std::max(minimum_viscosity, std::min(maximum_viscosity, out.viscosities[i])); } @@ -323,11 +323,12 @@ namespace aspect "Units: \\si{\\pascal\\per\\kelvin}."); prm.declare_entry ("Viscosity prefactors", "all:1", Patterns::Anything(), - "A list of prefactors for the viscosity for each phase. The reference " - "viscosity (modified by any compositional prefactors) will be multiplied " - "by this factor to get the corresponding viscosity for each phase. " + "A list of prefactors for the viscosity for each phase. The ``Viscosity'' " + "parameter (modified by the ``Composition viscosity prefactor'', depending on " + "composition) will be multiplied by this factor to get the corresponding " + "viscosity for each phase. " "List must have the same number of entries as there are phases, that is one " - "more than Phase transition depths for each composition that is used in the " + "more than ``Phase transition depths'' for each composition that is used in the " "model. " "Units: non-dimensional."); prm.declare_entry ("Minimum viscosity", "1e19", diff --git a/tests/latent_heat_viscosity/screen-output b/tests/latent_heat_viscosity/screen-output index 790e1021e37..1f93e93a380 100644 --- a/tests/latent_heat_viscosity/screen-output +++ b/tests/latent_heat_viscosity/screen-output @@ -24,7 +24,7 @@ Number of degrees of freedom: 10,656 (6,528+864+3,264) *** Timestep 0: t=0 years, dt=0 years Solving temperature system... 0 iterations. - Solving Stokes system (GMG)... 141+0 iterations. + Solving Stokes system (GMG)... 142+0 iterations. Number of active cells: 2,535 (on 5 levels) Number of degrees of freedom: 35,586 (21,840+2,826+10,920) @@ -37,18 +37,18 @@ Number of degrees of freedom: 35,586 (21,840+2,826+10,920) Writing graphical output: output-latent_heat_viscosity/solution/solution-00000 RMS, max velocity: 0.0181 m/year, 0.077 m/year Temperature min/avg/max: 273 K, 1587 K, 2600 K - Heat fluxes through boundary parts: -6.56e+05 W, 9.45e+05 W + Heat fluxes through boundary parts: -6.56e+05 W, 9.512e+05 W Writing depth average: output-latent_heat_viscosity/depth_average *** Timestep 1: t=1.15653e+06 years, dt=1.15653e+06 years Solving temperature system... 21 iterations. - Solving Stokes system (GMG)... 140+0 iterations. + Solving Stokes system (GMG)... 139+0 iterations. Postprocessing: Writing graphical output: output-latent_heat_viscosity/solution/solution-00001 RMS, max velocity: 0.018 m/year, 0.077 m/year Temperature min/avg/max: 273 K, 1587 K, 2600 K - Heat fluxes through boundary parts: -7.877e+05 W, 1.459e+06 W + Heat fluxes through boundary parts: -7.877e+05 W, 1.464e+06 W Writing depth average: output-latent_heat_viscosity/depth_average Termination requested by criterion: end step