Loading [MathJax]/extensions/tex2jax.js
OGS
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim > Class Template Reference

Detailed Description

template<typename ShapeFunction, int GlobalDim>
class ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >

Definition at line 36 of file ThermoRichardsFlowFEM.h.

#include <ThermoRichardsFlowFEM.h>

Inheritance diagram for ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >:
[legend]
Collaboration diagram for ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >:
[legend]

Public Types

using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
using IpData = IntegrationPointData<ShapeMatricesType>

Public Member Functions

 ThermoRichardsFlowLocalAssembler (ThermoRichardsFlowLocalAssembler const &)=delete
 ThermoRichardsFlowLocalAssembler (ThermoRichardsFlowLocalAssembler &&)=delete
 ThermoRichardsFlowLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, ThermoRichardsFlowProcessData &process_data)
std::size_t setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
void setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override
void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_rhs_data) override
void initializeConcrete () override
void postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
void computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
std::vector< double > const & getIntPtDarcyVelocity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > getSaturation () const override
std::vector< double > const & getIntPtSaturation (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > getPorosity () const override
std::vector< double > const & getIntPtPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const & getIntPtDryDensitySolid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
virtual void setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void preAssemble (double const, double const, std::vector< double > const &)
virtual void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
virtual std::optional< VectorSegment > getVectorDeformationSegment () const
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default

Private Member Functions

unsigned getNumberOfIntegrationPoints () const override

Private Attributes

ThermoRichardsFlowProcessData & _process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
NumLib::GenericIntegrationMethod const & _integration_method
MeshLib::Element const & _element
bool const _is_axially_symmetric

Static Private Attributes

static const int temperature_index = 0
static const int temperature_size = ShapeFunction::NPOINTS
static const int pressure_index = temperature_size
static const int pressure_size = ShapeFunction::NPOINTS

Member Typedef Documentation

◆ GlobalDimMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

Definition at line 43 of file ThermoRichardsFlowFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType

Definition at line 44 of file ThermoRichardsFlowFEM.h.

◆ IpData

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::IpData = IntegrationPointData<ShapeMatricesType>

Definition at line 46 of file ThermoRichardsFlowFEM.h.

◆ ShapeMatricesType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>

Definition at line 41 of file ThermoRichardsFlowFEM.h.

Constructor & Destructor Documentation

◆ ThermoRichardsFlowLocalAssembler() [1/3]

template<typename ShapeFunction, int GlobalDim>
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler ( ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim > const & )
delete

References ThermoRichardsFlowLocalAssembler().

Referenced by ThermoRichardsFlowLocalAssembler(), and ThermoRichardsFlowLocalAssembler().

◆ ThermoRichardsFlowLocalAssembler() [2/3]

template<typename ShapeFunction, int GlobalDim>
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler ( ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim > && )
delete

References ThermoRichardsFlowLocalAssembler().

◆ ThermoRichardsFlowLocalAssembler() [3/3]

template<typename ShapeFunction, int GlobalDim>
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler ( MeshLib::Element const & e,
std::size_t const ,
NumLib::GenericIntegrationMethod const & integration_method,
bool const is_axially_symmetric,
ThermoRichardsFlowProcessData & process_data )

Definition at line 34 of file ThermoRichardsFlowFEM-impl.h.

41 : _process_data(process_data),
42 _integration_method(integration_method),
43 _element(e),
44 _is_axially_symmetric(is_axially_symmetric)
45{
46 unsigned const n_integration_points =
47 _integration_method.getNumberOfPoints();
48
49 _ip_data.reserve(n_integration_points);
50
51 auto const shape_matrices =
52 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim>(
53 e, is_axially_symmetric, _integration_method);
54
55 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
56
57 for (unsigned ip = 0; ip < n_integration_points; ip++)
58 {
59 auto const& sm = shape_matrices[ip];
60 _ip_data.emplace_back();
61 auto& ip_data = _ip_data[ip];
62 _ip_data[ip].integration_weight =
63 _integration_method.getWeightedPoint(ip).getWeight() *
64 sm.integralMeasure * sm.detJ;
65
66 ip_data.N = sm.N;
67 ip_data.dNdx = sm.dNdx;
68
69 ParameterLib::SpatialPosition const x_position{
70 std::nullopt, _element.getID(),
71 MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
72 ShapeMatricesType>(
73 _element, sm.N))};
74 // Initial porosity. Could be read from integration point data or mesh.
75 ip_data.porosity = medium[MPL::porosity].template initialValue<double>(
76 x_position,
77 std::numeric_limits<double>::quiet_NaN() /* t independent */);
78 }
79}
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler::ShapeMatricesType
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
Definition ThermoRichardsFlowFEM.h:41
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition ThermoRichardsFlowFEM.h:154
NumLib::initShapeMatrices
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Definition InitShapeMatrices.h:54

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, NumLib::initShapeMatrices(), NumLib::interpolateCoordinates(), and MaterialPropertyLib::porosity.

Member Function Documentation

◆ assemble()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::assemble ( double const t,
double const dt,
std::vector< double > const & local_x,
std::vector< double > const & local_x_prev,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > & local_rhs_data )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 687 of file ThermoRichardsFlowFEM-impl.h.

691{
692 auto const local_matrix_dim = pressure_size + temperature_size;
693 assert(local_x.size() == local_matrix_dim);
694
695 auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
696 temperature_size> const>(local_x.data() + temperature_index,
697 temperature_size);
698 auto const p_L = Eigen::Map<
699 typename ShapeMatricesType::template VectorType<pressure_size> const>(
700 local_x.data() + pressure_index, pressure_size);
701
702 auto const p_L_prev = Eigen::Map<
703 typename ShapeMatricesType::template VectorType<pressure_size> const>(
704 local_x_prev.data() + pressure_index, pressure_size);
705
706 auto local_K = MathLib::createZeroedMatrix<
707 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
708 local_matrix_dim>>(
709 local_K_data, local_matrix_dim, local_matrix_dim);
710
711 auto local_M = MathLib::createZeroedMatrix<
712 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
713 local_matrix_dim>>(
714 local_M_data, local_matrix_dim, local_matrix_dim);
715
716 auto local_rhs = MathLib::createZeroedVector<
717 typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
718 local_rhs_data, local_matrix_dim);
719
720 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
721 auto const& liquid_phase = medium.phase("AqueousLiquid");
722 auto const& solid_phase = medium.phase("Solid");
723 MPL::Phase const* gas_phase =
724 medium.hasPhase("Gas") ? &medium.phase("Gas") : nullptr;
725 MPL::VariableArray variables;
726 MPL::VariableArray variables_prev;
727
728 unsigned const n_integration_points =
729 _integration_method.getNumberOfPoints();
730 for (unsigned ip = 0; ip < n_integration_points; ip++)
731 {
732 auto const& w = _ip_data[ip].integration_weight;
733
734 auto const& N = _ip_data[ip].N;
735 auto const& dNdx = _ip_data[ip].dNdx;
736
737 ParameterLib::SpatialPosition const x_position{
738 std::nullopt, _element.getID(),
739 MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
740 ShapeMatricesType>(
741 _element, N))};
742
743 double T_ip;
744 NumLib::shapeFunctionInterpolate(T, N, T_ip);
745
746 double p_cap_ip;
747 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
748
749 double p_cap_prev_ip;
750 NumLib::shapeFunctionInterpolate(-p_L_prev, N, p_cap_prev_ip);
751
752 variables.capillary_pressure = p_cap_ip;
753 variables.liquid_phase_pressure = -p_cap_ip;
754 // setting pG to 1 atm
755 // TODO : rewrite equations s.t. p_L = pG-p_cap
756 variables.gas_phase_pressure = 1.0e5;
757 variables.temperature = T_ip;
758
759 auto& S_L = _ip_data[ip].saturation;
760 auto const S_L_prev = _ip_data[ip].saturation_prev;
761 auto const alpha =
762 medium[MPL::PropertyType::biot_coefficient].template value<double>(
763 variables, x_position, t, dt);
764
765 auto& solid_elasticity = *_process_data.simplified_elasticity;
766 // TODO (buchwaldj)
767 // is bulk_modulus good name for bulk modulus of solid skeleton?
768 auto const beta_S =
769 solid_elasticity.bulkCompressibilityFromYoungsModulus(
770 solid_phase, variables, x_position, t, dt);
771 auto const beta_SR = (1 - alpha) * beta_S;
772 variables.grain_compressibility = beta_SR;
773
774 auto const rho_LR =
775 liquid_phase[MPL::PropertyType::density].template value<double>(
776 variables, x_position, t, dt);
777 auto const& b = _process_data.specific_body_force;
778
779 double const drho_LR_dp =
780 liquid_phase[MPL::PropertyType::density].template dValue<double>(
781 variables, MPL::Variable::liquid_phase_pressure, x_position, t,
782 dt);
783 auto const beta_LR = drho_LR_dp / rho_LR;
784
785 S_L = medium[MPL::PropertyType::saturation].template value<double>(
786 variables, x_position, t, dt);
787 variables.liquid_saturation = S_L;
788 variables_prev.liquid_saturation = S_L_prev;
789
790 // tangent derivative for Jacobian
791 double const dS_L_dp_cap =
792 medium[MPL::PropertyType::saturation].template dValue<double>(
793 variables, MPL::Variable::capillary_pressure, x_position, t,
794 dt);
795 // secant derivative from time discretization for storage
796 // use tangent, if secant is not available
797 double const DeltaS_L_Deltap_cap =
798 (p_cap_ip == p_cap_prev_ip)
799 ? dS_L_dp_cap
800 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
801
802 auto chi_S_L = S_L;
803 auto chi_S_L_prev = S_L_prev;
804 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
805 {
806 auto const chi = [&medium, x_position, t, dt](double const S_L)
807 {
808 MPL::VariableArray variables;
809 variables.liquid_saturation = S_L;
810 return medium[MPL::PropertyType::bishops_effective_stress]
811 .template value<double>(variables, x_position, t, dt);
812 };
813 chi_S_L = chi(S_L);
814 chi_S_L_prev = chi(S_L_prev);
815 }
816 // TODO (buchwaldj)
817 // should solid_grain_pressure or effective_pore_pressure remain?
818 // double const p_FR = -chi_S_L * p_cap_ip;
819 // variables.solid_grain_pressure = p_FR;
820
821 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
822 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
823
824 auto& phi = _ip_data[ip].porosity;
825 { // Porosity update
826
827 variables_prev.porosity = _ip_data[ip].porosity_prev;
828 phi = medium[MPL::PropertyType::porosity].template value<double>(
829 variables, variables_prev, x_position, t, dt);
830 variables.porosity = phi;
831 }
832
833 if (alpha < phi)
834 {
835 OGS_FATAL(
836 "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
837 "porosity {} in element/integration point {}/{}.",
838 alpha, phi, _element.getID(), ip);
839 }
840
841 auto const K_pT_thermal_osmosis =
842 (solid_phase.hasProperty(
843 MaterialPropertyLib::PropertyType::thermal_osmosis_coefficient)
844 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
845 solid_phase
846 [MPL::PropertyType::thermal_osmosis_coefficient]
847 .value(variables, x_position, t, dt))
848 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
849
850 double const k_rel =
851 medium[MPL::PropertyType::relative_permeability]
852 .template value<double>(variables, x_position, t, dt);
853 auto const mu =
854 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
855 variables, x_position, t, dt);
856
857 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
858 medium[MPL::PropertyType::permeability].value(variables, x_position,
859 t, dt));
860
861 GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
862 GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
863
864 // Consider anisotropic thermal expansion.
865 // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
866 // component.
867 Eigen::Matrix<double, 3, 3> const
868 solid_linear_thermal_expansion_coefficient =
869 MaterialPropertyLib::formEigenTensor<3>(
870 solid_phase
871 [MaterialPropertyLib::PropertyType::thermal_expansivity]
872 .value(variables, x_position, t, dt));
873
874 auto const rho_SR =
875 solid_phase[MPL::PropertyType::density].template value<double>(
876 variables, x_position, t, dt);
877
878 //
879 // pressure equation, pressure part.
880 //
881 local_K
882 .template block<pressure_size, pressure_size>(pressure_index,
883 pressure_index)
884 .noalias() += dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
885
886 const double alphaB_minus_phi = alpha - phi;
887 double const a0 = alphaB_minus_phi * beta_SR;
888 double const specific_storage_a_p =
889 S_L * (phi * beta_LR + S_L * a0 +
890 chi_S_L * alpha * alpha *
891 solid_elasticity.storageContribution(
892 solid_phase, variables, x_position, t, dt));
893 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
894
895 local_M
896 .template block<pressure_size, pressure_size>(pressure_index,
897 pressure_index)
898 .noalias() += N.transpose() * rho_LR *
899 (specific_storage_a_p -
900 specific_storage_a_S * DeltaS_L_Deltap_cap) *
901 N * w;
902
903 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
904 dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
905
906 //
907 // pressure equation, temperature part.
908 //
909 double const fluid_volumetric_thermal_expansion_coefficient =
910 MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
911 x_position, t, dt);
912 const double eff_thermal_expansion =
913 S_L * (alphaB_minus_phi *
914 solid_linear_thermal_expansion_coefficient.trace() +
915 phi * fluid_volumetric_thermal_expansion_coefficient +
916 alpha * solid_elasticity.thermalExpansivityContribution(
917 solid_linear_thermal_expansion_coefficient,
918 solid_phase, variables, x_position, t, dt));
919
920 local_K
921 .template block<pressure_size, temperature_size>(pressure_index,
922 temperature_index)
923 .noalias() +=
924 dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
925
926 local_M
927 .template block<pressure_size, temperature_size>(pressure_index,
928 temperature_index)
929 .noalias() -=
930 N.transpose() * rho_LR * eff_thermal_expansion * N * w;
931
932 //
933 // temperature equation.
934 //
935 {
936 auto const specific_heat_capacity_fluid =
937 liquid_phase[MaterialPropertyLib::specific_heat_capacity]
938 .template value<double>(variables, x_position, t, dt);
939
940 auto const specific_heat_capacity_solid =
941 solid_phase
942 [MaterialPropertyLib::PropertyType::specific_heat_capacity]
943 .template value<double>(variables, x_position, t, dt);
944
945 local_M
946 .template block<temperature_size, temperature_size>(
947 temperature_index, temperature_index)
948 .noalias() +=
949 w *
950 (rho_SR * specific_heat_capacity_solid * (1 - phi) +
951 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
952 N.transpose() * N;
953
954 auto const thermal_conductivity =
955 MaterialPropertyLib::formEigenTensor<GlobalDim>(
956 medium[MaterialPropertyLib::PropertyType::
957 thermal_conductivity]
958 .value(variables, x_position, t, dt));
959
960 GlobalDimVectorType const velocity_L = GlobalDimVectorType(
961 -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
962 K_pT_thermal_osmosis * dNdx * T);
963
964 local_K
965 .template block<temperature_size, temperature_size>(
966 temperature_index, temperature_index)
967 .noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
968 N.transpose() * velocity_L.transpose() * dNdx *
969 rho_LR * specific_heat_capacity_fluid) *
970 w;
971 local_K
972 .template block<temperature_size, pressure_size>(
973 temperature_index, pressure_index)
974 .noalias() +=
975 dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
976 }
977 if (gas_phase && S_L < 1.0)
978 {
979 variables.density = rho_LR;
980
981 double const rho_wv =
982 gas_phase->property(MaterialPropertyLib::density)
983 .template value<double>(variables, x_position, t, dt);
984
985 double const drho_wv_dT =
986 gas_phase->property(MaterialPropertyLib::density)
987 .template dValue<double>(variables,
988 MPL::Variable::temperature,
989 x_position, t, dt);
990 double const drho_wv_dp =
991 gas_phase->property(MaterialPropertyLib::density)
992 .template dValue<double>(
993 variables, MPL::Variable::liquid_phase_pressure,
994 x_position, t, dt);
995 auto const f_Tv =
996 gas_phase
997 ->property(
998 MPL::PropertyType::thermal_diffusion_enhancement_factor)
999 .template value<double>(variables, x_position, t, dt);
1000
1001 variables.porosity = phi;
1002 auto const tortuosity =
1003 medium.property(MPL::PropertyType::tortuosity)
1004 .template value<double>(variables, x_position, t, dt);
1005 double const D_v =
1006 phi * (1.0 - S_L) * tortuosity *
1007 gas_phase->property(MPL::PropertyType::diffusion)
1008 .template value<double>(variables, x_position, t, dt);
1009
1010 double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
1011 double const D_pv = D_v * drho_wv_dp;
1012
1013 GlobalDimVectorType const grad_T = dNdx * T;
1014 GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
1015 GlobalDimVectorType const vapour_flux =
1016 -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
1017 double const specific_heat_capacity_vapour =
1018 gas_phase->property(MaterialPropertyLib::specific_heat_capacity)
1019 .template value<double>(variables, x_position, t, dt);
1020
1021 local_M
1022 .template block<temperature_size, temperature_size>(
1023 temperature_index, temperature_index)
1024 .noalias() +=
1025 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
1026 N.transpose() * N;
1027
1028 local_K
1029 .template block<temperature_size, temperature_size>(
1030 temperature_index, temperature_index)
1031 .noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
1032 specific_heat_capacity_vapour * w;
1033
1034 double const storage_coefficient_by_water_vapor =
1035 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
1036 local_M
1037 .template block<pressure_size, pressure_size>(pressure_index,
1038 pressure_index)
1039 .noalias() +=
1040 N.transpose() * storage_coefficient_by_water_vapor * N * w;
1041
1042 double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
1043 local_M
1044 .template block<pressure_size, temperature_size>(
1045 pressure_index, temperature_index)
1046 .noalias() += N.transpose() * vapor_expansion_factor * N * w;
1047
1048 local_rhs.template segment<pressure_size>(pressure_index)
1049 .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
1050
1051 local_K
1052 .template block<pressure_size, pressure_size>(pressure_index,
1053 pressure_index)
1054 .noalias() += dNdx.transpose() * D_pv * dNdx * w;
1055
1056 //
1057 // Latent heat term
1058 //
1059 if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
1060 {
1061 double const factor = phi * (1 - S_L) / rho_LR;
1062 // The volumetric latent heat of vaporization of liquid water
1063 double const L0 =
1064 gas_phase->property(MPL::PropertyType::specific_latent_heat)
1065 .template value<double>(variables, x_position, t, dt) *
1066 rho_LR;
1067
1068 double const drho_LR_dT =
1069 liquid_phase.property(MPL::PropertyType::density)
1070 .template dValue<double>(variables,
1071 MPL::Variable::temperature,
1072 x_position, t, dt);
1073
1074 double const rho_wv_over_rho_L = rho_wv / rho_LR;
1075 local_M
1076 .template block<temperature_size, temperature_size>(
1077 temperature_index, temperature_index)
1078 .noalias() +=
1079 factor * L0 *
1080 (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
1081 N.transpose() * N * w;
1082
1083 local_M
1084 .template block<temperature_size, pressure_size>(
1085 temperature_index, pressure_index)
1086 .noalias() +=
1087 (factor * L0 *
1088 (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
1089 L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
1090 N.transpose() * N * w;
1091
1092 // temperature equation, temperature part
1093 local_K
1094 .template block<temperature_size, temperature_size>(
1095 temperature_index, temperature_index)
1096 .noalias() +=
1097 L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
1098 // temperature equation, pressure part
1099 local_K
1100 .template block<temperature_size, pressure_size>(
1101 temperature_index, pressure_index)
1102 .noalias() +=
1103 L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
1104 }
1105 }
1106 }
1107
1108 if (_process_data.apply_mass_lumping)
1109 {
1110 auto Mpp = local_M.template block<pressure_size, pressure_size>(
1111 pressure_index, pressure_index);
1112 Mpp = Mpp.colwise().sum().eval().asDiagonal();
1113 }
1114}
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler::GlobalDimVectorType
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
Definition ThermoRichardsFlowFEM.h:44
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler::GlobalDimMatrixType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
Definition ThermoRichardsFlowFEM.h:43
MaterialPropertyLib::formEigenTensor< 3 >
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
MaterialPropertyLib::formEigenTensor
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.cpp:126
MaterialPropertyLib::getLiquidThermalExpansivity
double getLiquidThermalExpansivity(Phase const &phase, VariableArray const &vars, const double density, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Definition GetLiquidThermalExpansivity.cpp:21
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:27

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::VariableArray::capillary_pressure, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::diffusion, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::formEigenTensor< 3 >(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getLiquidThermalExpansivity(), MaterialPropertyLib::VariableArray::grain_compressibility, MaterialPropertyLib::Phase::hasProperty(), NumLib::interpolateCoordinates(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, OGS_FATAL, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, MaterialPropertyLib::Phase::property(), MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::specific_heat_capacity, MaterialPropertyLib::specific_latent_heat, MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, temperature_index, temperature_size, MaterialPropertyLib::thermal_diffusion_enhancement_factor, MaterialPropertyLib::thermal_expansivity, MaterialPropertyLib::thermal_osmosis_coefficient, MaterialPropertyLib::tortuosity, and MaterialPropertyLib::viscosity.

◆ assembleWithJacobian()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::assembleWithJacobian ( double const t,
double const dt,
std::vector< double > const & local_x,
std::vector< double > const & local_x_prev,
std::vector< double > & local_rhs_data,
std::vector< double > & local_Jac_data )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 154 of file ThermoRichardsFlowFEM-impl.h.

160{
161 auto const local_matrix_dim = pressure_size + temperature_size;
162 assert(local_x.size() == local_matrix_dim);
163
164 auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
165 temperature_size> const>(local_x.data() + temperature_index,
166 temperature_size);
167 auto const p_L = Eigen::Map<
168 typename ShapeMatricesType::template VectorType<pressure_size> const>(
169 local_x.data() + pressure_index, pressure_size);
170
171 auto const T_prev =
172 Eigen::Map<typename ShapeMatricesType::template VectorType<
173 temperature_size> const>(local_x_prev.data() + temperature_index,
174 temperature_size);
175 auto const p_L_prev = Eigen::Map<
176 typename ShapeMatricesType::template VectorType<pressure_size> const>(
177 local_x_prev.data() + pressure_index, pressure_size);
178
179 auto local_Jac = MathLib::createZeroedMatrix<
180 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
181 local_matrix_dim>>(
182 local_Jac_data, local_matrix_dim, local_matrix_dim);
183
184 auto local_rhs = MathLib::createZeroedVector<
185 typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
186 local_rhs_data, local_matrix_dim);
187
188 typename ShapeMatricesType::NodalMatrixType M_TT =
189 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
190 temperature_size);
191 typename ShapeMatricesType::NodalMatrixType M_Tp =
192 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
193 pressure_size);
194 typename ShapeMatricesType::NodalMatrixType K_TT =
195 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
196 temperature_size);
197 typename ShapeMatricesType::NodalMatrixType K_Tp =
198 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
199 pressure_size);
200 typename ShapeMatricesType::NodalMatrixType dK_TT_dp =
201 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
202 pressure_size);
203 typename ShapeMatricesType::NodalMatrixType M_pT =
204 ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
205 temperature_size);
206 typename ShapeMatricesType::NodalMatrixType laplace_p =
207 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
208 typename ShapeMatricesType::NodalMatrixType laplace_T =
209 ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
210 temperature_size);
211 typename ShapeMatricesType::NodalMatrixType storage_p_a_p =
212 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
213
214 typename ShapeMatricesType::NodalMatrixType storage_p_a_S_Jpp =
215 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
216
217 typename ShapeMatricesType::NodalMatrixType storage_p_a_S =
218 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
219
220 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
221 auto const& liquid_phase = medium.phase("AqueousLiquid");
222 auto const& solid_phase = medium.phase("Solid");
223 MPL::Phase const* gas_phase =
224 medium.hasPhase("Gas") ? &medium.phase("Gas") : nullptr;
225 MPL::VariableArray variables;
226 MPL::VariableArray variables_prev;
227
228 unsigned const n_integration_points =
229 _integration_method.getNumberOfPoints();
230 for (unsigned ip = 0; ip < n_integration_points; ip++)
231 {
232 auto const& w = _ip_data[ip].integration_weight;
233
234 auto const& N = _ip_data[ip].N;
235 auto const& dNdx = _ip_data[ip].dNdx;
236
237 ParameterLib::SpatialPosition const x_position{
238 std::nullopt, _element.getID(),
239 MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
240 ShapeMatricesType>(
241 _element, N))};
242
243 double T_ip;
244 NumLib::shapeFunctionInterpolate(T, N, T_ip);
245
246 double p_cap_ip;
247 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
248
249 double p_cap_prev_ip;
250 NumLib::shapeFunctionInterpolate(-p_L_prev, N, p_cap_prev_ip);
251
252 variables.capillary_pressure = p_cap_ip;
253 variables.liquid_phase_pressure = -p_cap_ip;
254 // setting pG to 1 atm
255 // TODO : rewrite equations s.t. p_L = pG-p_cap
256 variables.gas_phase_pressure = 1.0e5;
257 variables.temperature = T_ip;
258
259 auto& S_L = _ip_data[ip].saturation;
260 auto const S_L_prev = _ip_data[ip].saturation_prev;
261 auto const alpha =
262 medium[MPL::PropertyType::biot_coefficient].template value<double>(
263 variables, x_position, t, dt);
264
265 auto& solid_elasticity = *_process_data.simplified_elasticity;
266 // TODO (buchwaldj)
267 // is bulk_modulus good name for bulk modulus of solid skeleton?
268 auto const beta_S =
269 solid_elasticity.bulkCompressibilityFromYoungsModulus(
270 solid_phase, variables, x_position, t, dt);
271 auto const beta_SR = (1 - alpha) * beta_S;
272 variables.grain_compressibility = beta_SR;
273
274 auto const rho_LR =
275 liquid_phase[MPL::PropertyType::density].template value<double>(
276 variables, x_position, t, dt);
277 variables.density = rho_LR;
278 auto const& b = _process_data.specific_body_force;
279
280 double const drho_LR_dp =
281 liquid_phase[MPL::PropertyType::density].template dValue<double>(
282 variables, MPL::Variable::liquid_phase_pressure, x_position, t,
283 dt);
284 auto const beta_LR = drho_LR_dp / rho_LR;
285
286 S_L = medium[MPL::PropertyType::saturation].template value<double>(
287 variables, x_position, t, dt);
288 variables.liquid_saturation = S_L;
289 variables_prev.liquid_saturation = S_L_prev;
290
291 // tangent derivative for Jacobian
292 double const dS_L_dp_cap =
293 medium[MPL::PropertyType::saturation].template dValue<double>(
294 variables, MPL::Variable::capillary_pressure, x_position, t,
295 dt);
296 // secant derivative from time discretization for storage
297 // use tangent, if secant is not available
298 double const DeltaS_L_Deltap_cap =
299 (p_cap_ip == p_cap_prev_ip)
300 ? dS_L_dp_cap
301 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
302
303 auto chi_S_L = S_L;
304 auto chi_S_L_prev = S_L_prev;
305 auto dchi_dS_L = 1.0;
306 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
307 {
308 auto const chi = [&medium, x_position, t, dt](double const S_L)
309 {
310 MPL::VariableArray variables;
311 variables.liquid_saturation = S_L;
312 return medium[MPL::PropertyType::bishops_effective_stress]
313 .template value<double>(variables, x_position, t, dt);
314 };
315 chi_S_L = chi(S_L);
316 chi_S_L_prev = chi(S_L_prev);
317
318 dchi_dS_L = medium[MPL::PropertyType::bishops_effective_stress]
319 .template dValue<double>(
320 variables, MPL::Variable::liquid_saturation,
321 x_position, t, dt);
322 }
323 // TODO (buchwaldj)
324 // should solid_grain_pressure or effective_pore_pressure remain?
325 // double const p_FR = -chi_S_L * p_cap_ip;
326 // variables.solid_grain_pressure = p_FR;
327
328 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
329 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
330
331 auto& phi = _ip_data[ip].porosity;
332 { // Porosity update
333
334 variables_prev.porosity = _ip_data[ip].porosity_prev;
335 phi = medium[MPL::PropertyType::porosity].template value<double>(
336 variables, variables_prev, x_position, t, dt);
337 variables.porosity = phi;
338 }
339
340 if (alpha < phi)
341 {
342 OGS_FATAL(
343 "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
344 "porosity {} in element/integration point {}/{}.",
345 alpha, phi, _element.getID(), ip);
346 }
347
348 double const k_rel =
349 medium[MPL::PropertyType::relative_permeability]
350 .template value<double>(variables, x_position, t, dt);
351 auto const mu =
352 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
353 variables, x_position, t, dt);
354
355 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
356 medium[MPL::PropertyType::permeability].value(variables, x_position,
357 t, dt));
358
359 GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
360 GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
361
362 auto const K_pT_thermal_osmosis =
363 (solid_phase.hasProperty(
364 MaterialPropertyLib::PropertyType::thermal_osmosis_coefficient)
365 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
366 solid_phase
367 [MPL::PropertyType::thermal_osmosis_coefficient]
368 .value(variables, x_position, t, dt))
369 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
370
371 // Consider anisotropic thermal expansion.
372 // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
373 // component.
374 Eigen::Matrix<double, 3, 3> const
375 solid_linear_thermal_expansion_coefficient =
376 MaterialPropertyLib::formEigenTensor<3>(
377 solid_phase
378 [MaterialPropertyLib::PropertyType::thermal_expansivity]
379 .value(variables, x_position, t, dt));
380
381 auto const rho_SR =
382 solid_phase[MPL::PropertyType::density].template value<double>(
383 variables, x_position, t, dt);
384
385 //
386 // pressure equation, pressure part.
387 //
388 laplace_p.noalias() +=
389 dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
390 laplace_T.noalias() +=
391 dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
392 const double alphaB_minus_phi = alpha - phi;
393 double const a0 = alphaB_minus_phi * beta_SR;
394 double const specific_storage_a_p =
395 S_L * (phi * beta_LR + S_L * a0 +
396 chi_S_L * alpha * alpha *
397 solid_elasticity.storageContribution(
398 solid_phase, variables, x_position, t, dt));
399 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
400
401 double const dspecific_storage_a_p_dp_cap =
402 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0 +
403 alpha * alpha *
404 solid_elasticity.storageContribution(
405 solid_phase, variables, x_position, t, dt) *
406 (chi_S_L + dchi_dS_L * S_L));
407 double const dspecific_storage_a_S_dp_cap =
408 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
409
410 storage_p_a_p.noalias() +=
411 N.transpose() * rho_LR * specific_storage_a_p * N * w;
412
413 storage_p_a_S.noalias() -= N.transpose() * rho_LR *
414 specific_storage_a_S * DeltaS_L_Deltap_cap *
415 N * w;
416
417 local_Jac
418 .template block<pressure_size, pressure_size>(pressure_index,
419 pressure_index)
420 .noalias() += N.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
421 rho_LR * dspecific_storage_a_p_dp_cap * N * w;
422
423 storage_p_a_S_Jpp.noalias() -=
424 N.transpose() * rho_LR *
425 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
426 specific_storage_a_S * dS_L_dp_cap) /
427 dt * N * w;
428
429 double const dk_rel_dS_L =
430 medium[MPL::PropertyType::relative_permeability]
431 .template dValue<double>(variables,
432 MPL::Variable::liquid_saturation,
433 x_position, t, dt);
434 GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
435 local_Jac
436 .template block<pressure_size, pressure_size>(pressure_index,
437 pressure_index)
438 .noalias() += dNdx.transpose() * rho_Ki_over_mu * grad_p_cap *
439 dk_rel_dS_L * dS_L_dp_cap * N * w;
440
441 local_Jac
442 .template block<pressure_size, pressure_size>(pressure_index,
443 pressure_index)
444 .noalias() += dNdx.transpose() * rho_LR * rho_Ki_over_mu * b *
445 dk_rel_dS_L * dS_L_dp_cap * N * w;
446
447 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
448 dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
449
450 //
451 // pressure equation, temperature part.
452 //
453 double const fluid_volumetric_thermal_expansion_coefficient =
454 MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
455 x_position, t, dt);
456 const double eff_thermal_expansion =
457 S_L * (alphaB_minus_phi *
458 solid_linear_thermal_expansion_coefficient.trace() +
459 phi * fluid_volumetric_thermal_expansion_coefficient +
460 alpha * solid_elasticity.thermalExpansivityContribution(
461 solid_linear_thermal_expansion_coefficient,
462 solid_phase, variables, x_position, t, dt));
463 M_pT.noalias() -=
464 N.transpose() * rho_LR * eff_thermal_expansion * N * w;
465
466 //
467 // temperature equation.
468 //
469 {
470 auto const specific_heat_capacity_fluid =
471 liquid_phase[MaterialPropertyLib::specific_heat_capacity]
472 .template value<double>(variables, x_position, t, dt);
473
474 auto const specific_heat_capacity_solid =
475 solid_phase
476 [MaterialPropertyLib::PropertyType::specific_heat_capacity]
477 .template value<double>(variables, x_position, t, dt);
478
479 M_TT.noalias() +=
480 w *
481 (rho_SR * specific_heat_capacity_solid * (1 - phi) +
482 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
483 N.transpose() * N;
484
485 auto const thermal_conductivity =
486 MaterialPropertyLib::formEigenTensor<GlobalDim>(
487 medium[MaterialPropertyLib::PropertyType::
488 thermal_conductivity]
489 .value(variables, x_position, t, dt));
490
491 GlobalDimVectorType const velocity_L = GlobalDimVectorType(
492 -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
493 K_pT_thermal_osmosis * dNdx * T);
494
495 K_TT.noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
496 N.transpose() * velocity_L.transpose() * dNdx *
497 rho_LR * specific_heat_capacity_fluid) *
498 w;
499
500 //
501 // temperature equation, pressure part
502 //
503 K_Tp.noalias() +=
504 dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
505 dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
506 N.transpose() * (dNdx * T).transpose() *
507 k_rel * Ki_over_mu * dNdx * w;
508
509 dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
510 N.transpose() * velocity_L.dot(dNdx * T) /
511 k_rel * dk_rel_dS_L * dS_L_dp_cap * N * w;
512 }
513 if (gas_phase && S_L < 1.0)
514 {
515 variables.density = rho_LR;
516
517 double const rho_wv =
518 gas_phase->property(MPL::PropertyType::density)
519 .template value<double>(variables, x_position, t, dt);
520
521 double const drho_wv_dT =
522 gas_phase->property(MPL::PropertyType::density)
523 .template dValue<double>(variables,
524 MPL::Variable::temperature,
525 x_position, t, dt);
526 double const drho_wv_dp =
527 gas_phase->property(MPL::PropertyType::density)
528 .template dValue<double>(
529 variables, MPL::Variable::liquid_phase_pressure,
530 x_position, t, dt);
531 auto const f_Tv =
532 gas_phase
533 ->property(
534 MPL::PropertyType::thermal_diffusion_enhancement_factor)
535 .template value<double>(variables, x_position, t, dt);
536
537 variables.porosity = phi;
538 auto const tortuosity =
539 medium.property(MPL::PropertyType::tortuosity)
540 .template value<double>(variables, x_position, t, dt);
541 double const D_v =
542 phi * (1.0 - S_L) * tortuosity *
543 gas_phase->property(MPL::PropertyType::diffusion)
544 .template value<double>(variables, x_position, t, dt);
545
546 double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
547 double const D_pv = D_v * drho_wv_dp;
548
549 GlobalDimVectorType const grad_T = dNdx * T;
550 GlobalDimVectorType const vapour_flux =
551 -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
552 double const specific_heat_capacity_vapour =
553 gas_phase->property(MaterialPropertyLib::specific_heat_capacity)
554 .template value<double>(variables, x_position, t, dt);
555
556 M_TT.noalias() +=
557 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
558 N.transpose() * N;
559
560 K_TT.noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
561 specific_heat_capacity_vapour * w;
562
563 double const storage_coefficient_by_water_vapor =
564 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
565
566 storage_p_a_p.noalias() +=
567 N.transpose() * storage_coefficient_by_water_vapor * N * w;
568
569 double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
570 M_pT.noalias() += N.transpose() * vapor_expansion_factor * N * w;
571
572 local_Jac
573 .template block<pressure_size, temperature_size>(
574 pressure_index, temperature_index)
575 .noalias() += dNdx.transpose() * f_Tv_D_Tv * dNdx * w;
576
577 local_rhs.template segment<pressure_size>(pressure_index)
578 .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
579
580 laplace_p.noalias() += dNdx.transpose() * D_pv * dNdx * w;
581
582 //
583 // Latent heat term
584 //
585 if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
586 {
587 double const factor = phi * (1 - S_L) / rho_LR;
588 // The volumetric latent heat of vaporization of liquid water
589 double const L0 =
590 gas_phase->property(MPL::PropertyType::specific_latent_heat)
591 .template value<double>(variables, x_position, t, dt) *
592 rho_LR;
593
594 double const drho_LR_dT =
595 liquid_phase.property(MPL::PropertyType::density)
596 .template dValue<double>(variables,
597 MPL::Variable::temperature,
598 x_position, t, dt);
599
600 double const rho_wv_over_rho_L = rho_wv / rho_LR;
601 M_TT.noalias() +=
602 factor * L0 *
603 (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
604 N.transpose() * N * w;
605
606 M_Tp.noalias() +=
607 (factor * L0 *
608 (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
609 L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
610 N.transpose() * N * w;
611
612 // temperature equation, temperature part
613 K_TT.noalias() +=
614 L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
615 // temperature equation, pressure part
616 K_Tp.noalias() +=
617 L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
618 }
619 }
620 }
621
622 if (_process_data.apply_mass_lumping)
623 {
624 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
625 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
626 storage_p_a_S_Jpp =
627 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
628 }
629
630 //
631 // -- Jacobian
632 //
633 // temperature equation.
634 local_Jac
635 .template block<temperature_size, temperature_size>(temperature_index,
636 temperature_index)
637 .noalias() += M_TT / dt + K_TT;
638 // temperature equation, pressure part
639 local_Jac
640 .template block<temperature_size, pressure_size>(temperature_index,
641 pressure_index)
642 .noalias() += K_Tp + dK_TT_dp;
643
644 // pressure equation, pressure part.
645 local_Jac
646 .template block<pressure_size, pressure_size>(pressure_index,
647 pressure_index)
648 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
649
650 // pressure equation, temperature part (contributed by thermal expansion).
651 local_Jac
652 .template block<pressure_size, temperature_size>(pressure_index,
653 temperature_index)
654 .noalias() += M_pT / dt + laplace_T;
655
656 //
657 // -- Residual
658 //
659 // temperature equation
660 local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
661 M_TT * (T - T_prev) / dt + K_TT * T;
662 local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
663 K_Tp * p_L;
664
665 // pressure equation
666 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
667 laplace_p * p_L + laplace_T * T +
668 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
669 M_pT * (T - T_prev) / dt;
670 if (gas_phase)
671 {
672 if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
673 {
674 // Jacobian: temperature equation, pressure part
675 local_Jac
676 .template block<temperature_size, pressure_size>(
677 temperature_index, pressure_index)
678 .noalias() += M_Tp / dt;
679 // RHS: temperature part
680 local_rhs.template segment<temperature_size>(temperature_index)
681 .noalias() -= M_Tp * (p_L - p_L_prev) / dt;
682 }
683 }
684}
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition ShapeMatrixPolicy.h:74

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::VariableArray::capillary_pressure, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::diffusion, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::formEigenTensor< 3 >(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getLiquidThermalExpansivity(), MaterialPropertyLib::VariableArray::grain_compressibility, MaterialPropertyLib::Phase::hasProperty(), NumLib::interpolateCoordinates(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::liquid_saturation, MaterialPropertyLib::VariableArray::liquid_saturation, OGS_FATAL, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, MaterialPropertyLib::Phase::property(), MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::specific_heat_capacity, MaterialPropertyLib::specific_latent_heat, MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, temperature_index, temperature_size, MaterialPropertyLib::thermal_diffusion_enhancement_factor, MaterialPropertyLib::thermal_expansivity, MaterialPropertyLib::thermal_osmosis_coefficient, MaterialPropertyLib::tortuosity, and MaterialPropertyLib::viscosity.

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::computeSecondaryVariableConcrete ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1197 of file ThermoRichardsFlowFEM-impl.h.

1201{
1202 auto const T =
1203 local_x.template segment<temperature_size>(temperature_index);
1204
1205 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
1206
1207 auto p_L_prev =
1208 local_x_prev.template segment<pressure_size>(pressure_index);
1209
1210 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
1211 auto const& liquid_phase = medium.phase("AqueousLiquid");
1212 auto const& solid_phase = medium.phase("Solid");
1213 MPL::VariableArray variables;
1214 MPL::VariableArray variables_prev;
1215
1216 unsigned const n_integration_points =
1217 _integration_method.getNumberOfPoints();
1218
1219 double saturation_avg = 0;
1220 double porosity_avg = 0;
1221
1222 for (unsigned ip = 0; ip < n_integration_points; ip++)
1223 {
1224 auto const& N = _ip_data[ip].N;
1225
1226 ParameterLib::SpatialPosition const x_position{
1227 std::nullopt, _element.getID(),
1228 MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
1229 ShapeMatricesType>(
1230 _element, N))};
1231
1232 double T_ip;
1233 NumLib::shapeFunctionInterpolate(T, N, T_ip);
1234
1235 double p_cap_ip;
1236 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
1237
1238 double p_cap_prev_ip;
1239 NumLib::shapeFunctionInterpolate(-p_L_prev, N, p_cap_prev_ip);
1240
1241 variables.capillary_pressure = p_cap_ip;
1242 variables.liquid_phase_pressure = -p_cap_ip;
1243 // setting pG to 1 atm
1244 // TODO : rewrite equations s.t. p_L = pG-p_cap
1245 variables.gas_phase_pressure = 1.0e5;
1246
1247 variables.temperature = T_ip;
1248
1249 auto& S_L = _ip_data[ip].saturation;
1250 auto const S_L_prev = _ip_data[ip].saturation_prev;
1251 S_L = medium[MPL::PropertyType::saturation].template value<double>(
1252 variables, x_position, t, dt);
1253 variables.liquid_saturation = S_L;
1254 variables_prev.liquid_saturation = S_L_prev;
1255
1256 auto chi_S_L = S_L;
1257 auto chi_S_L_prev = S_L_prev;
1258 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
1259 {
1260 auto const chi = [&medium, x_position, t, dt](double const S_L)
1261 {
1262 MPL::VariableArray variables;
1263 variables.liquid_saturation = S_L;
1264 return medium[MPL::PropertyType::bishops_effective_stress]
1265 .template value<double>(variables, x_position, t, dt);
1266 };
1267 chi_S_L = chi(S_L);
1268 chi_S_L_prev = chi(S_L_prev);
1269 }
1270 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
1271 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
1272
1273 auto const alpha =
1274 medium[MPL::PropertyType::biot_coefficient].template value<double>(
1275 variables, x_position, t, dt);
1276
1277 auto& solid_elasticity = *_process_data.simplified_elasticity;
1278 auto const beta_S =
1279 solid_elasticity.bulkCompressibilityFromYoungsModulus(
1280 solid_phase, variables, x_position, t, dt);
1281 auto const beta_SR = (1 - alpha) * beta_S;
1282 variables.grain_compressibility = beta_SR;
1283
1284 auto& phi = _ip_data[ip].porosity;
1285 { // Porosity update
1286 variables_prev.porosity = _ip_data[ip].porosity_prev;
1287 phi = medium[MPL::PropertyType::porosity].template value<double>(
1288 variables, variables_prev, x_position, t, dt);
1289 variables.porosity = phi;
1290 }
1291
1292 auto const mu =
1293 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
1294 variables, x_position, t, dt);
1295 auto const rho_LR =
1296 liquid_phase[MPL::PropertyType::density].template value<double>(
1297 variables, x_position, t, dt);
1298
1299 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
1300 medium[MPL::PropertyType::permeability].value(variables, x_position,
1301 t, dt));
1302
1303 double const k_rel =
1304 medium[MPL::PropertyType::relative_permeability]
1305 .template value<double>(variables, x_position, t, dt);
1306
1307 GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1308
1309 auto const rho_SR =
1310 solid_phase[MPL::PropertyType::density].template value<double>(
1311 variables, x_position, t, dt);
1312 _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
1313
1314 auto const& b = _process_data.specific_body_force;
1315
1316 auto const K_pT_thermal_osmosis =
1317 (solid_phase.hasProperty(
1318 MaterialPropertyLib::PropertyType::thermal_osmosis_coefficient)
1319 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
1320 solid_phase
1321 [MPL::PropertyType::thermal_osmosis_coefficient]
1322 .value(variables, x_position, t, dt))
1323 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
1324
1325 // Compute the velocity
1326 auto const& dNdx = _ip_data[ip].dNdx;
1327 _ip_data[ip].v_darcy.noalias() = -K_over_mu * dNdx * p_L -
1328 K_pT_thermal_osmosis * dNdx * T +
1329 rho_LR * K_over_mu * b;
1330
1331 saturation_avg += S_L;
1332 porosity_avg += phi;
1333 }
1334 saturation_avg /= n_integration_points;
1335 porosity_avg /= n_integration_points;
1336
1337 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
1338 (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
1339}

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::temperature, temperature_index, MaterialPropertyLib::thermal_osmosis_coefficient, and MaterialPropertyLib::viscosity.

◆ getIntPtDarcyVelocity()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1118 of file ThermoRichardsFlowFEM-impl.h.

1124{
1125 unsigned const n_integration_points =
1126 _integration_method.getNumberOfPoints();
1127
1128 cache.clear();
1129 auto cache_matrix = MathLib::createZeroedMatrix<
1130 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1131 cache, GlobalDim, n_integration_points);
1132
1133 for (unsigned ip = 0; ip < n_integration_points; ip++)
1134 {
1135 cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
1136 }
1137
1138 return cache;
1139}

