From bc3cc0f9fb47d5d2c63d0286a56d4f5418228d62 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Fri, 27 Sep 2019 17:51:06 +0200
Subject: [PATCH] [PL] H2PP; Cleanup FEM implementation a little.
Fix wrong comment.
Move media and phases out of integration loop.
Add aliases for integration weight, dNdx, mass operator.
Set position to current integration point.
Remove Sw and temperature local variables.
---
.../TwoPhaseFlowWithPPLocalAssembler-impl.h | 48 ++++++++-----------
1 file changed, 19 insertions(+), 29 deletions(-)
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
index d904dbf6b52..da2261143df 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -14,7 +14,6 @@
* pn_int_pt pressure for nonwetting phase at each integration point
* pc_int_pt capillary pressure at each integration point
* --------------secondary variable--------------------
- * temperature capillary pressure
* Sw wetting phase saturation
* dSw_dpc derivative of wetting phase saturation with respect
* to capillary pressure
@@ -92,6 +91,10 @@ void TwoPhaseFlowWithPPLocalAssembler<
auto Bl =
local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
+ auto const& medium = *_process_data.media_map->getMedium(_element.getID());
+ auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& gas_phase = medium.phase("Gas");
+
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -100,37 +103,33 @@ void TwoPhaseFlowWithPPLocalAssembler<
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
+ pos.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+ auto const& dNdx = _ip_data[ip].dNdx;
+ auto const& M = _ip_data[ip].massOperator;
+
double pc_int_pt = 0.;
double pn_int_pt = 0.;
NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, pn_int_pt,
pc_int_pt);
_pressure_wet[ip] = pn_int_pt - pc_int_pt;
-
- const double temperature = _process_data.temperature(t, pos)[0];
-
- auto const& medium =
- *_process_data.media_map->getMedium(_element.getID());
- auto const& liquid_phase = medium.phase("AqueousLiquid");
- auto const& gas_phase = medium.phase("Gas");
-
MPL::VariableArray variables;
variables[static_cast<int>(MPL::Variable::phase_pressure)] = pn_int_pt;
variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
pc_int_pt;
- variables[static_cast<int>(MPL::Variable::temperature)] = temperature;
+ variables[static_cast<int>(MPL::Variable::temperature)] =
+ _process_data.temperature(t, pos)[0];
auto const rho_nonwet = gas_phase.property(MPL::PropertyType::density)
.template value<double>(variables, pos, t);
auto const rho_wet = liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, pos, t);
-
- auto const Sw = medium.property(MPL::PropertyType::saturation)
- .template value<double>(variables, pos, t);
-
- _saturation[ip] = Sw;
+ auto& Sw = _saturation[ip];
+ Sw = medium.property(MPL::PropertyType::saturation)
+ .template value<double>(variables, pos, t);
auto const dSw_dpc =
medium.property(MPL::PropertyType::saturation)
@@ -166,20 +165,13 @@ void TwoPhaseFlowWithPPLocalAssembler<
medium.property(MPL::PropertyType::permeability)
.template value<double>(variables, pos, t));
- Mgp.noalias() +=
- porosity * (1 - Sw) * drhononwet_dpn * _ip_data[ip].massOperator;
- Mgpc.noalias() +=
- -porosity * rho_nonwet * dSw_dpc * _ip_data[ip].massOperator;
+ Mgp.noalias() += porosity * (1 - Sw) * drhononwet_dpn * M;
+ Mgpc.noalias() += -porosity * rho_nonwet * dSw_dpc * M;
+ Mlpc.noalias() += porosity * dSw_dpc * rho_wet * M;
- Mlpc.noalias() +=
- porosity * dSw_dpc * rho_wet * _ip_data[ip].massOperator;
-
- laplace_operator.noalias() = _ip_data[ip].dNdx.transpose() * K *
- _ip_data[ip].dNdx *
- _ip_data[ip].integration_weight;
+ laplace_operator.noalias() = dNdx.transpose() * K * dNdx * w;
Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator;
-
Klp.noalias() += rho_wet * lambda_wet * laplace_operator;
Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator;
@@ -187,9 +179,7 @@ void TwoPhaseFlowWithPPLocalAssembler<
{
auto const& b = _process_data.specific_body_force;
- NodalVectorType gravity_operator = _ip_data[ip].dNdx.transpose() *
- K * b *
- _ip_data[ip].integration_weight;
+ NodalVectorType gravity_operator = dNdx.transpose() * K * b * w;
Bg.noalias() +=
rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator;
Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator;
--
GitLab