Merge pull request #2304 from renchao-lu/MultiComponentTransportProcess

[PL] Multi-component transport process

Applications/ApplicationsLib/ProjectData.cpp

@@ -517,7 +517,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
                 ProcessLib::ComponentTransport::createComponentTransportProcess(
                     *_mesh_vec[0], std::move(jacobian_assembler),
                     _process_variables, _parameters, integration_order,
-                    process_config, _mesh_vec, output_directory);
+                    process_config, _mesh_vec, output_directory, _media);
         }
         else
 #endif

MaterialLib/MPL/PropertyType.h

@@ -49,9 +49,11 @@ enum PropertyType : int
     entry_pressure,
     fredlund_parameters,
     heat_capacity,
+    /// used to compute the hydrodynamic dispersion tensor.
     longitudinal_dispersivity,
     molar_mass,
     mole_fraction,
+    /// used to compute the hydrodynamic dispersion tensor.
     molecular_diffusion,
     name,
     permeability,
@@ -67,6 +69,7 @@ enum PropertyType : int
     specific_heat_capacity,
     thermal_conductivity,
     thermal_expansivity,
+    /// used to compute the hydrodynamic dispersion tensor.
     transversal_dispersivity,
     viscosity,
     number_of_properties

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -50,7 +50,8 @@ void ComponentTransportProcess::initializeConcreteProcess(
     ProcessLib::createLocalAssemblers<LocalAssemblerData>(
         mesh.getDimension(), mesh.getElements(), dof_table,
         pv.getShapeFunctionOrder(), _local_assemblers,
-        mesh.isAxiallySymmetric(), integration_order, _process_data);
+        mesh.isAxiallySymmetric(), integration_order, _process_data,
+        _process_variables[process_id]);
 
     _secondary_variables.addSecondaryVariable(
         "darcy_velocity",
@@ -144,4 +145,3 @@ void ComponentTransportProcess::postTimestepConcreteProcess(
 
 }  // namespace ComponentTransport
 }  // namespace ProcessLib
-

ProcessLib/ComponentTransport/ComponentTransportProcessData.h

@@ -12,10 +12,16 @@
 #include <memory>
 
 #include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
 #include "MaterialLib/PorousMedium/PorousMediaProperties.h"
 
+namespace MaterialPropertyLib
+{
+class Medium;
+}
+
 namespace ProcessLib
 {
 template <typename ReturnType>
@@ -31,41 +37,24 @@ struct ComponentTransportProcessData
         ProcessLib::Parameter<double> const& fluid_reference_density_,
         std::unique_ptr<MaterialLib::Fluid::FluidProperties>&&
             fluid_properties_,
-        ProcessLib::Parameter<double> const& molecular_diffusion_coefficient_,
-        ProcessLib::Parameter<double> const& solute_dispersivity_longitudinal_,
-        ProcessLib::Parameter<double> const& solute_dispersivity_transverse_,
+        std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>&&
+            media_map_,
         ProcessLib::Parameter<double> const& retardation_factor_,
         ProcessLib::Parameter<double> const& decay_rate_,
         Eigen::VectorXd const& specific_body_force_,
         bool const has_gravity_)
         : porous_media_properties(std::move(porous_media_properties_)),
           fluid_reference_density(fluid_reference_density_),
           fluid_properties(std::move(fluid_properties_)),
-          molecular_diffusion_coefficient(molecular_diffusion_coefficient_),
-          solute_dispersivity_longitudinal(solute_dispersivity_longitudinal_),
-          solute_dispersivity_transverse(solute_dispersivity_transverse_),
+          media_map(std::move(media_map_)),
           retardation_factor(retardation_factor_),
           decay_rate(decay_rate_),
           specific_body_force(specific_body_force_),
           has_gravity(has_gravity_)
     {
     }
 
-    ComponentTransportProcessData(ComponentTransportProcessData&& other)
-        : porous_media_properties(std::move(other.porous_media_properties)),
-          fluid_reference_density(other.fluid_reference_density),
-          fluid_properties(other.fluid_properties.release()),
-          molecular_diffusion_coefficient(
-              other.molecular_diffusion_coefficient),
-          solute_dispersivity_longitudinal(
-              other.solute_dispersivity_longitudinal),
-          solute_dispersivity_transverse(other.solute_dispersivity_transverse),
-          retardation_factor(other.retardation_factor),
-          decay_rate(other.decay_rate),
-          specific_body_force(other.specific_body_force),
-          has_gravity(other.has_gravity)
-    {
-    }
+    ComponentTransportProcessData(ComponentTransportProcessData&&) = default;
 
     //! Copies are forbidden.
     ComponentTransportProcessData(ComponentTransportProcessData const&) =
@@ -80,9 +69,8 @@ struct ComponentTransportProcessData
     MaterialLib::PorousMedium::PorousMediaProperties porous_media_properties;
     Parameter<double> const& fluid_reference_density;
     std::unique_ptr<MaterialLib::Fluid::FluidProperties> fluid_properties;
-    Parameter<double> const& molecular_diffusion_coefficient;
-    Parameter<double> const& solute_dispersivity_longitudinal;
-    Parameter<double> const& solute_dispersivity_transverse;
+    std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
+        media_map;
     Parameter<double> const& retardation_factor;
     Parameter<double> const& decay_rate;
     Eigen::VectorXd const specific_body_force;

ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp

@@ -10,6 +10,7 @@
 #include "CreateComponentTransportProcess.h"
 
 #include "MaterialLib/Fluid/FluidProperties/CreateFluidProperties.h"
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MaterialLib/PorousMedium/CreatePorousMediaProperties.h"
 
 #include "MeshLib/IO/readMeshFromFile.h"
@@ -21,7 +22,6 @@
 
 #include "ComponentTransportProcess.h"
 #include "ComponentTransportProcessData.h"
-
 namespace ProcessLib
 {
 namespace ComponentTransport
@@ -34,7 +34,8 @@ std::unique_ptr<Process> createComponentTransportProcess(
     unsigned const integration_order,
     BaseLib::ConfigTree const& config,
     std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
-    std::string const& output_directory)
+    std::string const& output_directory,
+    std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media)
 {
     //! \ogs_file_param{prj__processes__process__type}
     config.checkConfigParameter("type", "ComponentTransport");
@@ -53,28 +54,46 @@ std::unique_ptr<Process> createComponentTransportProcess(
 
     std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
         process_variables;
+
+    // Collect all process variables in a vector before allocation
+    // pressure first, concentration then
+    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
+    auto it = std::find_if(
+        collected_process_variables.cbegin(),
+        collected_process_variables.cend(),
+        [](std::reference_wrapper<ProcessLib::ProcessVariable> const& pv) {
+            return pv.get().getNumberOfComponents() != 1;
+        });
+
+    if (it != collected_process_variables.end())
+        OGS_FATAL(
+            "Number of components for process variable '%s' should be 1 rather "
+            "than %d.",
+            it->get().getName().c_str(),
+            it->get().getNumberOfComponents());
+
+    // Allocate the collected process variables into a two-dimensional vector,
+    // depending on what scheme is adopted
     if (use_monolithic_scheme)  // monolithic scheme.
     {
-        auto per_process_variables = findProcessVariables(
-            variables, pv_config,
-            {
-            //! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__concentration}
-             "concentration",
-             //! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__pressure}
-             "pressure"});
-        process_variables.push_back(std::move(per_process_variables));
+        process_variables.push_back(std::move(collected_process_variables));
     }
     else  // staggered scheme.
     {
-        std::array<std::string, 2> variable_names = {
-            {"concentration",
-             "pressure"}};  // double-braces required in C++11 (not in C++14)
+        std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
+            per_process_variable;
 
-        for (int i = 0; i < 2; i++)
+        for (auto& pv : collected_process_variables)
         {
-            auto per_process_variables =
-                findProcessVariables(variables, pv_config, {variable_names[i]});
-            process_variables.push_back(std::move(per_process_variables));
+            per_process_variable.emplace_back(pv);
+            process_variables.push_back(std::move(per_process_variable));
         }
     }
 
@@ -96,30 +115,6 @@ std::unique_ptr<Process> createComponentTransportProcess(
     DBUG("Use '%s' as fluid_reference_density parameter.",
          fluid_reference_density.name.c_str());
 
-    // Parameter for the longitudinal molecular diffusion coefficient.
-    auto const& molecular_diffusion_coefficient = findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__ComponentTransport__molecular_diffusion_coefficient}
-        "molecular_diffusion_coefficient", parameters, 1);
-    DBUG("Use '%s' as molecular diffusion coefficient parameter.",
-         molecular_diffusion_coefficient.name.c_str());
-
-    // Parameter for the longitudinal solute dispersivity.
-    auto const& solute_dispersivity_longitudinal = findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__ComponentTransport__solute_dispersivity_longitudinal}
-        "solute_dispersivity_longitudinal", parameters, 1);
-    DBUG("Use '%s' as longitudinal solute dispersivity parameter.",
-         solute_dispersivity_longitudinal.name.c_str());
-
-    // Parameter for the transverse solute dispersivity.
-    auto const& solute_dispersivity_transverse = findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__ComponentTransport__solute_dispersivity_transverse}
-        "solute_dispersivity_transverse", parameters, 1);
-    DBUG("Use '%s' as transverse solute dispersivity parameter.",
-         solute_dispersivity_transverse.name.c_str());
-
     // Parameter for the retardation factor.
     auto const& retardation_factor =
         findParameter<double>(config,
@@ -150,13 +145,14 @@ std::unique_ptr<Process> createComponentTransportProcess(
         std::copy_n(b.data(), b.size(), specific_body_force.data());
     }
 
+    auto media_map =
+        MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
+
     ComponentTransportProcessData process_data{
         std::move(porous_media_properties),
         fluid_reference_density,
         std::move(fluid_properties),
-        molecular_diffusion_coefficient,
-        solute_dispersivity_longitudinal,
-        solute_dispersivity_transverse,
+        std::move(media_map),
         retardation_factor,
         decay_rate,
         specific_body_force,

ProcessLib/ComponentTransport/CreateComponentTransportProcess.h

@@ -12,6 +12,11 @@
 #include <memory>
 #include "ProcessLib/Process.h"
 
+namespace MaterialPropertyLib
+{
+class Medium;
+}
+
 namespace ProcessLib
 {
 namespace ComponentTransport
@@ -24,7 +29,8 @@ std::unique_ptr<Process> createComponentTransportProcess(
     unsigned const integration_order,
     BaseLib::ConfigTree const& config,
     std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
-    std::string const& output_directory);
+    std::string const& output_directory,
+    std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
 
 }  // namespace ComponentTransport
 }  // namespace ProcessLib