References _integration_method, _ip_data, and MathLib::createZeroedMatrix().

◆ getIntPtDryDensitySolid()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtDryDensitySolid ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1185 of file ThermoRichardsFlowFEM-impl.h.

1191{
1192 return ProcessLib::getIntegrationPointScalarData(
1193 _ip_data, &IpData::dry_density_solid, cache);
1194}
ProcessLib::getIntegrationPointScalarData
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:176

References _ip_data, ProcessLib::ThermoRichardsFlow::IntegrationPointData< ShapeMatricesType >::dry_density_solid, and ProcessLib::getIntegrationPointScalarData().

◆ getIntPtPorosity()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtPorosity ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1173 of file ThermoRichardsFlowFEM-impl.h.

1178{
1179 return ProcessLib::getIntegrationPointScalarData(_ip_data,
1180 &IpData::porosity, cache);
1181}

References _ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::ThermoRichardsFlow::IntegrationPointData< ShapeMatricesType >::porosity.

Referenced by getPorosity().

◆ getIntPtSaturation()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtSaturation ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1152 of file ThermoRichardsFlowFEM-impl.h.

1157{
1158 return ProcessLib::getIntegrationPointScalarData(
1159 _ip_data, &IpData::saturation, cache);
1160}

References _ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::ThermoRichardsFlow::IntegrationPointData< ShapeMatricesType >::saturation.

