[PL/TH2M] Extract data and model of W1 equation

ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h

@@ -41,6 +41,7 @@
 #include "TotalStress.h"
 #include "VapourPartialPressure.h"
 #include "Viscosity.h"
+#include "WEquation.h"
 
 namespace ProcessLib::TH2M
 {
@@ -166,6 +167,8 @@ struct ConstitutiveTempData
     FC4MCTData fC_4_MCT;
     FC4MCuData fC_4_MCu;
 
+    FW1Data<DisplacementDim> fW_1;
+
     using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
     using DisplacementDimMatrix =
         Eigen::Matrix<double, DisplacementDim, DisplacementDim>;

ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h

@@ -31,6 +31,7 @@
 #include "ThermalConductivity.h"
 #include "TotalStress.h"
 #include "Viscosity.h"
+#include "WEquation.h"
 
 namespace ProcessLib::TH2M
 {
@@ -91,6 +92,8 @@ struct ConstitutiveModels
     FC4MCpCModel fC_4_MCpC_model;
     FC4MCTModel<DisplacementDim> fC_4_MCT_model;
     FC4MCuModel fC_4_MCu_model;
+
+    FW1Model<DisplacementDim> fW_1_model;
 };
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M

ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "WEquation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void FW1Model<DisplacementDim>::eval(
+    AdvectionData<DisplacementDim> const& advection_data,
+    FluidDensityData const& fluid_density_data,
+    FW1Data<DisplacementDim>& fW_1) const
+{
+    fW_1.A = advection_data.advection_W_G * fluid_density_data.rho_GR +
+             advection_data.advection_W_L * fluid_density_data.rho_LR;
+}
+
+template struct FW1Model<2>;
+template struct FW1Model<3>;
+}  // namespace ConstitutiveRelations
+}  // namespace ProcessLib::TH2M

ProcessLib/TH2M/ConstitutiveRelations/WEquation.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include "Advection.h"
+#include "Base.h"
+#include "Biot.h"
+#include "ConstitutiveDensity.h"
+#include "FluidDensity.h"
+#include "PermeabilityData.h"
+#include "PhaseTransitionData.h"
+#include "Porosity.h"
+#include "Saturation.h"
+#include "SolidCompressibility.h"
+#include "Viscosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct FW1Data
+{
+    GlobalDimMatrix<DisplacementDim> A;
+};
+
+template <int DisplacementDim>
+struct FW1Model
+{
+    void eval(AdvectionData<DisplacementDim> const& advection_data,
+              FluidDensityData const& fluid_density_data,
+              FW1Data<DisplacementDim>& fW_1) const;
+};
+
+extern template struct FW1Model<2>;
+extern template struct FW1Model<3>;
+}  // namespace ConstitutiveRelations
+}  // namespace ProcessLib::TH2M

ProcessLib/TH2M/TH2MFEM-impl.h

@@ -351,6 +351,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                                    current_state.S_L_data,
                                    ip_cv.fC_4_MCu);
 
+        models.fW_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
+                               ip_cv.fW_1);
+
         // for variable output
         auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
 
@@ -1271,10 +1274,7 @@ void TH2MLocalAssembler<
 
         LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
 
-        fW.noalias() += gradNpT *
-                        (ip_cv.advection_data.advection_W_G * rhoGR +
-                         ip_cv.advection_data.advection_W_L * rhoLR) *
-                        b * w;
+        fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
 
         if (!this->process_data_.apply_mass_lumping)
         {
@@ -1811,10 +1811,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
 
         // fW_1
-        fW.noalias() += gradNpT *
-                        (ip_cv.advection_data.advection_W_G * rhoGR +
-                         ip_cv.advection_data.advection_W_L * rhoLR) *
-                        b * w;
+        fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
 
         // fW_2 = \int (f - g) * s_L_dot
         if (!this->process_data_.apply_mass_lumping)