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 29 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 int getNumberOfVectorElementsForDeformation () 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 36 of file ThermoRichardsFlowFEM.h.

◆ GlobalDimVectorType

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

Definition at line 37 of file ThermoRichardsFlowFEM.h.

◆ IpData

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

Definition at line 39 of file ThermoRichardsFlowFEM.h.

◆ ShapeMatricesType

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

Definition at line 34 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 30 of file ThermoRichardsFlowFEM-impl.h.

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

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 680 of file ThermoRichardsFlowFEM-impl.h.

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, 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::Gas, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getLiquidThermalExpansivity(), ProcessLib::getThermoOsmoticCoefficient(), 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::Solid, 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::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 150 of file ThermoRichardsFlowFEM-impl.h.

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, 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::Gas, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getLiquidThermalExpansivity(), ProcessLib::getThermoOsmoticCoefficient(), 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::Solid, 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::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 1187 of file ThermoRichardsFlowFEM-impl.h.

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, 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, ProcessLib::getThermoOsmoticCoefficient(), 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::Solid, MaterialPropertyLib::VariableArray::temperature, temperature_index, 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 1108 of file ThermoRichardsFlowFEM-impl.h.

1114{
1115 unsigned const n_integration_points =
1116 _integration_method.getNumberOfPoints();
1117
1118 cache.clear();
1119 auto cache_matrix = MathLib::createZeroedMatrix<
1120 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1121 cache, GlobalDim, n_integration_points);
1122
1123 for (unsigned ip = 0; ip < n_integration_points; ip++)
1124 {
1125 cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
1126 }
1127
1128 return cache;
1129}

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 1175 of file ThermoRichardsFlowFEM-impl.h.

1181{
1182 return ProcessLib::getIntegrationPointScalarData(
1183 _ip_data, &IpData::dry_density_solid, cache);
1184}
ProcessLib::getIntegrationPointScalarData
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:168

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 1163 of file ThermoRichardsFlowFEM-impl.h.

1168{
1169 return ProcessLib::getIntegrationPointScalarData(_ip_data,
1170 &IpData::porosity, cache);
1171}

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 1142 of file ThermoRichardsFlowFEM-impl.h.

1147{
1148 return ProcessLib::getIntegrationPointScalarData(
1149 _ip_data, &IpData::saturation, cache);
1150}

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 1330 of file ThermoRichardsFlowFEM-impl.h.

1331{
1332 return _integration_method.getNumberOfPoints();
1333}

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 1154 of file ThermoRichardsFlowFEM-impl.h.

1155{
1156 std::vector<double> result;
1157 getIntPtPorosity(0, {}, {}, result);
1158 return result;
1159}
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:1163

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 1133 of file ThermoRichardsFlowFEM-impl.h.

1134{
1135 std::vector<double> result;
1136 getIntPtSaturation(0, {}, {}, result);
1137 return result;
1138}
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:1142

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 106 of file ThermoRichardsFlowFEM.h.

108 {
109 auto const& N = _ip_data[integration_point].N;
110
111 // assumes N is stored contiguously in memory
112 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
113 }

References _ip_data.

◆ initializeConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 76 of file ThermoRichardsFlowFEM.h.

77 {
78 unsigned const n_integration_points =
79 _integration_method.getNumberOfPoints();
80
81 for (unsigned ip = 0; ip < n_integration_points; ip++)
82 {
83 auto& ip_data = _ip_data[ip];
84 ip_data.pushBackState();
85 }
86 }

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 88 of file ThermoRichardsFlowFEM.h.

92 {
93 unsigned const n_integration_points =
94 _integration_method.getNumberOfPoints();
95
96 for (unsigned ip = 0; ip < n_integration_points; ip++)
97 {
98 _ip_data[ip].pushBackState();
99 }
100 }

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 107 of file ThermoRichardsFlowFEM-impl.h.

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

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 78 of file ThermoRichardsFlowFEM-impl.h.

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

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 150 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 149 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 147 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 151 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 145 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 155 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 156 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 153 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 154 of file ThermoRichardsFlowFEM.h.

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

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