Referenced by getSaturation().

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunction, int GlobalDim>
unsigned ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getNumberOfIntegrationPoints ( ) const
overrideprivatevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1343 of file ThermoRichardsFlowFEM-impl.h.

1344{
1345 return _integration_method.getNumberOfPoints();
1346}

References _integration_method, and getNumberOfIntegrationPoints().

Referenced by getNumberOfIntegrationPoints().

◆ getPorosity()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getPorosity ( ) const
overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1164 of file ThermoRichardsFlowFEM-impl.h.

1165{
1166 std::vector<double> result;
1167 getIntPtPorosity(0, {}, {}, result);
1168 return result;
1169}
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler::getIntPtPorosity
std::vector< double > const & getIntPtPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Definition ThermoRichardsFlowFEM-impl.h:1173

References getIntPtPorosity().

◆ getSaturation()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getSaturation ( ) const
overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1143 of file ThermoRichardsFlowFEM-impl.h.

1144{
1145 std::vector<double> result;
1146 getIntPtSaturation(0, {}, {}, result);
1147 return result;
1148}
ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler::getIntPtSaturation
std::vector< double > const & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Definition ThermoRichardsFlowFEM-impl.h:1152

References getIntPtSaturation(), and getSaturation().

Referenced by getSaturation().

◆ getShapeMatrix()

