Merge branch 'StoreChemicalSystemOnIntPt' into 'master'

[PL/CT] Store chemical system id at the gauss integration point

See merge request ogs/ogs!3345

ChemistryLib/ChemicalSolverInterface.h

@@ -33,9 +33,15 @@ public:
         return {};
     }
 
+    virtual void computeSecondaryVariable(
+        std::size_t const /*ele_id*/,
+        std::vector<GlobalIndexType> const& /*chemical_system_indices*/)
+    {
+    }
+
     virtual ~ChemicalSolverInterface() = default;
 
 public:
-    std::vector<std::vector<GlobalIndexType>> chemical_system_index_map;
+    std::vector<GlobalIndexType> chemical_system_index_map;
 };
 }  // namespace ChemistryLib

ChemistryLib/PhreeqcIO.h

@@ -69,6 +69,10 @@ public:
 
     friend std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io);
 
+    void computeSecondaryVariable(
+        std::size_t const ele_id,
+        std::vector<GlobalIndexType> const& chemical_system_indices) override;
+
     std::vector<std::string> const getComponentList() const override;
 
     std::string const _phreeqc_input_file;
@@ -92,7 +96,7 @@ private:
     Knobs const _knobs;
     double _dt = std::numeric_limits<double>::quiet_NaN();
     const int phreeqc_instance_id = 0;
-    int _num_chemical_systems = -1;
+    std::size_t _num_chemical_systems = -1;
 };
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib

ProcessLib/ComponentTransport/ChemicalProcessData.h

-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <vector>
-
-#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
-
-namespace ProcessLib
-{
-namespace ComponentTransport
-{
-struct ChemicalProcessData
-{
-    ChemicalProcessData(
-        std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map_)
-        : chemical_system_index_map(chemical_system_index_map_)
-    {
-    }
-
-    std::vector<std::vector<GlobalIndexType>>& chemical_system_index_map;
-};
-}  // namespace ComponentTransport
-}  // namespace ProcessLib

ProcessLib/ComponentTransport/ComponentTransportFEM.h

@@ -48,6 +48,8 @@ struct IntegrationPointData final
     GlobalDimNodalMatrixType const dNdx;
     double const integration_weight;
 
+    GlobalIndexType chemical_system_id = 0;
+
     double porosity = std::numeric_limits<double>::quiet_NaN();
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
@@ -66,6 +68,8 @@ public:
         _coupled_solutions = coupling_term;
     }
 
+    virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
+
     virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double t,
         std::vector<GlobalVector*> const& x,
@@ -159,23 +163,23 @@ public:
                 medium->property(MaterialPropertyLib::PropertyType::porosity)
                     .template initialValue<double>(pos, 0.);
         }
+    }
 
-        if (_process_data.chemical_process_data)
-        {
-            // chemical system index map
-            auto& chemical_system_index_map =
-                _process_data.chemical_process_data->chemical_system_index_map;
-
-            GlobalIndexType const start_value =
-                chemical_system_index_map.empty()
-                    ? 0
-                    : chemical_system_index_map.back().back() + 1;
-
-            std::vector<GlobalIndexType> indices(
-                _integration_method.getNumberOfPoints());
-            std::iota(indices.begin(), indices.end(), start_value);
+    void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
+    {
+        assert(_process_data.chemical_solver_interface);
+        // chemical system index map
+        auto& chemical_system_index_map =
+            _process_data.chemical_solver_interface->chemical_system_index_map;
 
-            chemical_system_index_map.push_back(std::move(indices));
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data[ip].chemical_system_id = chemical_system_index_map.empty()
+                                     ? 0
+                                     : chemical_system_index_map.back() + 1;
+            chemical_system_index_map.push_back(_ip_data[ip].chemical_system_id);
         }
     }
 
@@ -945,14 +949,19 @@ public:
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
         std::vector<double>& cache) const override
     {
-        assert(_process_data.chemical_process_data);
+        assert(_process_data.chemical_solver_interface);
         assert(int_pt_x.size() == 1);
 
         cache.clear();
-        auto const ele_id = _element.getID();
-        auto const& indices = _process_data.chemical_process_data
-                                  ->chemical_system_index_map[ele_id];
-        cache = int_pt_x[0]->get(indices);
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+        for (unsigned ip(0); ip < n_integration_points; ++ip)
+        {
+            auto const& chemical_system_id = _ip_data[ip].chemical_system_id;
+            auto const c_int_pt = int_pt_x[0]->get(chemical_system_id);
+            cache.push_back(c_int_pt);
+        }
 
         return cache;
     }
