diff --git a/CMakeLists.txt b/CMakeLists.txt index 5606709c6aa..47c7ee450a5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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.") diff --git a/include/aspect/material_model/rheology/compositional_viscosity_prefactors.h b/include/aspect/material_model/rheology/compositional_viscosity_prefactors.h index 7af8672c14a..7c68d30faca 100644 --- a/include/aspect/material_model/rheology/compositional_viscosity_prefactors.h +++ b/include/aspect/material_model/rheology/compositional_viscosity_prefactors.h @@ -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 { @@ -70,11 +70,14 @@ namespace aspect * Compute the viscosity. */ double - compute_viscosity (const MaterialModel::MaterialModelInputs &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: /** diff --git a/source/material_model/rheology/compositional_viscosity_prefactors.cc b/source/material_model/rheology/compositional_viscosity_prefactors.cc index 9fc84f47254..408fea3d4a1 100644 --- a/source/material_model/rheology/compositional_viscosity_prefactors.cc +++ b/source/material_model/rheology/compositional_viscosity_prefactors.cc @@ -43,11 +43,14 @@ namespace aspect template double - CompositionalViscosityPrefactors::compute_viscosity (const MaterialModel::MaterialModelInputs &in, - const double base_viscosity, - const unsigned int composition_index, - const unsigned int q, - const ModifiedFlowLaws &modified_flow_laws) const + CompositionalViscosityPrefactors::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) @@ -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] diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc index 94f91c241ea..b98b8b52f0d 100644 --- a/source/material_model/rheology/visco_plastic.cc +++ b/source/material_model/rheology/visco_plastic.cc @@ -190,27 +190,40 @@ 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 viscosity_diffusion = numbers::signaling_nan(); + 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::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 viscosity_dislocation = numbers::signaling_nan(); + 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::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. @@ -218,14 +231,13 @@ namespace aspect { case diffusion: { - non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \ - CompositionalViscosityPrefactors::ModifiedFlowLaws::diffusion); + + non_yielding_viscosity = viscosity_diffusion; break; } case dislocation: { - non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \ - CompositionalViscosityPrefactors::ModifiedFlowLaws::dislocation); + non_yielding_viscosity = viscosity_dislocation; break; } case frank_kamenetskii: @@ -238,21 +250,13 @@ namespace aspect } case composite: { - const double scaled_viscosity_diffusion = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i, \ - CompositionalViscosityPrefactors::ModifiedFlowLaws::diffusion); - const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \ - CompositionalViscosityPrefactors::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::ModifiedFlowLaws::diffusion); - const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i, \ - CompositionalViscosityPrefactors::ModifiedFlowLaws::dislocation); - non_yielding_viscosity = std::min(scaled_viscosity_diffusion, scaled_viscosity_dislocation); + non_yielding_viscosity = std::min(viscosity_diffusion, viscosity_dislocation); break; } default: