heat transport process added

Applications/ApplicationsLib/ProjectData.cpp

@@ -33,6 +33,7 @@
 
 #include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
 #include "ProcessLib/TES/CreateTESProcess.h"
+#include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
 
 namespace detail
 {
@@ -263,6 +264,11 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config)
             process = ProcessLib::TES::createTESProcess(
                 *_mesh_vec[0], _process_variables, _parameters, process_config);
         }
+        else if (type == "HEAT_CONDUCTION")
+        {
+            process = ProcessLib::HeatConduction::createHeatConductionProcess(
+                *_mesh_vec[0], _process_variables, _parameters, process_config);
+        }
         else
         {
             OGS_FATAL("Unknown process type: %s", type.c_str());

ProcessLib/CMakeLists.txt

@@ -15,6 +15,9 @@ APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
 add_subdirectory(TES)
 APPEND_SOURCE_FILES(SOURCES TES)
 
+add_subdirectory(HeatConduction)
+APPEND_SOURCE_FILES(SOURCES HeatConduction)
+
 APPEND_SOURCE_FILES(SOURCES Utils)
 
 add_library(ProcessLib ${SOURCES})

ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "CreateHeatConductionProcess.h"
+
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "HeatConductionProcess.h"
+#include "HeatConductionProcessData.h"
+
+namespace ProcessLib
+{
+namespace HeatConduction
+{
+std::unique_ptr<Process> createHeatConductionProcess(
+    MeshLib::Mesh& mesh,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{process__type}
+    config.checkConfigParameter("type", "HEAT_CONDUCTION");
+
+    DBUG("Create HeatConductionProcess.");
+
+    // Process variable.
+    auto process_variables = findProcessVariables(
+        variables, config,
+        {//! \ogs_file_param_special{process__HEAT_CONDUCTION__process_variables__process_variable}
+         "process_variable"});
+
+    // thermal conductivity parameter.
+    auto& thermal_conductivity = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HEAT_CONDUCTION__thermal_conductivity}
+        "thermal_conductivity", parameters, 1);
+
+    DBUG("Use \'%s\' as thermal conductivity parameter.",
+         thermal_conductivity.name.c_str());
+
+    // heat capacity parameter.
+    auto& heat_capacity = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HEAT_CONDUCTION__heat_capacity}
+        "heat_capacity", parameters, 1);
+
+    DBUG("Use \'%s\' as heat capacity parameter.", heat_capacity.name.c_str());
+
+    // density parameter.
+    auto& density = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HEAT_CONDUCTION__density}
+        "density", parameters, 1);
+
+    DBUG("Use \'%s\' as density parameter.", density.name.c_str());
+
+    HeatConductionProcessData process_data{thermal_conductivity, heat_capacity,
+                                           density};
+
+    SecondaryVariableCollection secondary_variables;
+
+    NumLib::NamedFunctionCaller named_function_caller(
+        {"HeatConduction_temperature"});
+
+    ProcessLib::parseSecondaryVariables(config, secondary_variables,
+                                        named_function_caller);
+
+    return std::unique_ptr<Process>{new HeatConductionProcess{
+        mesh, parameters, std::move(process_variables), std::move(process_data),
+        std::move(secondary_variables), std::move(named_function_caller)}};
+}
+
+}  // namespace HeatConduction
+}  // namespace ProcessLib

ProcessLib/HeatConduction/CreateHeatConductionProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_CREATE_HEATCONDUCTIONPROCESS_H_
+#define PROCESS_LIB_CREATE_HEATCONDUCTIONPROCESS_H_
+
+#include <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace HeatConduction
+{
+std::unique_ptr<Process> createHeatConductionProcess(
+    MeshLib::Mesh& mesh,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace HeatConduction
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_CREATE_HEATCONDUCTIONPROCESS_H_

ProcessLib/HeatConduction/HeatConductionFEM.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_HEATCONDUCTION_FEM_H_
+#define PROCESS_LIB_HEATCONDUCTION_FEM_H_
+
+#include <vector>
+
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "HeatConductionProcessData.h"
+
+namespace ProcessLib
+{
+namespace HeatConduction
+{
+const unsigned NUM_NODAL_DOF = 1;
+
+class HeatConductionLocalAssemblerInterface
+    : public ProcessLib::LocalAssemblerInterface,
+      public NumLib::ExtrapolatableElement
+{
+public:
+    virtual std::vector<double> const& getIntPtHeatFluxX(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtHeatFluxY(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtHeatFluxZ(
+        std::vector<double>& /*cache*/) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+class LocalAssemblerData : public HeatConductionLocalAssemblerInterface
+{
+    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+    using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+    using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
+        ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
+
+    using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
+    using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
+    using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+
+public:
+    /// The thermal_conductivity factor is directly integrated into the local
+    /// element matrix.
+    LocalAssemblerData(MeshLib::Element const& element,
+                       std::size_t const local_matrix_size,
+                       unsigned const integration_order,
+                       HeatConductionProcessData const& process_data)
+        : _element(element),
+          _process_data(process_data),
+          _localK(local_matrix_size,
+                  local_matrix_size),  // TODO narrowing conversion
+          _localM(local_matrix_size, local_matrix_size),
+          _localRhs(local_matrix_size),
+          _integration_method(integration_order),
+          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                            IntegrationMethod, GlobalDim>(
+              element, _integration_method))
+    {
+        // This assertion is valid only if all nodal d.o.f. use the same shape
+        // matrices.
+        assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+    }
+
+    void assembleConcrete(
+        double const t, std::vector<double> const& local_x,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
+    {
+        _localK.setZero();
+        _localM.setZero();
+        _localRhs.setZero();
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            pos.setIntegrationPoint(ip);
+            auto const& sm = _shape_matrices[ip];
+            auto const& wp = _integration_method.getWeightedPoint(ip);
+            auto const k = _process_data.thermal_conductivity(t, pos)[0];
+            auto const heat_capacity = _process_data.heat_capacity(t, pos)[0];
+            auto const density = _process_data.density(t, pos)[0];
+
+            _localK.noalias() +=
+                sm.dNdx.transpose() * k * sm.dNdx * sm.detJ * wp.getWeight();
+            _localM.noalias() += sm.N.transpose() * density * heat_capacity *
+                                 sm.N * sm.detJ * wp.getWeight();
+            // heat flux only computed for output.
+            GlobalDimVectorType const heat_flux =
+                -k * sm.dNdx * Eigen::Map<const NodalVectorType>(
+                                    local_x.data(), ShapeFunction::NPOINTS);
+
+            for (unsigned d = 0; d < GlobalDim; ++d)
+            {
+                _heat_fluxes[d][ip] = heat_flux[d];
+            }
+        }
+
+        K.add(indices, _localK);
+        M.add(indices, _localM);
+        b.add(indices.rows, _localRhs);
+    }
+
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
+    {
+        auto const& N = _shape_matrices[integration_point].N;
+
+        // assumes N is stored contiguously in memory
+        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+    }
+
+    std::vector<double> const& getIntPtHeatFluxX(
+        std::vector<double>& /*cache*/) const override
+    {
+        assert(_heat_fluxes.size() > 0);
+        return _heat_fluxes[0];
+    }
+
+    std::vector<double> const& getIntPtHeatFluxY(
+        std::vector<double>& /*cache*/) const override
+    {
+        assert(_heat_fluxes.size() > 1);
+        return _heat_fluxes[1];
+    }
+
+    std::vector<double> const& getIntPtHeatFluxZ(
+        std::vector<double>& /*cache*/) const override
+    {
+        assert(_heat_fluxes.size() > 2);
+        return _heat_fluxes[2];
+    }
+
+private:
+    MeshLib::Element const& _element;
+    HeatConductionProcessData const& _process_data;
+
+    NodalMatrixType _localK;
+    NodalMatrixType _localM;
+    NodalVectorType _localRhs;
+
+    IntegrationMethod const _integration_method;
+    std::vector<ShapeMatrices> _shape_matrices;
+
+    std::vector<std::vector<double>> _heat_fluxes =
+        std::vector<std::vector<double>>(
+            GlobalDim, std::vector<double>(ShapeFunction::NPOINTS));
+};
+
+}  // namespace HeatConduction
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_HEATCONDUCTION_FEM_H_

ProcessLib/HeatConduction/HeatConductionProcess.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "HeatConductionProcess.h"
+
+#include <cassert>
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+namespace HeatConduction
+{
+HeatConductionProcess::HeatConductionProcess(
+    MeshLib::Mesh& mesh,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+    HeatConductionProcessData&& process_data,
+    SecondaryVariableCollection&& secondary_variables,
+    NumLib::NamedFunctionCaller&& named_function_caller)
+    : Process(mesh, parameters, std::move(process_variables),
+              std::move(secondary_variables), std::move(named_function_caller)),
+      _process_data(std::move(process_data))
+{
+}
+
+void HeatConductionProcess::initializeConcreteProcess(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh,
+    unsigned const integration_order)
+{
+    ProcessLib::createLocalAssemblers<LocalAssemblerData>(
+        mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
+        _local_assemblers, _process_data);
+
+    _secondary_variables.addSecondaryVariable(
+        "heat_flux_x", 1,
+        makeExtrapolator(
+            getExtrapolator(), _local_assemblers,
+            &HeatConductionLocalAssemblerInterface::getIntPtHeatFluxX));
+
+    if (mesh.getDimension() > 1)
+    {
+        _secondary_variables.addSecondaryVariable(
+            "heat_flux_y", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &HeatConductionLocalAssemblerInterface::getIntPtHeatFluxY));
+    }
+    if (mesh.getDimension() > 2)
+    {
+        _secondary_variables.addSecondaryVariable(
+            "heat_flux_z", 1,
+            makeExtrapolator(
+                getExtrapolator(), _local_assemblers,
+                &HeatConductionLocalAssemblerInterface::getIntPtHeatFluxZ));
+    }
+}
+
+void HeatConductionProcess::assembleConcreteProcess(const double t,
+                                                    GlobalVector const& x,
+                                                    GlobalMatrix& M,
+                                                    GlobalMatrix& K,
+                                                    GlobalVector& b)
+{
+    DBUG("Assemble HeatConductionProcess.");
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberOnDereferenced(
+        &HeatConductionLocalAssemblerInterface::assemble, _local_assemblers,
+        *_local_to_global_index_map, t, x, M, K, b);
+}
+
+}  // namespace HeatConduction
+}  // namespace ProcessLib

ProcessLib/HeatConduction/HeatConductionProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_HEATCONDUCTIONPROCESS_H_
+#define PROCESS_LIB_HEATCONDUCTIONPROCESS_H_
+
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
+#include "ProcessLib/Process.h"
+#include "HeatConductionFEM.h"
+#include "HeatConductionProcessData.h"
+
+namespace ProcessLib
+{
+namespace HeatConduction
+{
+class HeatConductionProcess final : public Process
+{
+
+public:
+    HeatConductionProcess(
+        MeshLib::Mesh& mesh,
+        std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+        std::vector<std::reference_wrapper<ProcessVariable>>&&
+            process_variables,
+        HeatConductionProcessData&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller);
+
+    //! \name ODESystem interface
+    //! @{
+
+    bool isLinear() const override { return true; }
+    //! @}
+
+private:
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override;
+
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
+
+    HeatConductionProcessData _process_data;
+
+    std::vector<std::unique_ptr<HeatConductionLocalAssemblerInterface>>
+        _local_assemblers;
+};
+
+}  // namespace HeatConduction
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_HEATCONDUCTIONPROCESS_H_

ProcessLib/HeatConduction/HeatConductionProcessData.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_HEATCONDUCTION_HEATCONDUCTIONPROCESSDATA_H
+#define PROCESSLIB_HEATCONDUCTION_HEATCONDUCTIONPROCESSDATA_H
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+template <typename T>
+struct Parameter;
+
+namespace HeatConduction
+{
+
+struct HeatConductionProcessData
+{
+    HeatConductionProcessData(Parameter<double> const& thermal_conductivity_,
+                              Parameter<double> const& heat_capacity_,
+                              Parameter<double> const& density_)
+        : thermal_conductivity(thermal_conductivity_),
+          heat_capacity(heat_capacity_),
+          density(density_)
+    {
+    }
+
+    HeatConductionProcessData(HeatConductionProcessData&& other)
+        : thermal_conductivity(other.thermal_conductivity),
+          heat_capacity(other.heat_capacity),
+          density(other.density)
+    {
+    }
+
+    //! Copies are forbidden.
+    HeatConductionProcessData(HeatConductionProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(HeatConductionProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(HeatConductionProcessData&&) = delete;
+
+    Parameter<double> const& thermal_conductivity;
+    Parameter<double> const& heat_capacity;
+    Parameter<double> const& density;
+};
+
+}  // namespace HeatConduction
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_HEATCONDUCTION_HEATCONDUCTIONPROCESSDATA_H