template<typename ShapeFunction, int GlobalDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getShapeMatrix ( const unsigned integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 113 of file ThermoRichardsFlowFEM.h.

115 {
116 auto const& N = _ip_data[integration_point].N;
117
118 // assumes N is stored contiguously in memory
119 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
120 }

References _ip_data.

◆ initializeConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::initializeConcrete ( )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 83 of file ThermoRichardsFlowFEM.h.

84 {
85 unsigned const n_integration_points =
86 _integration_method.getNumberOfPoints();
87
88 for (unsigned ip = 0; ip < n_integration_points; ip++)
89 {
90 auto& ip_data = _ip_data[ip];
91 ip_data.pushBackState();
92 }
93 }

References _integration_method, and _ip_data.

◆ postTimestepConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::postTimestepConcrete ( Eigen::VectorXd const & ,
Eigen::VectorXd const & ,
double const ,
double const ,
int const  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 95 of file ThermoRichardsFlowFEM.h.

99 {
100 unsigned const n_integration_points =
101 _integration_method.getNumberOfPoints();
102
103 for (unsigned ip = 0; ip < n_integration_points; ip++)
104 {
105 _ip_data[ip].pushBackState();
106 }
107 }

References _integration_method, and _ip_data.

◆ setInitialConditionsConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::setInitialConditionsConcrete ( Eigen::VectorXd const local_x,
double const t,
int const process_id )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 111 of file ThermoRichardsFlowFEM-impl.h.

115{
116 assert(local_x.size() == temperature_size + pressure_size);
117
118 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
119
120 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
121 MPL::VariableArray variables;
122
123 unsigned const n_integration_points =
124 _integration_method.getNumberOfPoints();
125 for (unsigned ip = 0; ip < n_integration_points; ip++)
126 {
127 auto const& N = _ip_data[ip].N;
128
129 ParameterLib::SpatialPosition const x_position{
130 std::nullopt, _element.getID(),
131 MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
132 ShapeMatricesType>(
133 _element, N))};
134
135 double p_cap_ip;
136 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
137
138 variables.capillary_pressure = p_cap_ip;
139 variables.liquid_phase_pressure = -p_cap_ip;
140 // setting pG to 1 atm
141 // TODO : rewrite equations s.t. p_L = pG-p_cap
142 variables.gas_phase_pressure = 1.0e5;
143
144 // Note: temperature dependent saturation model is not considered so
145 // far.
146 _ip_data[ip].saturation_prev =
147 medium[MPL::PropertyType::saturation].template value<double>(
148 variables, x_position, t,
149 std::numeric_limits<double>::quiet_NaN());
150 }
151}

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, pressure_index, pressure_size, MaterialPropertyLib::saturation, NumLib::detail::shapeFunctionInterpolate(), and temperature_size.

◆ setIPDataInitialConditions()

template<typename ShapeFunction, int GlobalDim>
std::size_t ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::setIPDataInitialConditions ( std::string_view const name,
double const * values,
int const integration_order )
overridevirtual
Returns
the number of read integration points.

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 82 of file ThermoRichardsFlowFEM-impl.h.

86{
87 if (integration_order !=
88 static_cast<int>(_integration_method.getIntegrationOrder()))
89 {
90 OGS_FATAL(
91 "Setting integration point initial conditions; The integration "
92 "order of the local assembler for element {:d} is different "
93 "from the integration order in the initial condition.",
94 _element.getID());
95 }
96
97 if (name == "saturation")
98 {
99 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
100 &IpData::saturation);
101 }
102 if (name == "porosity")
103 {
104 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
105 &IpData::porosity);
106 }
107 return 0;
108}
ProcessLib::setIntegrationPointScalarData
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
Definition SetOrGetIntegrationPointData.h:215

