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

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

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

◆ assembleWithJacobian()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ computeSecondaryVariableConcrete()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

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

◆ getIntPtDarcyVelocity()

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

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

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

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

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

1185{
1186 return ProcessLib::getIntegrationPointScalarData(
1187 _ip_data, &IpData::dry_density_solid, cache);
1188}
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 1167 of file ThermoRichardsFlowFEM-impl.h.

1172{
1173 return ProcessLib::getIntegrationPointScalarData(_ip_data,
1174 &IpData::porosity, cache);
1175}

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

1151{
1152 return ProcessLib::getIntegrationPointScalarData(
1153 _ip_data, &IpData::saturation, cache);
1154}

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

1338{
1339 return _integration_method.getNumberOfPoints();
1340}

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

1159{
1160 std::vector<double> result;
1161 getIntPtPorosity(0, {}, {}, result);
1162 return result;
1163}
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:1167

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

1138{
1139 std::vector<double> result;
1140 getIntPtSaturation(0, {}, {}, result);
1141 return result;
1142}
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:1146

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

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

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

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