diff --git a/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.cpp b/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.cpp new file mode 100644 index 00000000000..c77f68c4010 --- /dev/null +++ b/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.cpp @@ -0,0 +1,53 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + * Created on October 9, 2024, 6:03 PM + */ + +#include "CheckMPLPhasesForSinglePhaseFlow.h" + +#include +#include +#include + +#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h" +#include "MaterialLib/MPL/Medium.h" +#include "MeshLib/Elements/Element.h" +#include "MeshLib/Mesh.h" + +namespace MaterialPropertyLib +{ +void checkMPLPhasesForSinglePhaseFlow( + MeshLib::Mesh const& mesh, + MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map) +{ + // Check all of the elements have a medium defined. + ranges::for_each(mesh.getElements() | MeshLib::views::ids, + [&](auto const& element_id) + { media_map.checkElementHasMedium(element_id); }); + + // Collect phases of all elements... + auto all_phases = + media_map.media() | + ranges::views::transform([&](auto const& medium) + { return &fluidPhase(*medium); }) | + ranges::to_vector; + + assert(!all_phases.empty()); + + // ... and check if any of the phases are different by name. + if (ranges::any_of(all_phases, + [p0 = all_phases.front()](auto const* const p) + { return p->name != p0->name; })) + { + OGS_FATAL( + "You are mixing liquid and gas phases in your model domain. OGS " + "does not yet know how to handle this."); + } +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h b/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h new file mode 100644 index 00000000000..e719ccb5226 --- /dev/null +++ b/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h @@ -0,0 +1,26 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + * Created on October 9, 2024, 6:03 PM + */ + +#pragma once + +namespace MeshLib +{ +class Mesh; +} + +namespace MaterialPropertyLib +{ +class MaterialSpatialDistributionMap; + +void checkMPLPhasesForSinglePhaseFlow( + MeshLib::Mesh const& mesh, + MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map); +} // namespace MaterialPropertyLib diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp index 5e55893938a..bf6823b7aca 100644 --- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp +++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp @@ -11,6 +11,7 @@ #include "CreateHydroMechanicsProcess.h" #include +#include #include #include "HydroMechanicsProcess.h" @@ -18,6 +19,7 @@ #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h" #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h" #include "MaterialLib/MPL/Medium.h" +#include "MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h" #include "MaterialLib/SolidModels/CreateConstitutiveRelation.h" #include "MaterialLib/SolidModels/MechanicsBase.h" #include "MeshLib/Utils/Is2DMeshOnRotatedVerticalPlane.h" @@ -215,36 +217,20 @@ std::unique_ptr createHydroMechanicsProcess( MaterialPropertyLib::density}; std::array const requiredSolidProperties = {MaterialPropertyLib::density}; - // setting default to suppress -Wmaybe-uninitialized warning - MaterialPropertyLib::Variable phase_pressure = - MaterialPropertyLib::Variable::liquid_phase_pressure; - for (auto const& element_id : mesh.getElements() | MeshLib::views::ids) + MaterialPropertyLib::checkMPLPhasesForSinglePhaseFlow(mesh, media_map); + + // The uniqueness of phase has already been checked in + // `checkMPLPhasesForSinglePhaseFlow`. + MaterialPropertyLib::Variable const phase_variable = + (*ranges::begin(media_map.media()))->hasPhase("Gas") + ? MaterialPropertyLib::Variable::gas_phase_pressure + : MaterialPropertyLib::Variable::liquid_phase_pressure; + for (auto const& medium : media_map.media()) { - media_map.checkElementHasMedium(element_id); - auto const& medium = *media_map.getMedium(element_id); - if (element_id == 0) - { - phase_pressure = - medium.hasPhase("Gas") - ? MaterialPropertyLib::Variable::gas_phase_pressure - : MaterialPropertyLib::Variable::liquid_phase_pressure; - } - else - { - auto const phase_pressure_tmp = - medium.hasPhase("Gas") - ? MaterialPropertyLib::Variable::gas_phase_pressure - : MaterialPropertyLib::Variable::liquid_phase_pressure; - if (phase_pressure != phase_pressure_tmp) - { - OGS_FATAL( - "You are mixing liquid and gas phases in your model domain." - "OGS does not yet know how to handle this."); - } - } - checkRequiredProperties(medium, requiredMediumProperties); - checkRequiredProperties(fluidPhase(medium), requiredFluidProperties); - checkRequiredProperties(medium.phase("Solid"), requiredSolidProperties); + checkRequiredProperties(*medium, requiredMediumProperties); + checkRequiredProperties(fluidPhase(*medium), requiredFluidProperties); + checkRequiredProperties(medium->phase("Solid"), + requiredSolidProperties); } DBUG("Media properties verified."); @@ -266,7 +252,7 @@ std::unique_ptr createHydroMechanicsProcess( hydraulic_process_id, mechanics_related_process_id, use_taylor_hood_elements, - phase_pressure}; + phase_variable}; SecondaryVariableCollection secondary_variables; diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h index 867359cb049..4bf6dacfc92 100644 --- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h +++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h @@ -182,8 +182,9 @@ void HydroMechanicsLocalAssemblerphase("Solid"); auto const& fluid = fluidPhase(*medium); - auto const& phase_pressure = _process_data.phase_pressure; MPL::VariableArray vars; + auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure + : vars.liquid_phase_pressure; auto const T_ref = medium->property(MPL::PropertyType::reference_temperature) @@ -219,9 +220,7 @@ void HydroMechanicsLocalAssembler(vars, x_position, t, dt); - auto const beta_p = fluid.property(MPL::PropertyType::density) - .template dValue(vars, phase_pressure, - x_position, t, dt) / - rho_fr; + auto const beta_p = + fluid.property(MPL::PropertyType::density) + .template dValue(vars, _process_data.phase_variable, + x_position, t, dt) / + rho_fr; // Set mechanical variables for the intrinsic permeability model // For stress dependent permeability. @@ -411,6 +411,8 @@ std::vector const& HydroMechanicsLocalAssembler< auto const& medium = _process_data.media_map.getMedium(_element.getID()); auto const& fluid = fluidPhase(*medium); MPL::VariableArray vars; + auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure + : vars.liquid_phase_pressure; // TODO (naumov) Temporary value not used by current material models. Need // extension of secondary variables interface. @@ -427,9 +429,7 @@ std::vector const& HydroMechanicsLocalAssembler< double const p_int_pt = _ip_data[ip].N_p.dot(p); - // setting both vars equal to enable MPL access for gas and liquid - // properties - vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt; + phase_pressure = p_int_pt; auto const alpha = medium->property(MPL::PropertyType::biot_coefficient) .template value(vars, x_position, t, dt); @@ -535,8 +535,9 @@ void HydroMechanicsLocalAssemblerhasPhase("Gas") ? vars.gas_phase_pressure + : vars.liquid_phase_pressure; auto const T_ref = medium->property(MPL::PropertyType::reference_temperature) @@ -563,9 +564,7 @@ void HydroMechanicsLocalAssembler(vars, x_position, t, dt); - auto const beta_p = fluid.property(MPL::PropertyType::density) - .template dValue(vars, phase_pressure, - x_position, t, dt) / - rho_fr; + auto const beta_p = + fluid.property(MPL::PropertyType::density) + .template dValue(vars, _process_data.phase_variable, + x_position, t, dt) / + rho_fr; auto const K_over_mu = K / mu; @@ -701,6 +701,8 @@ void HydroMechanicsLocalAssemblerphase("Solid"); auto const& fluid = fluidPhase(*medium); MPL::VariableArray vars; + auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure + : vars.liquid_phase_pressure; auto const T_ref = medium->property(MPL::PropertyType::reference_temperature) @@ -731,9 +733,7 @@ void HydroMechanicsLocalAssemblerproperty(MPL::PropertyType::biot_coefficient) .template value(vars, x_position, t, dt); @@ -1159,6 +1159,8 @@ void HydroMechanicsLocalAssemblerhasPhase("Gas") ? vars.gas_phase_pressure + : vars.liquid_phase_pressure; SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize); auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin( @@ -1177,9 +1179,7 @@ void HydroMechanicsLocalAssembler(vars, x_position, t, dt); double const p_int_pt = _ip_data[ip].N_p.dot(p); - // setting both vars equal to enable MPL access for gas and liquid - // properties - vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt; + phase_pressure = p_int_pt; // Set mechanical variables for the intrinsic permeability model // For stress dependent permeability. diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h index 5b530e0162c..86a532c6267 100644 --- a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h +++ b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h @@ -90,7 +90,7 @@ struct HydroMechanicsProcessData const bool use_taylor_hood_elements; - MaterialPropertyLib::Variable phase_pressure; + MaterialPropertyLib::Variable const phase_variable; MeshLib::PropertyVector* pressure_interpolated = nullptr; std::array*, 3> principal_stress_vector = { diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp index cca70eec3a1..31042261674 100644 --- a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp +++ b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp @@ -12,17 +12,13 @@ #include "CreateLiquidFlowProcess.h" #include -#include -#include -#include -#include #include #include "LiquidFlowProcess.h" #include "MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h" #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h" #include "MaterialLib/MPL/Properties/Constant.h" -#include "MaterialLib/PhysicalConstant.h" +#include "MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h" #include "MeshLib/Utils/GetElementRotationMatrices.h" #include "MeshLib/Utils/GetSpaceDimension.h" #include "ParameterLib/Utils.h" @@ -38,29 +34,7 @@ void checkMPLProperties( MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map, bool const is_equation_type_volume) { - // Check all of the elements have a medium defined. - ranges::for_each(mesh.getElements() | MeshLib::views::ids, - [&](auto const& element_id) - { media_map.checkElementHasMedium(element_id); }); - - // Collect phases of all elements... - auto all_phases = - media_map.media() | - ranges::views::transform([&](auto const& medium) - { return &fluidPhase(*medium); }) | - ranges::to_vector; - - assert(!all_phases.empty()); - - // ... and check if any of the phases are different by name. - if (ranges::any_of(all_phases, - [p0 = all_phases.front()](auto const* const p) - { return p->name != p0->name; })) - { - OGS_FATAL( - "You are mixing liquid and gas phases in your model domain. OGS " - "does not yet know how to handle this."); - } + MaterialPropertyLib::checkMPLPhasesForSinglePhaseFlow(mesh, media_map); std::array const required_medium_properties = { MaterialPropertyLib::reference_temperature, @@ -214,6 +188,13 @@ std::unique_ptr createLiquidFlowProcess( *aperture_config, "parameter", parameters, 1); } + // The uniqueness of phase has already been checked in + // `checkMPLProperties`. + MaterialPropertyLib::Variable const phase_variable = + (*ranges::begin(media_map.media()))->hasPhase("Gas") + ? MaterialPropertyLib::Variable::gas_phase_pressure + : MaterialPropertyLib::Variable::liquid_phase_pressure; + LiquidFlowData process_data{ covertEquationBalanceTypeFromString(equation_balance_type_str), std::move(media_map), @@ -223,7 +204,8 @@ std::unique_ptr createLiquidFlowProcess( std::move(specific_body_force), has_gravity, *aperture_size_parameter, - NumLib::ShapeMatrixCache{integration_order}}; + NumLib::ShapeMatrixCache{integration_order}, + phase_variable}; return std::make_unique( std::move(name), mesh, std::move(jacobian_assembler), parameters, diff --git a/ProcessLib/LiquidFlow/LiquidFlowData.h b/ProcessLib/LiquidFlow/LiquidFlowData.h index c0f5b8640c0..9d63e2faef1 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowData.h +++ b/ProcessLib/LiquidFlow/LiquidFlowData.h @@ -13,6 +13,7 @@ #include #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h" +#include "MaterialLib/MPL/VariableType.h" #include "NumLib/Fem/ShapeMatrixCache.h" #include "ParameterLib/Parameter.h" @@ -54,6 +55,8 @@ struct LiquidFlowData final /// caches for each mesh element type the shape matrix NumLib::ShapeMatrixCache shape_matrix_cache; + + MaterialPropertyLib::Variable const phase_variable; }; } // namespace LiquidFlow diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h index 6774e9ba19c..0bbf1b0aeca 100644 --- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h +++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h @@ -134,14 +134,10 @@ void LiquidFlowLocalAssembler:: auto const& medium = *_process_data.media_map.getMedium(_element.getID()); auto const& fluid_phase = fluidPhase(medium); - bool is_gas_phase = medium.hasPhase("Gas"); - auto const fluid_pressure_variable = - is_gas_phase ? MaterialPropertyLib::Variable::gas_phase_pressure - : MaterialPropertyLib::Variable::liquid_phase_pressure; MaterialPropertyLib::VariableArray vars; - auto& phase_pressue = - is_gas_phase ? vars.gas_phase_pressure : vars.liquid_phase_pressure; + auto& phase_pressure = medium.hasPhase("Gas") ? vars.gas_phase_pressure + : vars.liquid_phase_pressure; vars.temperature = medium[MaterialPropertyLib::PropertyType::reference_temperature] @@ -160,7 +156,7 @@ void LiquidFlowLocalAssembler:: auto const& ip_data = _ip_data[ip]; auto const& N = Ns[ip]; - phase_pressue = N.dot(local_p_vec); + phase_pressure = N.dot(local_p_vec); auto const [fluid_density, viscosity] = getFluidDensityAndViscosity(t, dt, pos, fluid_phase, vars); @@ -175,8 +171,8 @@ void LiquidFlowLocalAssembler:: auto const get_drho_dp = [&]() { return fluid_phase[MaterialPropertyLib::PropertyType::density] - .template dValue(vars, fluid_pressure_variable, pos, t, - dt); + .template dValue(vars, _process_data.phase_variable, + pos, t, dt); }; auto const storage = _process_data.equation_balance_type == EquationBalanceType::volume @@ -295,8 +291,8 @@ void LiquidFlowLocalAssembler:: auto const& fluid_phase = fluidPhase(medium); MaterialPropertyLib::VariableArray vars; - auto& phase_pressue = medium.hasPhase("Gas") ? vars.gas_phase_pressure - : vars.liquid_phase_pressure; + auto& phase_pressure = medium.hasPhase("Gas") ? vars.gas_phase_pressure + : vars.liquid_phase_pressure; vars.temperature = medium[MaterialPropertyLib::PropertyType::reference_temperature] @@ -315,7 +311,7 @@ void LiquidFlowLocalAssembler:: auto const& ip_data = _ip_data[ip]; auto const& N = Ns[ip]; - phase_pressue = N.dot(local_p_vec); + phase_pressure = N.dot(local_p_vec); auto const [fluid_density, viscosity] = getFluidDensityAndViscosity(t, dt, pos, fluid_phase, vars);