From 52289b1f79dbc0b2ce33df21f9fba3d1ea0a3e2e Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Tue, 12 Jul 2022 23:46:57 +0200
Subject: [PATCH] [PL/TH2M] Add/correct residuum output
The residuum consists of heat flux, gas and liquid mass flow rates,
and nodal forces.
While the residuum for the Newton scheme is readily available,
some computations are necessary in the Picard scheme.
---
ProcessLib/TH2M/TH2MProcess.cpp | 49 +++++++++++++++++----------------
ProcessLib/TH2M/TH2MProcess.h | 4 ++-
2 files changed, 29 insertions(+), 24 deletions(-)
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index 08719c89922..bf2a5004746 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -15,6 +15,7 @@
#include "MeshLib/Elements/Utils.h"
#include "NumLib/DOF/ComputeSparsityPattern.h"
#include "ProcessLib/Process.h"
+#include "ProcessLib/Utils/ComputeResiduum.h"
#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
#include "ProcessLib/Utils/SetIPDataInitialConditions.h"
#include "TH2MFEM.h"
@@ -41,12 +42,15 @@ TH2MProcess<DisplacementDim>::TH2MProcess(
std::move(secondary_variables), use_monolithic_scheme),
_process_data(std::move(process_data))
{
+ _heat_flow_rate = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "HeatFlowRate", MeshLib::MeshItemType::Node, 1);
+ _gas_mass_flow_rate = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "GasMassFlowRate", MeshLib::MeshItemType::Node, 1);
+ _liquid_mass_flow_rate = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "LiquidMassFlowRate", MeshLib::MeshItemType::Node, 1);
_nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
- _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
- mesh, "MassFlowRate", MeshLib::MeshItemType::Node, 1);
-
// TODO (naumov) remove ip suffix. Probably needs modification of the mesh
// properties, s.t. there is no "overlapping" with cell/point data.
// See getOrCreateMeshProperty.
@@ -316,6 +320,18 @@ void TH2MProcess<DisplacementDim>::assembleConcreteProcess(
_global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
b);
+
+ auto const residuum = computeResiduum(*x[0], *xdot[0], M, K, b);
+ auto copyRhs = [&](int const variable_id, auto& output_vector)
+ {
+ transformVariableFromGlobalVector(residuum, variable_id, dof_table[0],
+ output_vector, std::negate<double>());
+ };
+
+ copyRhs(0, *_gas_mass_flow_rate);
+ copyRhs(1, *_liquid_mass_flow_rate);
+ copyRhs(2, *_heat_flow_rate);
+ copyRhs(3, *_nodal_forces);
}
template <int DisplacementDim>
@@ -347,27 +363,14 @@ void TH2MProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
auto copyRhs = [&](int const variable_id, auto& output_vector)
{
- if (_use_monolithic_scheme)
- {
- transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
- output_vector,
- std::negate<double>());
- }
- else
- {
- transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
- output_vector,
- std::negate<double>());
- }
+ transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
+ output_vector, std::negate<double>());
};
- if (_use_monolithic_scheme || process_id == 1)
- {
- copyRhs(0, *_hydraulic_flow);
- }
- if (_use_monolithic_scheme || process_id == 2)
- {
- copyRhs(1, *_nodal_forces);
- }
+
+ copyRhs(0, *_gas_mass_flow_rate);
+ copyRhs(1, *_liquid_mass_flow_rate);
+ copyRhs(2, *_heat_flow_rate);
+ copyRhs(3, *_nodal_forces);
}
template <int DisplacementDim>
diff --git a/ProcessLib/TH2M/TH2MProcess.h b/ProcessLib/TH2M/TH2MProcess.h
index 388bd4855fb..e94b0bfd685 100644
--- a/ProcessLib/TH2M/TH2MProcess.h
+++ b/ProcessLib/TH2M/TH2MProcess.h
@@ -138,8 +138,10 @@ private:
return _use_monolithic_scheme || process_id == deformation_process_id;
}
+ MeshLib::PropertyVector<double>* _heat_flow_rate = nullptr;
+ MeshLib::PropertyVector<double>* _gas_mass_flow_rate = nullptr;
+ MeshLib::PropertyVector<double>* _liquid_mass_flow_rate = nullptr;
MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
- MeshLib::PropertyVector<double>* _hydraulic_flow = nullptr;
static constexpr int monolithic_process_id = 0;
static constexpr int deformation_process_id = 3;
--
GitLab