From 066b01f0c02872e1b6232612fd8d7ba7afebd5d2 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 28 May 2024 10:07:16 +0200 Subject: [PATCH 1/2] [PL/RM] Added empty constitutive setting code --- ProcessLib/RichardsMechanics/CMakeLists.txt | 1 + .../ConstitutiveRelations/Base.h | 81 +++++++++++++++++++ .../ConstitutiveRelations/ConstitutiveData.h | 38 +++++++++ .../ConstitutiveModels.h | 34 ++++++++ .../ConstitutiveSetting.cpp | 35 ++++++++ .../ConstitutiveSetting.h | 42 ++++++++++ .../CreateConstitutiveSetting.cpp | 32 ++++++++ .../CreateConstitutiveSetting.h | 33 ++++++++ .../RichardsMechanicsFEM-impl.h | 5 ++ 9 files changed, 301 insertions(+) create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveModels.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.cpp create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.h diff --git a/ProcessLib/RichardsMechanics/CMakeLists.txt b/ProcessLib/RichardsMechanics/CMakeLists.txt index 34d366a89f5..6b5f50fb71a 100644 --- a/ProcessLib/RichardsMechanics/CMakeLists.txt +++ b/ProcessLib/RichardsMechanics/CMakeLists.txt @@ -1,4 +1,5 @@ get_source_files(SOURCES) +append_source_files(SOURCES ConstitutiveRelations) ogs_add_library(RichardsMechanics ${SOURCES}) target_link_libraries(RichardsMechanics PUBLIC ProcessLib PRIVATE ParameterLib) diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h new file mode 100644 index 00000000000..be08c171856 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h @@ -0,0 +1,81 @@ +/** + * \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 + * + */ + +#pragma once + +#include "MaterialLib/MPL/Medium.h" +#include "MathLib/KelvinVector.h" +#include "ParameterLib/SpatialPosition.h" +#include "ProcessLib/ConstitutiveRelations/Base.h" +#include "ProcessLib/Reflection/ReflectionData.h" + +namespace ProcessLib::RichardsMechanics +{ + +using namespace ProcessLib::ConstitutiveRelations; +namespace KV = MathLib::KelvinVector; + +template +using KelvinVector = KV::KelvinVectorType; + +template +using KelvinMatrix = KV::KelvinMatrixType; + +template +using GlobalDimVector = Eigen::Vector; + +template +using GlobalDimMatrix = + Eigen::Matrix; + +/// Used to set a D dimensional vector to all not-a-number. +template +constexpr GlobalDimVector DVnan() +{ + return GlobalDimVector::Constant(nan); +} + +/// Used to set a D x D matrix to all not-a-number. +template +constexpr GlobalDimMatrix DMnan() +{ + return GlobalDimMatrix::Constant(nan); +} + +struct MediaData +{ + explicit MediaData(MaterialPropertyLib::Medium const& medium) + : medium{medium}, + liquid{medium.phase("AqueousLiquid")}, + solid{medium.phase("Solid")} + { + } + + MaterialPropertyLib::Medium const& medium; + MaterialPropertyLib::Phase const& liquid; + MaterialPropertyLib::Phase const& solid; +}; + +template +struct TemperatureData +{ + double T; + double T_prev; + Eigen::Vector grad_T; +}; + +template +struct CapillaryPressureData +{ + double p_cap; + double p_cap_prev; + Eigen::Vector grad_p_cap; +}; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h new file mode 100644 index 00000000000..8b750271bdb --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -0,0 +1,38 @@ +/** + * \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 + * + */ + +#pragma once + +#include "ProcessLib/ConstitutiveRelations/Base.h" + +namespace ProcessLib::RichardsMechanics +{ +// TODO directly declare these type aliases in Traits.h +/// Data whose state must be tracked by the TRM process. +template +using StatefulData = std::tuple<>; + +template +using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf< + StatefulData>; + +/// Data that is needed for output purposes, but not directly for the assembly. +template +using OutputData = std::tuple<>; + +/// Data that is needed for the equation system assembly. +template +using ConstitutiveData = std::tuple<>; + +/// Data that stores intermediate values, which are not needed outside the +/// constitutive setting. +template +using ConstitutiveTempData = std::tuple<>; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveModels.h new file mode 100644 index 00000000000..901b64a83e0 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveModels.h @@ -0,0 +1,34 @@ +/** + * \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 + * + */ + +#pragma once + +#include "MaterialLib/SolidModels/MechanicsBase.h" +#include "ProcessLib/Graph/ConstructModels.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SpecificBodyForceData.h" + +namespace ProcessLib::RichardsMechanics +{ +/// Constitutive models used for assembly. +template +using ConstitutiveModels = std::tuple<>; + +template +ConstitutiveModels createConstitutiveModels( + TRMProcessData const& process_data, + MaterialLib::Solids::MechanicsBase const& solid_material) +{ + return ProcessLib::Graph::constructModels< + ConstitutiveModels>( + ProcessLib::ThermoRichardsMechanics::SpecificBodyForceData< + DisplacementDim>{process_data.specific_body_force}, + solid_material); +} +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp new file mode 100644 index 00000000000..836ae0d3b7f --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp @@ -0,0 +1,35 @@ +/** + * \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 + * + */ + +#include "ConstitutiveSetting.h" + +namespace ProcessLib::RichardsMechanics +{ +template +void ConstitutiveSetting::eval( + ConstitutiveModels& models, double const t, + double const dt, ParameterLib::SpatialPosition const& x_position, + MaterialPropertyLib::Medium const& medium, + TemperatureData const& T_data, + CapillaryPressureData const& p_cap_data, + KelvinVector const& eps_arg, + StatefulData& state, + StatefulDataPrev const& prev_state, + ProcessLib::ThermoRichardsMechanics::MaterialStateData& + mat_state, + ConstitutiveTempData& tmp, + OutputData& out, + ConstitutiveData& cd) const +{ +} + +template struct ConstitutiveSetting<2>; +template struct ConstitutiveSetting<3>; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h new file mode 100644 index 00000000000..3dbae0f6195 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h @@ -0,0 +1,42 @@ +/** + * \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 + * + */ + +#pragma once + +#include "Base.h" +#include "ConstitutiveData.h" +#include "ConstitutiveModels.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/MaterialState.h" + +namespace ProcessLib::RichardsMechanics +{ +template +struct ConstitutiveSetting +{ + /// Evaluate the constitutive setting. + void eval( + ConstitutiveModels& models, double const t, + double const dt, ParameterLib::SpatialPosition const& x_position, + MaterialPropertyLib::Medium const& medium, + TemperatureData const& T_data, + CapillaryPressureData const& p_cap_data, + KelvinVector const& eps_arg, + StatefulData& state, + StatefulDataPrev const& prev_state, + ProcessLib::ThermoRichardsMechanics::MaterialStateData& + mat_state, + ConstitutiveTempData& tmp, + OutputData& out, + ConstitutiveData& cd) const; +}; + +extern template struct ConstitutiveSetting<2>; +extern template struct ConstitutiveSetting<3>; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.cpp b/ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.cpp new file mode 100644 index 00000000000..5ed1a9226ea --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.cpp @@ -0,0 +1,32 @@ +/** + * \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 + * + */ + +#include "CreateConstitutiveSetting.h" + +#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h" + +namespace ProcessLib::RichardsMechanics +{ +template +std::map>> +CreateConstitutiveSetting::createSolidConstitutiveRelations( + std::vector> const& parameters, + std::optional const& + local_coordinate_system, + BaseLib::ConfigTree const& config) +{ + return MaterialLib::Solids::createConstitutiveRelations( + parameters, local_coordinate_system, config); +} + +template struct CreateConstitutiveSetting<2>; +template struct CreateConstitutiveSetting<3>; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.h new file mode 100644 index 00000000000..3a1d56fb416 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/CreateConstitutiveSetting.h @@ -0,0 +1,33 @@ +/** + * \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 + * + */ + +#pragma once + +#include + +#include "MaterialLib/SolidModels/MechanicsBase.h" +#include "ParameterLib/Parameter.h" + +namespace ProcessLib::RichardsMechanics +{ +template +struct CreateConstitutiveSetting +{ + static std::map< + int, + std::unique_ptr>> + createSolidConstitutiveRelations( + std::vector> const& + parameters, + std::optional const& + local_coordinate_system, + BaseLib::ConfigTree const& config); +}; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index ca741dea3d3..87647190529 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -15,6 +15,7 @@ #include #include "ComputeMicroPorosity.h" +#include "ConstitutiveRelations/ConstitutiveModels.h" #include "IntegrationPointData.h" #include "MaterialLib/MPL/Medium.h" #include "MaterialLib/MPL/Utils/FormEigenTensor.h" @@ -774,6 +775,10 @@ void RichardsMechanicsLocalAssembler CD; + [[maybe_unused]] auto models = createConstitutiveModels( + _process_data, _ip_data[ip].solid_material); + x_position.setIntegrationPoint(ip); auto const& w = _ip_data[ip].integration_weight; From 2fe391a69a66b78306a0738af01e372276faa0d6 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Tue, 28 May 2024 12:11:56 +0200 Subject: [PATCH 2/2] [PL/RM] Extracted some constitutive setting evaluations ... from assembleWithJacobian() into a separate method. I added some constitutive data structs and reused some from TRM. --- .../ConstitutiveRelations/Base.h | 9 +- .../ConstitutiveRelations/ConstitutiveData.h | 23 +- .../ConstitutiveSetting.cpp | 24 +- .../ConstitutiveSetting.h | 3 +- .../ConstitutiveRelations/Density.h | 19 + .../EffectivePorePressure.h | 19 + .../ConstitutiveRelations/LiquidDensity.h | 18 + .../SaturationSecantDerivative.h | 21 + .../ConstitutiveRelations/StiffnessTensor.h | 20 + .../RichardsMechanicsFEM-impl.h | 424 +++++++++++------- .../RichardsMechanics/RichardsMechanicsFEM.h | 10 + 11 files changed, 410 insertions(+), 180 deletions(-) create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h create mode 100644 ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h index be08c171856..a70b6227fc9 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Base.h @@ -10,6 +10,7 @@ #pragma once +#include "BaseLib/StrongType.h" #include "MaterialLib/MPL/Medium.h" #include "MathLib/KelvinVector.h" #include "ParameterLib/SpatialPosition.h" @@ -63,13 +64,7 @@ struct MediaData MaterialPropertyLib::Phase const& solid; }; -template -struct TemperatureData -{ - double T; - double T_prev; - Eigen::Vector grad_T; -}; +using TemperatureData = BaseLib::StrongType; template struct CapillaryPressureData diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h index 8b750271bdb..e909848e70b 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h @@ -10,7 +10,17 @@ #pragma once +#include "Density.h" +#include "LiquidDensity.h" #include "ProcessLib/ConstitutiveRelations/Base.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Biot.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Saturation.h" +#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidCompressibilityData.h" +#include "SaturationSecantDerivative.h" +#include "StiffnessTensor.h" namespace ProcessLib::RichardsMechanics { @@ -29,7 +39,18 @@ using OutputData = std::tuple<>; /// Data that is needed for the equation system assembly. template -using ConstitutiveData = std::tuple<>; +using ConstitutiveData = std::tuple< + // TODO (CL) check if all that data should stay here + StiffnessTensor, + ProcessLib::ThermoRichardsMechanics::PorosityData, Density, LiquidDensity, + ProcessLib::ThermoRichardsMechanics::BiotData, + ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv, + ProcessLib::ThermoRichardsMechanics::LiquidViscosityData, + ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData, + ProcessLib::ThermoRichardsMechanics::BishopsData, + PrevState, + ProcessLib::ThermoRichardsMechanics::PermeabilityData, + SaturationSecantDerivative>; /// Data that stores intermediate values, which are not needed outside the /// constitutive setting. diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp index 836ae0d3b7f..63dd576bd85 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.cpp @@ -14,19 +14,19 @@ namespace ProcessLib::RichardsMechanics { template void ConstitutiveSetting::eval( - ConstitutiveModels& models, double const t, - double const dt, ParameterLib::SpatialPosition const& x_position, - MaterialPropertyLib::Medium const& medium, - TemperatureData const& T_data, - CapillaryPressureData const& p_cap_data, - KelvinVector const& eps_arg, - StatefulData& state, - StatefulDataPrev const& prev_state, + ConstitutiveModels& /*models*/, double const /*t*/, + double const /*dt*/, ParameterLib::SpatialPosition const& /*x_position*/, + MaterialPropertyLib::Medium const& /*medium*/, + TemperatureData const /*T_data*/, + CapillaryPressureData const& /*p_cap_data*/, + KelvinVector const& /*eps_arg*/, + StatefulData& /*state*/, + StatefulDataPrev const& /*prev_state*/, ProcessLib::ThermoRichardsMechanics::MaterialStateData& - mat_state, - ConstitutiveTempData& tmp, - OutputData& out, - ConstitutiveData& cd) const + /*mat_state*/, + ConstitutiveTempData& /*tmp*/, + OutputData& /*out*/, + ConstitutiveData& /*cd*/) const { } diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h index 3dbae0f6195..8e10f44ec91 100644 --- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveSetting.h @@ -24,8 +24,7 @@ struct ConstitutiveSetting void eval( ConstitutiveModels& models, double const t, double const dt, ParameterLib::SpatialPosition const& x_position, - MaterialPropertyLib::Medium const& medium, - TemperatureData const& T_data, + MaterialPropertyLib::Medium const& medium, TemperatureData const T_data, CapillaryPressureData const& p_cap_data, KelvinVector const& eps_arg, StatefulData& state, diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h new file mode 100644 index 00000000000..f53cd357b24 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/Density.h @@ -0,0 +1,19 @@ +/** + * \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 + * + */ + +#pragma once + +#include "BaseLib/StrongType.h" + +namespace ProcessLib::RichardsMechanics +{ +// effective density of both the solid and fluid phases +using Density = BaseLib::StrongType; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h new file mode 100644 index 00000000000..5df96875c01 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/EffectivePorePressure.h @@ -0,0 +1,19 @@ +/** + * \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 + * + */ + +#pragma once + +#include "BaseLib/StrongType.h" + +namespace ProcessLib::RichardsMechanics +{ +using EffectivePorePressure = + BaseLib::StrongType; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h new file mode 100644 index 00000000000..05b11b32c31 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/LiquidDensity.h @@ -0,0 +1,18 @@ +/** + * \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 + * + */ + +#pragma once + +#include "BaseLib/StrongType.h" + +namespace ProcessLib::RichardsMechanics +{ +using LiquidDensity = BaseLib::StrongType; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h new file mode 100644 index 00000000000..0e70b4cbc0f --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/SaturationSecantDerivative.h @@ -0,0 +1,21 @@ +/** + * \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 + * + */ + +#pragma once + +#include "Base.h" + +namespace ProcessLib::RichardsMechanics +{ +struct SaturationSecantDerivative +{ + double DeltaS_L_Deltap_cap = nan; +}; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h new file mode 100644 index 00000000000..46b272c2d95 --- /dev/null +++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/StiffnessTensor.h @@ -0,0 +1,20 @@ +/** + * \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 + * + */ + +#pragma once + +#include "Base.h" + +namespace ProcessLib::RichardsMechanics +{ +template +using StiffnessTensor = BaseLib::StrongType, + struct StiffnessTensorTag>; +} // namespace ProcessLib::RichardsMechanics diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h index 87647190529..54a0d49d9be 100644 --- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h +++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h @@ -685,6 +685,216 @@ void RichardsMechanicsLocalAssembler< } } +template +void RichardsMechanicsLocalAssembler:: + assembleWithJacobianEvalConstitutiveSetting( + double const t, double const dt, + ParameterLib::SpatialPosition const& x_position, + RichardsMechanicsLocalAssembler::IpData& ip_data, + MPL::VariableArray& variables, MPL::VariableArray& variables_prev, + MPL::Medium const* const medium, TemperatureData const T_data, + CapillaryPressureData const& p_cap_data, + ConstitutiveData& CD) +{ + auto const& liquid_phase = medium->phase("AqueousLiquid"); + auto const& solid_phase = medium->phase("Solid"); + + auto const& identity2 = MathLib::KelvinVector::Invariants< + MathLib::KelvinVector::kelvin_vector_dimensions( + DisplacementDim)>::identity2; + + double const temperature = T_data(); + double const p_cap_ip = p_cap_data.p_cap; + double const p_cap_prev_ip = p_cap_data.p_cap_prev; + + auto& eps = ip_data.eps; + auto& eps_m = ip_data.eps_m; + auto& S_L = ip_data.saturation; + auto const S_L_prev = ip_data.saturation_prev; + auto const alpha = + medium->property(MPL::PropertyType::biot_coefficient) + .template value(variables, x_position, t, dt); + *std::get(CD) = alpha; + + auto const C_el = + ip_data.computeElasticTangentStiffness(t, x_position, dt, temperature); + + auto const beta_SR = (1 - alpha) / ip_data.solid_material.getBulkModulus( + t, x_position, &C_el); + variables.grain_compressibility = beta_SR; + std::get(CD) + .beta_SR = beta_SR; + + auto const rho_LR = + liquid_phase.property(MPL::PropertyType::density) + .template value(variables, x_position, t, dt); + variables.density = rho_LR; + *std::get(CD) = rho_LR; + + S_L = medium->property(MPL::PropertyType::saturation) + .template value(variables, x_position, t, dt); + variables.liquid_saturation = S_L; + variables_prev.liquid_saturation = S_L_prev; + + // tangent derivative for Jacobian + double const dS_L_dp_cap = + medium->property(MPL::PropertyType::saturation) + .template dValue(variables, + MPL::Variable::capillary_pressure, + x_position, t, dt); + std::get(CD) + .dS_L_dp_cap = dS_L_dp_cap; + // secant derivative from time discretization for storage + // use tangent, if secant is not available + double const DeltaS_L_Deltap_cap = + (p_cap_ip == p_cap_prev_ip) + ? dS_L_dp_cap + : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip); + std::get(CD).DeltaS_L_Deltap_cap = + DeltaS_L_Deltap_cap; + + auto const chi = [medium, x_position, t, dt](double const S_L) + { + MPL::VariableArray vs; + vs.liquid_saturation = S_L; + return medium->property(MPL::PropertyType::bishops_effective_stress) + .template value(vs, x_position, t, dt); + }; + double const chi_S_L = chi(S_L); + std::get(CD).chi_S_L = + chi_S_L; + double const chi_S_L_prev = chi(S_L_prev); + std::get>(CD) + ->chi_S_L = chi_S_L_prev; + + auto const dchi_dS_L = + medium->property(MPL::PropertyType::bishops_effective_stress) + .template dValue( + variables, MPL::Variable::liquid_saturation, x_position, t, dt); + std::get(CD).dchi_dS_L = + dchi_dS_L; + + double const p_FR = -chi_S_L * p_cap_ip; + variables.effective_pore_pressure = p_FR; + variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip; + + // Set volumetric strain rate for the general case without swelling. + variables.volumetric_strain = Invariants::trace(ip_data.eps); + // TODO (CL) changed that, using eps_prev for the moment, not B * u_prev + // variables_prev.volumetric_strain = Invariants::trace(B * u_prev); + variables_prev.volumetric_strain = Invariants::trace(ip_data.eps_prev); + + auto& phi = ip_data.porosity; + { // Porosity update + + variables_prev.porosity = ip_data.porosity_prev; + phi = medium->property(MPL::PropertyType::porosity) + .template value(variables, variables_prev, x_position, + t, dt); + variables.porosity = phi; + } + std::get(CD).phi = phi; + + if (alpha < phi) + { + OGS_FATAL( + "RichardsMechanics: Biot-coefficient {} is smaller than " + "porosity {} in element/integration point {}/{}.", + alpha, phi, _element.getID(), -1 /*ip*/ /* TODO (CL) re-enable */); + } + + auto const mu = liquid_phase.property(MPL::PropertyType::viscosity) + .template value(variables, x_position, t, dt); + std::get(CD) + .viscosity = mu; + + // Swelling and possibly volumetric strain rate update. + updateSwellingStressAndVolumetricStrain( + ip_data, *medium, solid_phase, C_el, rho_LR, mu, + _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip, + variables, variables_prev, x_position, t, dt); + + if (medium->hasProperty(MPL::PropertyType::transport_porosity)) + { + if (!medium->hasProperty(MPL::PropertyType::saturation_micro)) + { + variables_prev.transport_porosity = ip_data.transport_porosity_prev; + + ip_data.transport_porosity = + medium->property(MPL::PropertyType::transport_porosity) + .template value(variables, variables_prev, + x_position, t, dt); + variables.transport_porosity = ip_data.transport_porosity; + } + } + else + { + variables.transport_porosity = phi; + } + + // Set mechanical variables for the intrinsic permeability model + // For stress dependent permeability. + { + auto const sigma_total = + (ip_data.sigma_eff + alpha * p_FR * identity2).eval(); + // For stress dependent permeability. + variables.total_stress.emplace( + MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total)); + } + + variables.equivalent_plastic_strain = + ip_data.material_state_variables->getEquivalentPlasticStrain(); + + double const k_rel = + medium->property(MPL::PropertyType::relative_permeability) + .template value(variables, x_position, t, dt); + + auto const K_intrinsic = MPL::formEigenTensor( + medium->property(MPL::PropertyType::permeability) + .value(variables, x_position, t, dt)); + + std::get< + ProcessLib::ThermoRichardsMechanics::PermeabilityData>( + CD) + .k_rel = k_rel; + std::get< + ProcessLib::ThermoRichardsMechanics::PermeabilityData>( + CD) + .Ki = K_intrinsic; + + // + // displacement equation, displacement part + // + auto& sigma_sw = ip_data.sigma_sw; + + eps_m.noalias() = + solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) + ? eps + C_el.inverse() * sigma_sw + : eps; + variables.mechanical_strain + .emplace>( + eps_m); + + auto C = ip_data.updateConstitutiveRelation(variables, t, x_position, dt, + temperature); + + *std::get>(CD) = std::move(C); + + // p_SR + variables.solid_grain_pressure = + p_FR - ip_data.sigma_eff.dot(identity2) / (3 * (1 - phi)); + auto const rho_SR = + solid_phase.property(MPL::PropertyType::density) + .template value(variables, x_position, t, dt); + + double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR; + *std::get(CD) = rho; +} + template void RichardsMechanicsLocalAssembler CD; + ConstitutiveData CD; [[maybe_unused]] auto models = createConstitutiveModels( _process_data, _ip_data[ip].solid_material); @@ -816,166 +1026,25 @@ void RichardsMechanicsLocalAssemblerproperty(MPL::PropertyType::biot_coefficient) - .template value(variables, x_position, t, dt); - auto const C_el = _ip_data[ip].computeElasticTangentStiffness( - t, x_position, dt, temperature); - - auto const beta_SR = - (1 - alpha) / - _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el); - variables.grain_compressibility = beta_SR; - - auto const rho_LR = - liquid_phase.property(MPL::PropertyType::density) - .template value(variables, x_position, t, dt); - variables.density = rho_LR; - - auto const& b = _process_data.specific_body_force; - - S_L = medium->property(MPL::PropertyType::saturation) - .template value(variables, x_position, t, dt); - variables.liquid_saturation = S_L; - variables_prev.liquid_saturation = S_L_prev; - - // tangent derivative for Jacobian - double const dS_L_dp_cap = - medium->property(MPL::PropertyType::saturation) - .template dValue(variables, - MPL::Variable::capillary_pressure, - x_position, t, dt); - // secant derivative from time discretization for storage - // use tangent, if secant is not available - double const DeltaS_L_Deltap_cap = - (p_cap_ip == p_cap_prev_ip) - ? dS_L_dp_cap - : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip); - - auto const chi = [medium, x_position, t, dt](double const S_L) - { - MPL::VariableArray vs; - vs.liquid_saturation = S_L; - return medium->property(MPL::PropertyType::bishops_effective_stress) - .template value(vs, x_position, t, dt); - }; - double const chi_S_L = chi(S_L); - double const chi_S_L_prev = chi(S_L_prev); - - double const p_FR = -chi_S_L * p_cap_ip; - variables.effective_pore_pressure = p_FR; - variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip; - - // Set volumetric strain rate for the general case without swelling. - variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps); - variables_prev.volumetric_strain = Invariants::trace(B * u_prev); - - auto& phi = _ip_data[ip].porosity; - { // Porosity update - - variables_prev.porosity = _ip_data[ip].porosity_prev; - phi = medium->property(MPL::PropertyType::porosity) - .template value(variables, variables_prev, - x_position, t, dt); - variables.porosity = phi; - } - - if (alpha < phi) - { - OGS_FATAL( - "RichardsMechanics: Biot-coefficient {} is smaller than " - "porosity {} in element/integration point {}/{}.", - alpha, phi, _element.getID(), ip); - } - - auto const mu = - liquid_phase.property(MPL::PropertyType::viscosity) - .template value(variables, x_position, t, dt); - - // Swelling and possibly volumetric strain rate update. - updateSwellingStressAndVolumetricStrain( - _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu, - _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip, - variables, variables_prev, x_position, t, dt); - - if (medium->hasProperty(MPL::PropertyType::transport_porosity)) - { - if (!medium->hasProperty(MPL::PropertyType::saturation_micro)) - { - variables_prev.transport_porosity = - _ip_data[ip].transport_porosity_prev; - - _ip_data[ip].transport_porosity = - medium->property(MPL::PropertyType::transport_porosity) - .template value(variables, variables_prev, - x_position, t, dt); - variables.transport_porosity = _ip_data[ip].transport_porosity; - } - } - else - { - variables.transport_porosity = phi; - } - - double const k_rel = - medium->property(MPL::PropertyType::relative_permeability) - .template value(variables, x_position, t, dt); - - // Set mechanical variables for the intrinsic permeability model - // For stress dependent permeability. - { - auto const sigma_total = - (_ip_data[ip].sigma_eff + alpha * p_FR * identity2).eval(); - // For stress dependent permeability. - variables.total_stress.emplace( - MathLib::KelvinVector::kelvinVectorToSymmetricTensor( - sigma_total)); - } - - variables.equivalent_plastic_strain = - _ip_data[ip].material_state_variables->getEquivalentPlasticStrain(); - - auto const K_intrinsic = MPL::formEigenTensor( - medium->property(MPL::PropertyType::permeability) - .value(variables, x_position, t, dt)); - - GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu; - - // - // displacement equation, displacement part - // - auto& sigma_sw = _ip_data[ip].sigma_sw; - - eps_m.noalias() = - solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate) - ? eps + C_el.inverse() * sigma_sw - : eps; - variables.mechanical_strain - .emplace>( - eps_m); - - auto C = _ip_data[ip].updateConstitutiveRelation( - variables, t, x_position, dt, temperature); + assembleWithJacobianEvalConstitutiveSetting( + t, dt, x_position, _ip_data[ip], variables, variables_prev, medium, + TemperatureData{temperature}, + CapillaryPressureData{ + p_cap_ip, p_cap_prev_ip, + Eigen::Vector::Zero()}, + CD); + auto const& C = *std::get>(CD); local_Jac .template block( displacement_index, displacement_index) .noalias() += B.transpose() * C * B * w; - // p_SR - variables.solid_grain_pressure = - p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi)); - auto const rho_SR = - solid_phase.property(MPL::PropertyType::density) - .template value(variables, x_position, t, dt); - - double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR; + auto const& sigma_eff = _ip_data[ip].sigma_eff; + double const rho = *std::get(CD); + auto const& b = _process_data.specific_body_force; local_rhs.template segment(displacement_index) .noalias() -= (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w; @@ -983,13 +1052,21 @@ void RichardsMechanicsLocalAssembler(CD); + double const chi_S_L = + std::get(CD) + .chi_S_L; Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w; - auto const dchi_dS_L = - medium->property(MPL::PropertyType::bishops_effective_stress) - .template dValue(variables, - MPL::Variable::liquid_saturation, - x_position, t, dt); + double const dchi_dS_L = + std::get(CD) + .dchi_dS_L; + double const dS_L_dp_cap = + std::get( + CD) + .dS_L_dp_cap; local_Jac .template block( displacement_index, pressure_index) @@ -997,6 +1074,9 @@ void RichardsMechanicsLocalAssembler(CD).phi; + double const rho_LR = *std::get(CD); local_Jac .template block( displacement_index, pressure_index) @@ -1033,8 +1113,12 @@ void RichardsMechanicsLocalAssembler>(CD) + ->chi_S_L; Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha * identity2.transpose() * B * w; } @@ -1047,6 +1131,22 @@ void RichardsMechanicsLocalAssembler>(CD) + .k_rel; + auto const& K_intrinsic = + std::get>(CD) + .Ki; + double const mu = + std::get( + CD) + .viscosity; + + GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu; + laplace_p.noalias() += dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w; @@ -1057,6 +1157,11 @@ void RichardsMechanicsLocalAssembler( + CD) + .beta_SR; double const a0 = (alpha - phi) * beta_SR; double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0); double const specific_storage_a_S = phi - p_cap_ip * S_L * a0; @@ -1069,6 +1174,8 @@ void RichardsMechanicsLocalAssembler(CD).DeltaS_L_Deltap_cap; storage_p_a_S.noalias() -= N_p.transpose() * rho_LR * specific_storage_a_S * DeltaS_L_Deltap_cap * N_p * w; @@ -1079,6 +1186,7 @@ void RichardsMechanicsLocalAssembler #include +#include "ConstitutiveRelations/Base.h" +#include "ConstitutiveRelations/ConstitutiveData.h" #include "IntegrationPointData.h" #include "LocalAssemblerInterface.h" #include "MaterialLib/MPL/VariableType.h" @@ -328,6 +330,14 @@ class RichardsMechanicsLocalAssembler getMaterialStateVariablesAt(unsigned integration_point) const override; private: + void assembleWithJacobianEvalConstitutiveSetting( + double const t, double const dt, + ParameterLib::SpatialPosition const& x_position, IpData& ip_data, + MPL::VariableArray& variables, MPL::VariableArray& variables_prev, + MPL::Medium const* const medium, TemperatureData const T_data, + CapillaryPressureData const& p_cap_data, + ConstitutiveData& CD); + RichardsMechanicsProcessData& _process_data; std::vector> _ip_data;