References _element, _integration_method, _ip_data, OGS_FATAL, ProcessLib::ThermoRichardsFlow::IntegrationPointData< ShapeMatricesType >::porosity, ProcessLib::ThermoRichardsFlow::IntegrationPointData< ShapeMatricesType >::saturation, and ProcessLib::setIntegrationPointScalarData().

Member Data Documentation

◆ _element

template<typename ShapeFunction, int GlobalDim>
MeshLib::Element const& ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_element
private

Definition at line 157 of file ThermoRichardsFlowFEM.h.

Referenced by ThermoRichardsFlowLocalAssembler(), assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), setInitialConditionsConcrete(), and setIPDataInitialConditions().

◆ _integration_method

template<typename ShapeFunction, int GlobalDim>
NumLib::GenericIntegrationMethod const& ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_integration_method
private

Definition at line 156 of file ThermoRichardsFlowFEM.h.

Referenced by ThermoRichardsFlowLocalAssembler(), assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), getIntPtDarcyVelocity(), getNumberOfIntegrationPoints(), initializeConcrete(), postTimestepConcrete(), setInitialConditionsConcrete(), and setIPDataInitialConditions().

◆ _ip_data

template<typename ShapeFunction, int GlobalDim>
std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_ip_data
private

Definition at line 154 of file ThermoRichardsFlowFEM.h.

