Skip to content

Commit

Permalink
address comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@jdannberg
jdannberg committed Jan 29, 2025
1 parent 525f2fe commit 2b2b703
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 19 deletions.
31 changes: 16 additions & 15 deletions source/material_model/latent_heat.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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];
Expand All @@ -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]));

Expand Down Expand Up @@ -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]));
}

Expand Down Expand Up @@ -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",
Expand Down
8 changes: 4 additions & 4 deletions tests/latent_heat_viscosity/screen-output
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down

0 comments on commit 2b2b703

Please sign in to comment.