From dfe77596f9817068ed4acd489716b4004da58d04 Mon Sep 17 00:00:00 2001
From: renchao_lu <renchao.lu@gmail.com>
Date: Tue, 22 Jan 2019 20:06:38 +0100
Subject: [PATCH] [PL] CT; Update review issues.
---
.../ComponentTransportFEM.h | 112 ++++++++----------
.../ComponentTransportProcess.cpp | 2 +-
.../CreateComponentTransportProcess.cpp | 23 ++--
.../CreateComponentTransportProcess.h | 5 +-
ProcessLib/Utils/ProcessUtils.cpp | 2 +-
5 files changed, 69 insertions(+), 75 deletions(-)
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 29f60032cc0..ccfacbe93ea 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -30,7 +30,6 @@
namespace ProcessLib
{
-class ProcessVariable;
namespace ComponentTransport
{
template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
@@ -68,7 +67,7 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
// Pressure always first
static const int pressure_index = 0;
static const int pressure_size = ShapeFunction::NPOINTS;
- // per component
+ // Size of concentration to each component
static const int concentration_size = ShapeFunction::NPOINTS;
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
@@ -93,24 +92,25 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
public:
- LocalAssemblerData(MeshLib::Element const& element,
- std::size_t const local_matrix_size,
- bool is_axially_symmetric,
- unsigned const integration_order,
- ComponentTransportProcessData const& process_data,
- std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const&
- process_variables)
+ LocalAssemblerData(
+ MeshLib::Element const& element,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ ComponentTransportProcessData const& process_data,
+ std::vector<std::reference_wrapper<ProcessVariable>> const&
+ per_process_variables)
: _element(element),
_process_data(process_data),
_integration_method(integration_order),
- _process_variables(process_variables)
+ _per_process_variables(per_process_variables)
{
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
- // _process_variables[0].size() can represents number of nodal DOFs in
- // either monolithic or staggered scheme.
+ // _per_process_variables.size() can represent number of nodal DOFs
+ // regardless of what scheme is adopted.
assert(local_matrix_size ==
- ShapeFunction::NPOINTS * _process_variables[0].size());
+ ShapeFunction::NPOINTS * _per_process_variables.size());
(void)local_matrix_size;
unsigned const n_integration_points =
@@ -139,7 +139,7 @@ public:
{
auto const local_matrix_size = local_x.size();
- int const num_nodal_dof = _process_variables[0].size();
+ int const num_nodal_dof = _per_process_variables.size();
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
@@ -152,9 +152,6 @@ public:
auto local_b = MathLib::createZeroedVector<LocalVectorType>(
local_b_data, local_matrix_size);
- // Nodal DOFs include pressure
- auto const num_comp = num_nodal_dof - 1;
-
// Get block matrices
auto Kpp = local_K.template block<pressure_size, pressure_size>(
pressure_index, pressure_index);
@@ -165,7 +162,10 @@ public:
auto local_p = Eigen::Map<const NodalVectorType>(
&local_x[pressure_index], pressure_size);
- for (int comp_id = 0; comp_id < num_comp; ++comp_id)
+ // Nodal DOFs include pressure
+ auto const number_of_components = num_nodal_dof - 1;
+ for (int component_id = 0; component_id < number_of_components;
+ ++component_id)
{
/* Partitioned assembler matrix
* | pp | pc1 | pc2 | pc3 |
@@ -177,7 +177,7 @@ public:
* | c3p | 0 | 0 | c3c3|
*/
auto concentration_index =
- pressure_size + comp_id * concentration_size;
+ pressure_size + component_id * concentration_size;
auto KCC =
local_K.template block<concentration_size, concentration_size>(
@@ -195,36 +195,25 @@ public:
auto local_C = Eigen::Map<const NodalVectorType>(
&local_x[concentration_index], concentration_size);
- // Prevent duplicate computation in loops
- Kpp.setZero();
- Mpp.setZero();
-
- assembleBlockMatrices(comp_id, t, local_C, local_p, KCC, MCC, MCp,
- MpC, Kpp, Mpp, Bp);
+ assembleBlockMatrices(component_id, t, local_C, local_p, KCC, MCC,
+ MCp, MpC, Kpp, Mpp, Bp);
}
}
void assembleBlockMatrices(
- int const comp_id,
- double const t,
+ int const component_id, double const t,
Eigen::Ref<const NodalVectorType> const& C_nodal_values,
Eigen::Ref<const NodalVectorType> const& p_nodal_values,
- Eigen::Ref<LocalBlockMatrixType>
- KCC,
- Eigen::Ref<LocalBlockMatrixType>
- MCC,
- Eigen::Ref<LocalBlockMatrixType>
- MCp,
- Eigen::Ref<LocalBlockMatrixType>
- MpC,
- Eigen::Ref<LocalBlockMatrixType>
- Kpp,
- Eigen::Ref<LocalBlockMatrixType>
- Mpp,
- Eigen::Ref<LocalSegmentVectorType>
- Bp)
+ Eigen::Ref<LocalBlockMatrixType> KCC,
+ Eigen::Ref<LocalBlockMatrixType> MCC,
+ Eigen::Ref<LocalBlockMatrixType> MCp,
+ Eigen::Ref<LocalBlockMatrixType> MpC,
+ Eigen::Ref<LocalBlockMatrixType> Kpp,
+ Eigen::Ref<LocalBlockMatrixType> Mpp,
+ Eigen::Ref<LocalSegmentVectorType> Bp)
{
- assert(comp_id <= static_cast<int>(_process_variables[0].size() - 1));
+ assert(component_id <=
+ static_cast<int>(_per_process_variables.size() - 1));
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -239,7 +228,7 @@ public:
GlobalDimMatrixType const& I(
GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
- // get material properties
+ // Get material properties
auto const& medium =
*_process_data.media_map->getMedium(_element.getID());
// Select the only valid for component transport liquid phase.
@@ -248,8 +237,8 @@ public:
// Assume that the component name is the same as the process variable
// name. Components are shifted by one because the first one is always
// pressure.
- auto const& component =
- phase.component(_process_variables[0][comp_id + 1].get().getName());
+ auto const& component = phase.component(
+ _per_process_variables[component_id + 1].get().getName());
for (std::size_t ip(0); ip < n_integration_points; ++ip)
{
@@ -283,8 +272,8 @@ public:
MaterialPropertyLib::longitudinal_dispersivity);
// Use the fluid density model to compute the density
- // TODO: concentration of which component as the argument for
- // calculation of fluid density
+ // TODO (renchao): concentration of which component as the argument
+ // for calculation of fluid density
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
vars[static_cast<int>(
@@ -353,11 +342,19 @@ public:
density * N) *
w;
- Mpp.noalias() += w * N.transpose() * porosity * drho_dp * N;
MpC.noalias() += w * N.transpose() * porosity * drho_dC * N;
- Kpp.noalias() += w * dNdx.transpose() * density * K_over_mu * dNdx;
- if (_process_data.has_gravity)
- Bp += w * density * density * dNdx.transpose() * K_over_mu * b;
+
+ // Calculate Mpp, Kpp, and bp in the first loop over components
+ if (component_id == 0)
+ {
+ Mpp.noalias() += w * N.transpose() * porosity * drho_dp * N;
+ Kpp.noalias() +=
+ w * dNdx.transpose() * density * K_over_mu * dNdx;
+
+ if (_process_data.has_gravity)
+ Bp += w * density * density * dNdx.transpose() * K_over_mu *
+ b;
+ }
}
}
@@ -374,10 +371,9 @@ public:
assert(!indices.empty());
auto const local_x = current_solution.get(indices);
- // concentration index of the component ranking first
+ // Assuming that fluid density always depends on the concentration of
+ // the component placed at the first.
auto const concentration_index = 1 * ShapeFunction::NPOINTS;
- // assuming that fluid density always depends on concentration of the
- // first component.
// get local_C and local_p
auto const C_nodal_values = Eigen::Map<const NodalVectorType>(
&local_x[concentration_index], concentration_size);
@@ -467,11 +463,7 @@ public:
return shape_matrices;
}();
- // concentration index of the component ranking first
auto const concentration_index = 1 * ShapeFunction::NPOINTS;
- // assuming that fluid density always depends on concentration of the
- // first component.
- // get local_C and local_p
auto const C_nodal_values = Eigen::Map<const NodalVectorType>(
&local_x[concentration_index], concentration_size);
auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
@@ -524,8 +516,8 @@ private:
ComponentTransportProcessData const& _process_data;
IntegrationMethod const _integration_method;
- std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const&
- _process_variables;
+ std::vector<std::reference_wrapper<ProcessVariable>> const&
+ _per_process_variables;
std::vector<
IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 8291831e17f..b4a7f72badb 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -51,7 +51,7 @@ void ComponentTransportProcess::initializeConcreteProcess(
mesh.getDimension(), mesh.getElements(), dof_table,
pv.getShapeFunctionOrder(), _local_assemblers,
mesh.isAxiallySymmetric(), integration_order, _process_data,
- _process_variables);
+ _process_variables[process_id]);
_secondary_variables.addSecondaryVariable(
"darcy_velocity",
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index e2684608558..db2ecd3bca8 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -55,41 +55,42 @@ std::unique_ptr<Process> createComponentTransportProcess(
std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
process_variables;
- // get all process variables stored in a vector before allocation
+ // Collect all process variables in a vector before allocation
// pressure first, concentration then
- auto const unalloc_process_variables = findProcessVariables(
+ auto const collected_process_variables = findProcessVariables(
variables, pv_config,
{//! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__pressure}
"pressure",
//! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__concentration}
"concentration"});
- // check number of components for each process variable
+ // Check number of components for each process variable
auto it = std::find_if(
- unalloc_process_variables.cbegin(), unalloc_process_variables.cend(),
+ collected_process_variables.cbegin(),
+ collected_process_variables.cend(),
[](std::reference_wrapper<ProcessLib::ProcessVariable> const& pv) {
return pv.get().getNumberOfComponents() != 1;
});
- if (it != unalloc_process_variables.end())
+ if (it != collected_process_variables.end())
OGS_FATAL(
- "Number of components for process variable \"%s\" should be 1 "
- "rather than %d.",
+ "Number of components for process variable '%s' should be 1 rather "
+ "than %d.",
it->get().getName().c_str(),
it->get().getNumberOfComponents());
- // allocated into a two-dimensional vector, depending on what scheme is
- // adopted
+ // Allocate the collected process variables into a two-dimensional vector,
+ // depending on what scheme is adopted
if (use_monolithic_scheme) // monolithic scheme.
{
- process_variables.push_back(std::move(unalloc_process_variables));
+ process_variables.push_back(std::move(collected_process_variables));
}
else // staggered scheme.
{
std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
per_process_variable;
- for (auto& pv : unalloc_process_variables)
+ for (auto& pv : collected_process_variables)
{
per_process_variable.emplace_back(pv);
process_variables.push_back(std::move(per_process_variable));
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
index 2c7d306756e..beb66dd3180 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.h
@@ -21,7 +21,8 @@ namespace ProcessLib
{
namespace ComponentTransport
{
-std::unique_ptr<Process> createComponentTransportProcess(MeshLib::Mesh& mesh,
+std::unique_ptr<Process> createComponentTransportProcess(
+ MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
@@ -29,7 +30,7 @@ std::unique_ptr<Process> createComponentTransportProcess(MeshLib::Mesh& mesh,
BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
std::string const& output_directory,
- std::map<int, std::unique_ptr<MaterialPropertyLib::Medium> > const& media);
+ std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
} // namespace ComponentTransport
} // namespace ProcessLib
diff --git a/ProcessLib/Utils/ProcessUtils.cpp b/ProcessLib/Utils/ProcessUtils.cpp
index 72ed3125d92..eee69ffeef4 100644
--- a/ProcessLib/Utils/ProcessUtils.cpp
+++ b/ProcessLib/Utils/ProcessUtils.cpp
@@ -71,7 +71,7 @@ std::vector<std::reference_wrapper<ProcessVariable>> findProcessVariables(
auto var_names = pv_config.getConfigParameterList<std::string>(tag);
if (var_names.empty())
- OGS_FATAL("Config tag <%s> is not found.", tag.c_str());
+ OGS_FATAL("No entity is found with config tag <%s>.", tag.c_str());
std::vector<std::string> cached_var_names;
--
GitLab