Referenced by ThermoRichardsFlowLocalAssembler(), assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), getIntPtDarcyVelocity(), getIntPtDryDensitySolid(), getIntPtPorosity(), getIntPtSaturation(), getShapeMatrix(), initializeConcrete(), postTimestepConcrete(), setInitialConditionsConcrete(), and setIPDataInitialConditions().

◆ _is_axially_symmetric

template<typename ShapeFunction, int GlobalDim>
bool const ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_is_axially_symmetric
private

Definition at line 158 of file ThermoRichardsFlowFEM.h.

Referenced by ThermoRichardsFlowLocalAssembler().

◆ _process_data

template<typename ShapeFunction, int GlobalDim>
ThermoRichardsFlowProcessData& ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_process_data
private

Definition at line 152 of file ThermoRichardsFlowFEM.h.

Referenced by ThermoRichardsFlowLocalAssembler(), assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), and setInitialConditionsConcrete().

◆ pressure_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::pressure_index = temperature_size
staticprivate

Definition at line 162 of file ThermoRichardsFlowFEM.h.

Referenced by assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), and setInitialConditionsConcrete().

◆ pressure_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::pressure_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 163 of file ThermoRichardsFlowFEM.h.

Referenced by assemble(), assembleWithJacobian(), and setInitialConditionsConcrete().

◆ temperature_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::temperature_index = 0
staticprivate

Definition at line 160 of file ThermoRichardsFlowFEM.h.

Referenced by assemble(), assembleWithJacobian(), and computeSecondaryVariableConcrete().

◆ temperature_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::temperature_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 161 of file ThermoRichardsFlowFEM.h.

Referenced by assemble(), assembleWithJacobian(), and setInitialConditionsConcrete().

The documentation for this class was generated from the following files: