From ebc7c1be825e1e7033647fe84a8ee86556ecd35a Mon Sep 17 00:00:00 2001
From: renchao_lu <renchao.lu@gmail.com>
Date: Tue, 9 Jun 2020 23:38:54 +0200
Subject: [PATCH] [PL] Interpolate process solutions to int_pts.
---
.../ComponentTransportFEM.h | 15 ++++
.../ComponentTransportProcess.cpp | 71 +++++++++++++++++++
.../ComponentTransportProcess.h | 3 +
ProcessLib/LocalAssemblerInterface.h | 6 ++
ProcessLib/Process.h | 6 ++
ProcessLib/TimeLoop.cpp | 10 +++
6 files changed, 111 insertions(+)
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index af86f986468..2425bfefe41 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -916,6 +916,21 @@ public:
return flux;
}
+ std::vector<double> interpolateNodalValuesToIntegrationPoints(
+ std::vector<double> const& local_x) override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ std::vector<double> interpolated_values(n_integration_points);
+ for (unsigned ip(0); ip < n_integration_points; ++ip)
+ {
+ NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N,
+ interpolated_values[ip]);
+ }
+ return interpolated_values;
+ }
+
std::vector<double> const& getInterpolatedLocalSolution(
const double /*t*/,
std::vector<GlobalVector*> const& int_pt_x,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index e7c0af35970..d824b58ebc6 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -174,6 +174,77 @@ void ComponentTransportProcess::
_local_assemblers, pv.getActiveElementIDs(), _coupled_solutions);
}
+std::vector<GlobalVector>
+ComponentTransportProcess::interpolateNodalValuesToIntegrationPoints(
+ std::vector<GlobalVector*> const& nodal_values_vectors) const
+{
+ // Result is for each process a vector of integration point values for each
+ // element stored consecutively.
+ auto interpolateNodalValuesToIntegrationPoints =
+ [](std::size_t mesh_item_id,
+ LocalAssemblerInterface& local_assembler,
+ std::vector<
+ std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+ dof_tables,
+ std::vector<GlobalVector*> const& x,
+ std::vector<std::vector<double>>& int_pt_x) {
+ for (unsigned process_id = 0; process_id < x.size(); ++process_id)
+ {
+ auto const& dof_table = dof_tables[process_id].get();
+ auto const indices =
+ NumLib::getIndices(mesh_item_id, dof_table);
+ auto const local_x = x[process_id]->get(indices);
+
+ std::vector<double> const interpolated_values =
+ local_assembler.interpolateNodalValuesToIntegrationPoints(
+ local_x);
+
+ // For each element (mesh_item_id) concatenate the integration
+ // point values.
+ int_pt_x[process_id].insert(int_pt_x[process_id].end(),
+ interpolated_values.begin(),
+ interpolated_values.end());
+ }
+ };
+
+ // Same dof table for each primary variable.
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_tables;
+ dof_tables.reserve(nodal_values_vectors.size());
+ std::generate_n(std::back_inserter(dof_tables), nodal_values_vectors.size(),
+ [&]() { return std::ref(*_local_to_global_index_map); });
+
+ std::vector<std::vector<double>> integration_point_values_vecs(
+ nodal_values_vectors.size());
+
+ GlobalExecutor::executeDereferenced(
+ interpolateNodalValuesToIntegrationPoints, _local_assemblers,
+ dof_tables, nodal_values_vectors, integration_point_values_vecs);
+
+ // Convert from std::vector<std::vector<double>> to
+ // std::vector<GlobalVector>
+ std::vector<GlobalVector> integration_point_values_vectors;
+ integration_point_values_vectors.reserve(nodal_values_vectors.size());
+ for (auto const& integration_point_values_vec :
+ integration_point_values_vecs)
+ {
+ GlobalIndexType const size = integration_point_values_vec.size();
+ // New vector of size (elements x integration_points_per_element)
+ GlobalVector integration_point_values_vector(size);
+
+ // Copy one by one.
+ for (GlobalIndexType i = 0; i < size; ++i)
+ {
+ integration_point_values_vector.set(
+ i, integration_point_values_vec[i]);
+ }
+ integration_point_values_vectors.push_back(
+ std::move(integration_point_values_vector));
+ }
+
+ return integration_point_values_vectors;
+}
+
void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes(
const double t,
std::vector<GlobalVector*> const& integration_point_values_vectors,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 374f771f5e7..d86b67547fa 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -117,6 +117,9 @@ public:
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
int const process_id) override;
+ std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
+ std::vector<GlobalVector*> const& nodal_values_vectors) const override;
+
void extrapolateIntegrationPointValuesToNodes(
const double t,
std::vector<GlobalVector*> const& integration_point_values_vectors,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 2b09c451092..253fa45fe52 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -98,6 +98,12 @@ public:
GlobalVector const& x, double const t,
double const dt, bool const use_monolithic_scheme);
+ virtual std::vector<double> interpolateNodalValuesToIntegrationPoints(
+ std::vector<double> const& /*local_x*/)
+ {
+ return {};
+ }
+
/// Computes the flux in the point \c p_local_coords that is given in local
/// coordinates using the values from \c local_x.
/// Fits to monolithic scheme.
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 5838c56e19c..43f04a22d21 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -101,6 +101,12 @@ public:
{
}
+ virtual std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
+ std::vector<GlobalVector*> const& /*nodal_values_vectors*/) const
+ {
+ return {};
+ }
+
virtual void extrapolateIntegrationPointValuesToNodes(
const double /*t*/,
std::vector<GlobalVector*> const& /*integration_point_values_vectors*/,
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 8afb0961cb4..35d7e748bab 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -453,6 +453,11 @@ void TimeLoop::initialize()
{
BaseLib::RunTime time_phreeqc;
time_phreeqc.start();
+
+ auto& pcs = _per_process_data[0]->process;
+ auto const interpolated_process_solutions =
+ pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions);
+
_chemical_solver_interface->executeInitialCalculation(
_process_solutions);
INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
@@ -815,6 +820,11 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
// space and localized chemical equilibrium between solutes.
BaseLib::RunTime time_phreeqc;
time_phreeqc.start();
+
+ auto& pcs = _per_process_data[0]->process;
+ auto const interpolated_process_solutions =
+ pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions);
+
_chemical_solver_interface->doWaterChemistryCalculation(
_process_solutions, dt);
INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
--
GitLab