From 0c48a4698d4f0579c9770b64dae9f5ff28e6f30d Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Fri, 13 Jan 2023 12:00:40 +0100
Subject: [PATCH] [PL/CT] Add Residuum output for all variables
---
.../ComponentTransportProcess.cpp | 27 ++++++++++++++-----
.../ComponentTransportProcess.h | 2 +-
2 files changed, 22 insertions(+), 7 deletions(-)
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 85d12f1cc0f..88022dfc031 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -21,6 +21,7 @@
#include "NumLib/DOF/ComputeSparsityPattern.h"
#include "ProcessLib/SurfaceFlux/SurfaceFlux.h"
#include "ProcessLib/SurfaceFlux/SurfaceFluxData.h"
+#include "ProcessLib/Utils/ComputeResiduum.h"
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
namespace ProcessLib
@@ -48,9 +49,14 @@ ComponentTransportProcess::ComponentTransportProcess(
_surfaceflux(std::move(surfaceflux)),
_chemical_solver_interface(std::move(chemical_solver_interface))
{
- // Todo(renchao): Need to adapt for the multi-component case.
- _molar_flow_rate = MeshLib::getOrCreateMeshProperty<double>(
- mesh, "MolarFlowRate", MeshLib::MeshItemType::Node, 1);
+ std::size_t n_processes = _process_variables.size();
+ _molar_flow_rate.resize(n_processes - 1);
+ for (int process_id = 0; process_id < n_processes - 1; ++process_id)
+ {
+ _molar_flow_rate[process_id] = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "MolarFlowRate_" + std::to_string(process_id),
+ MeshLib::MeshItemType::Node, 1);
+ }
}
void ComponentTransportProcess::initializeConcreteProcess(
@@ -170,6 +176,15 @@ void ComponentTransportProcess::assembleConcreteProcess(
_global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
b);
+
+ MathLib::finalizeMatrixAssembly(M);
+ MathLib::finalizeMatrixAssembly(K);
+ MathLib::finalizeVectorAssembly(b);
+
+ auto const residuum = computeResiduum(*x[0], *xdot[0], M, K, b);
+ transformVariableFromGlobalVector(residuum, 0, *_local_to_global_index_map,
+ *_molar_flow_rate[process_id],
+ std::negate<double>());
}
void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
@@ -201,9 +216,9 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
// b is the negated residumm used in the Newton's method.
// Here negating b is to recover the primitive residuum.
- transformVariableFromGlobalVector(b, 0 /*variable id*/,
- *_local_to_global_index_map,
- *_molar_flow_rate, std::negate<double>());
+ transformVariableFromGlobalVector(
+ b, process_id, *_local_to_global_index_map,
+ *_molar_flow_rate[process_id], std::negate<double>());
}
Eigen::Vector3d ComponentTransportProcess::getFlux(
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 8c48ca410fb..da1ac9fb79b 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -170,7 +170,7 @@ private:
std::unique_ptr<ChemistryLib::ChemicalSolverInterface>
_chemical_solver_interface;
- MeshLib::PropertyVector<double>* _molar_flow_rate = nullptr;
+ std::vector<MeshLib::PropertyVector<double>*> _molar_flow_rate;
};
} // namespace ComponentTransport
--
GitLab