@@ -976,11 +985,24 @@ public:
         auto const ele_velocity_mat =
             MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
 
-        auto ele_id = _element.getID();
+        auto const ele_id = _element.getID();
         Eigen::Map<LocalVectorType>(
             &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
             GlobalDim) =
             ele_velocity_mat.rowwise().sum() / n_integration_points;
+
+        if (_process_data.chemical_solver_interface)
+        {
+            std::vector<GlobalIndexType> chemical_system_indices;
+            chemical_system_indices.reserve(n_integration_points);
+            std::transform(
+                _ip_data.begin(), _ip_data.end(),
+                std::back_inserter(chemical_system_indices),
+                [](auto const& ip_data) { return ip_data.chemical_system_id; });
+
+            _process_data.chemical_solver_interface->computeSecondaryVariable(
+                ele_id, chemical_system_indices);
+        }
     }
 
     void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -86,6 +86,10 @@ void ComponentTransportProcess::initializeConcreteProcess(
 
     if (_chemical_solver_interface)
     {
+        GlobalExecutor::executeSelectedMemberOnDereferenced(
+            &ComponentTransportLocalAssemblerInterface::setChemicalSystemID,
+            _local_assemblers, pv.getActiveElementIDs());
+
         _chemical_solver_interface->initialize();
     }
 

ProcessLib/ComponentTransport/ComponentTransportProcessData.h

@@ -12,7 +12,7 @@
 
 #include <memory>
 
-#include "ChemicalProcessData.h"
+#include "ChemistryLib/ChemicalSolverInterface.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 
@@ -38,7 +38,7 @@ struct ComponentTransportProcessData
     Eigen::VectorXd const specific_body_force;
     bool const has_gravity;
     bool const non_advective_form;
-    std::unique_ptr<ChemicalProcessData> chemical_process_data;
+    ChemistryLib::ChemicalSolverInterface* const chemical_solver_interface;
 
     const int hydraulic_process_id;
     // TODO (renchao-lu): This variable is used in the calculation of the

ProcessLib/ComponentTransport/CreateChemicalProcessData.cpp

-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "CreateChemicalProcessData.h"
-#include "ChemicalProcessData.h"
-
-#include "ChemistryLib/ChemicalSolverInterface.h"
-
-namespace ProcessLib
-{
-namespace ComponentTransport
-{
-std::unique_ptr<ChemicalProcessData> createChemicalProcessData(
-    ChemistryLib::ChemicalSolverInterface* const chemical_solver_interface)
-{
-    if (!chemical_solver_interface)
-    {
-        return nullptr;
-    }
-
-    return std::make_unique<ChemicalProcessData>(
-        chemical_solver_interface->chemical_system_index_map);
-}
-}  // namespace ComponentTransport
-}  // namespace ProcessLib

ProcessLib/ComponentTransport/CreateChemicalProcessData.h

-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <memory>
-
-namespace ChemistryLib
-{
-class ChemicalSolverInterface;
-}
-
-namespace ProcessLib
-{
-namespace ComponentTransport
-{
-struct ChemicalProcessData;
-
-std::unique_ptr<ChemicalProcessData> createChemicalProcessData(
-    ChemistryLib::ChemicalSolverInterface* chemical_solver_interface);
-}  // namespace ComponentTransport
-}  // namespace ProcessLib

ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp

@@ -9,7 +9,6 @@
  */
 
 #include "CreateComponentTransportProcess.h"
-#include "CreateChemicalProcessData.h"
 
 #include "ChemistryLib/ChemicalSolverInterface.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
@@ -225,14 +224,11 @@ std::unique_ptr<Process> createComponentTransportProcess(
     checkMPLProperties(mesh, *media_map);
     DBUG("Media properties verified.");
 
-    auto chemical_process_data =
-        createChemicalProcessData(chemical_solver_interface.get());
-
     ComponentTransportProcessData process_data{std::move(media_map),
                                                specific_body_force,
                                                has_gravity,
                                                non_advective_form,
-                                               std::move(chemical_process_data),
+                                               chemical_solver_interface.get(),
                                                hydraulic_process_id,
                                                first_transport_process_id};