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

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

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

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(), 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::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 149 of file ThermoRichardsFlowFEM-impl.h.

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

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

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

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

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

1190{
1191 return ProcessLib::getIntegrationPointScalarData(
1192 _ip_data, &IpData::dry_density_solid, cache);
1193}
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 1172 of file ThermoRichardsFlowFEM-impl.h.

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

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

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

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

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

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

1164{
1165 std::vector<double> result;
1166 getIntPtPorosity(0, {}, {}, result);
1167 return result;
1168}
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:1172

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

1143{
1144 std::vector<double> result;
1145 getIntPtSaturation(0, {}, {}, result);
1146 return result;
1147}
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:1151

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

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

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

81{
82 if (integration_order !=
83 static_cast<int>(_integration_method.getIntegrationOrder()))
84 {
85 OGS_FATAL(
86 "Setting integration point initial conditions; The integration "
87 "order of the local assembler for element {:d} is different "
88 "from the integration order in the initial condition.",
89 _element.getID());
90 }
91
92 if (name == "saturation")
93 {
94 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
95 &IpData::saturation);
96 }
97 if (name == "porosity")
98 {
99 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
100 &IpData::porosity);
101 }
102 return 0;
103}
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: