From e50a0c7f653e94b8e25d0642b68f8dc5468b3f90 Mon Sep 17 00:00:00 2001
From: Boyan Meng <meng.boyan@ufz.de>
Date: Thu, 28 Apr 2022 19:47:14 +0200
Subject: [PATCH] [PL/TH2PP] Add a third component: PV, SV and local matrices
---
...CreateThermalTwoPhaseFlowWithPPProcess.cpp | 2 +
...malTwoPhaseFlowWithPPLocalAssembler-impl.h | 62 +++++++++++++++++-
.../ThermalTwoPhaseFlowWithPPLocalAssembler.h | 64 ++++++++++++++++++-
.../ThermalTwoPhaseFlowWithPPProcess.cpp | 18 ++++++
4 files changed, 142 insertions(+), 4 deletions(-)
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp
index 45aee43f5fe..52324bb9706 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/CreateThermalTwoPhaseFlowWithPPProcess.cpp
@@ -103,6 +103,8 @@ std::unique_ptr<Process> createThermalTwoPhaseFlowWithPPProcess(
"gas_pressure",
//! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_THERMAL__process_variables__capillary_pressure}
"capillary_pressure",
+ //! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_THERMAL__process_variables__total_molar_fraction_contaminant}
+ "total_molar_fraction_contaminant",
//! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_THERMAL__process_variables__temperature}
"temperature"});
std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
index 03d65fdf7a6..404bc834626 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -86,6 +86,9 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
nonwet_pressure_matrix_index, cap_pressure_matrix_index);
+ auto Mgx =
+ local_M.template block<nonwet_pressure_size, tot_mol_frac_contam_size>(
+ nonwet_pressure_matrix_index, tot_mol_frac_contam_matrix_index);
auto Mgt = local_M.template block<nonwet_pressure_size, temperature_size>(
nonwet_pressure_matrix_index, temperature_matrix_index);
@@ -93,13 +96,32 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
cap_pressure_matrix_index, nonwet_pressure_matrix_index);
auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
cap_pressure_matrix_index, cap_pressure_matrix_index);
+ auto Mlx =
+ local_M.template block<cap_pressure_size, tot_mol_frac_contam_size>(
+ cap_pressure_matrix_index, tot_mol_frac_contam_matrix_index);
auto Mlt = local_M.template block<cap_pressure_size, temperature_size>(
cap_pressure_matrix_index, temperature_matrix_index);
+ auto Mcp =
+ local_M.template block<tot_mol_frac_contam_size, nonwet_pressure_size>(
+ tot_mol_frac_contam_matrix_index, nonwet_pressure_matrix_index);
+ auto Mcpc =
+ local_M.template block<tot_mol_frac_contam_size, cap_pressure_size>(
+ tot_mol_frac_contam_matrix_index, cap_pressure_matrix_index);
+ auto Mcx = local_M.template block<tot_mol_frac_contam_size,
+ tot_mol_frac_contam_size>(
+ tot_mol_frac_contam_matrix_index, tot_mol_frac_contam_matrix_index);
+ auto Mct =
+ local_M.template block<tot_mol_frac_contam_size, temperature_size>(
+ tot_mol_frac_contam_matrix_index, temperature_matrix_index);
+
auto Mep = local_M.template block<temperature_size, nonwet_pressure_size>(
temperature_matrix_index, nonwet_pressure_matrix_index);
auto Mepc = local_M.template block<temperature_size, cap_pressure_size>(
temperature_matrix_index, cap_pressure_matrix_index);
+ auto Mex =
+ local_M.template block<temperature_size, tot_mol_frac_contam_size>(
+ temperature_matrix_index, tot_mol_frac_contam_matrix_index);
auto Met = local_M.template block<temperature_size, temperature_size>(
temperature_matrix_index, temperature_matrix_index);
@@ -111,6 +133,9 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
auto Kgpc = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
nonwet_pressure_matrix_index, cap_pressure_matrix_index);
+ auto Kgx =
+ local_K.template block<nonwet_pressure_size, tot_mol_frac_contam_size>(
+ nonwet_pressure_matrix_index, tot_mol_frac_contam_matrix_index);
auto Kgt = local_K.template block<nonwet_pressure_size, temperature_size>(
nonwet_pressure_matrix_index, temperature_matrix_index);
@@ -118,9 +143,25 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
cap_pressure_matrix_index, nonwet_pressure_matrix_index);
auto Klpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
cap_pressure_matrix_index, cap_pressure_matrix_index);
+ auto Klx =
+ local_K.template block<cap_pressure_size, tot_mol_frac_contam_size>(
+ cap_pressure_matrix_index, tot_mol_frac_contam_matrix_index);
auto Klt = local_K.template block<cap_pressure_size, temperature_size>(
cap_pressure_matrix_index, temperature_matrix_index);
+ auto Kcp =
+ local_K.template block<tot_mol_frac_contam_size, nonwet_pressure_size>(
+ tot_mol_frac_contam_matrix_index, nonwet_pressure_matrix_index);
+ auto Kcpc =
+ local_K.template block<tot_mol_frac_contam_size, cap_pressure_size>(
+ tot_mol_frac_contam_matrix_index, cap_pressure_matrix_index);
+ auto Kcx = local_K.template block<tot_mol_frac_contam_size,
+ tot_mol_frac_contam_size>(
+ tot_mol_frac_contam_matrix_index, tot_mol_frac_contam_matrix_index);
+ auto Kct =
+ local_K.template block<tot_mol_frac_contam_size, temperature_size>(
+ tot_mol_frac_contam_matrix_index, temperature_matrix_index);
+
auto Kep = local_K.template block<temperature_size, nonwet_pressure_size>(
temperature_matrix_index, nonwet_pressure_matrix_index);
auto Kepc = local_K.template block<temperature_size, cap_pressure_size>(
@@ -132,6 +173,8 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
nonwet_pressure_matrix_index);
auto Bl =
local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
+ auto Bc = local_b.template segment<tot_mol_frac_contam_size>(
+ tot_mol_frac_contam_matrix_index);
auto Be =
local_b.template segment<temperature_size>(temperature_matrix_index);
@@ -158,9 +201,10 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
auto const& diffusion_operator = ip_data.diffusion_operator;
double pg_int_pt = 0.;
double pc_int_pt = 0.;
- double T_int_pt = 0.0;
+ double Xc_int_pt = 0.;
+ double T_int_pt = 0.;
NumLib::shapeFunctionInterpolate(local_x, N, pg_int_pt, pc_int_pt,
- T_int_pt);
+ Xc_int_pt, T_int_pt);
double const ideal_gas_constant_times_T_int_pt =
IdealGasConstant * T_int_pt;
vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
@@ -557,16 +601,30 @@ void ThermalTwoPhaseFlowWithPPLocalAssembler<
Mgpc(row, column) = 0.0;
Mgt(row, row) += Mgt(row, column);
Mgt(row, column) = 0.0;
+ Mgx(row, row) += Mgx(row, column);
+ Mgx(row, column) = 0.0;
Mlp(row, row) += Mlp(row, column);
Mlp(row, column) = 0.0;
Mlpc(row, row) += Mlpc(row, column);
Mlpc(row, column) = 0.0;
+ Mlx(row, row) += Mlx(row, column);
+ Mlx(row, column) = 0.0;
Mlt(row, row) += Mlt(row, column);
Mlt(row, column) = 0.0;
+ Mcp(row, row) += Mcp(row, column);
+ Mcp(row, column) = 0.0;
+ Mcpc(row, row) += Mcpc(row, column);
+ Mcpc(row, column) = 0.0;
+ Mcx(row, row) += Mcx(row, column);
+ Mcx(row, column) = 0.0;
+ Mct(row, row) += Mct(row, column);
+ Mct(row, column) = 0.0;
Mep(row, row) += Mep(row, column);
Mep(row, column) = 0.0;
Mepc(row, row) += Mepc(row, column);
Mepc(row, column) = 0.0;
+ Mex(row, row) += Mex(row, column);
+ Mex(row, column) = 0.0;
Met(row, row) += Met(row, column);
Met(row, column) = 0.0;
}
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
index d934f0ceefb..2e5d7c0fe71 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
@@ -51,7 +51,7 @@ struct IntegrationPointData final
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
-const unsigned NUM_NODAL_DOF = 3;
+const unsigned NUM_NODAL_DOF = 4;
class ThermalTwoPhaseFlowWithPPLocalAssemblerInterface
: public ProcessLib::LocalAssemblerInterface,
@@ -69,6 +69,24 @@ public:
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtLiquidMolFracContaminant(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtGasMolFracWater(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtGasMolFracContaminant(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const = 0;
};
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
@@ -104,6 +122,12 @@ public:
_saturation(
std::vector<double>(_integration_method.getNumberOfPoints())),
_pressure_wetting(
+ std::vector<double>(_integration_method.getNumberOfPoints())),
+ _liquid_molar_fraction_contaminant(
+ std::vector<double>(_integration_method.getNumberOfPoints())),
+ _gas_molar_fraction_water(
+ std::vector<double>(_integration_method.getNumberOfPoints())),
+ _gas_molar_fraction_contaminant(
std::vector<double>(_integration_method.getNumberOfPoints()))
{
unsigned const n_integration_points =
@@ -164,6 +188,36 @@ public:
return _pressure_wetting;
}
+ std::vector<double> const& getIntPtLiquidMolFracContaminant(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(!_liquid_molar_fraction_contaminant.empty());
+ return _liquid_molar_fraction_contaminant;
+ }
+
+ std::vector<double> const& getIntPtGasMolFracWater(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(!_gas_molar_fraction_water.empty());
+ return _gas_molar_fraction_water;
+ }
+
+ std::vector<double> const& getIntPtGasMolFracContaminant(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(!_gas_molar_fraction_contaminant.empty());
+ return _gas_molar_fraction_contaminant;
+ }
+
private:
MeshLib::Element const& _element;
@@ -179,13 +233,19 @@ private:
std::vector<double> _saturation;
std::vector<double> _pressure_wetting;
+ std::vector<double> _liquid_molar_fraction_contaminant;
+ std::vector<double> _gas_molar_fraction_water;
+ std::vector<double> _gas_molar_fraction_contaminant;
static const int nonwet_pressure_matrix_index = 0;
static const int cap_pressure_matrix_index = ShapeFunction::NPOINTS;
- static const int temperature_matrix_index = 2 * ShapeFunction::NPOINTS;
+ static const int tot_mol_frac_contam_matrix_index =
+ 2 * ShapeFunction::NPOINTS;
+ static const int temperature_matrix_index = 3 * ShapeFunction::NPOINTS;
static const int nonwet_pressure_size = ShapeFunction::NPOINTS;
static const int cap_pressure_size = ShapeFunction::NPOINTS;
+ static const int tot_mol_frac_contam_size = ShapeFunction::NPOINTS;
static const int temperature_size = ShapeFunction::NPOINTS;
};
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
index 3625f50e574..0ee4920dbdb 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
@@ -64,6 +64,24 @@ void ThermalTwoPhaseFlowWithPPProcess::initializeConcreteProcess(
makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ThermalTwoPhaseFlowWithPPLocalAssemblerInterface::
getIntPtWettingPressure));
+
+ _secondary_variables.addSecondaryVariable(
+ "liquid_molar_fraction_contaminant",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &ThermalTwoPhaseFlowWithPPLocalAssemblerInterface::
+ getIntPtLiquidMolFracContaminant));
+
+ _secondary_variables.addSecondaryVariable(
+ "gas_molar_fraction_water",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &ThermalTwoPhaseFlowWithPPLocalAssemblerInterface::
+ getIntPtGasMolFracWater));
+
+ _secondary_variables.addSecondaryVariable(
+ "gas_molar_fraction_contaminant",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &ThermalTwoPhaseFlowWithPPLocalAssemblerInterface::
+ getIntPtGasMolFracContaminant));
}
void ThermalTwoPhaseFlowWithPPProcess::assembleConcreteProcess(
--
GitLab