Merge branch 'FCT' into 'master'

add the numerical stabilization type "FCT" and apply it to ComponentTransport process

See merge request ogs/ogs!4825

NumLib/NumericalStability/CreateNumericalStabilization.cpp

@@ -56,6 +56,10 @@ NumericalStabilization createNumericalStabilization(
 
         return FullUpwind{cutoff_velocity};
     }
+    if (type == "FluxCorrectedTransport")
+    {
+        return FluxCorrectedTransport();
+    }
 
     OGS_FATAL("The stabilization type {:s} is not available.", type);
 }

NumLib/NumericalStability/CreateNumericalStabilization.h

@@ -26,8 +26,11 @@ namespace NumLib
 struct NoStabilization;
 class IsotropicDiffusionStabilization;
 class FullUpwind;
+class FluxCorrectedTransport;
+
 using NumericalStabilization =
-    std::variant<NoStabilization, IsotropicDiffusionStabilization, FullUpwind>;
+    std::variant<NoStabilization, IsotropicDiffusionStabilization, FullUpwind,
+                 FluxCorrectedTransport>;
 }  // namespace NumLib
 
 namespace NumLib

NumLib/NumericalStability/FluxCorrectedTransport.h

@@ -0,0 +1,233 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 <Eigen/Core>
+#include <Eigen/Sparse>
+
+#include "MathLib/LinAlg/MatrixSpecifications.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "NumericalStabilization.h"
+
+namespace NumLib
+{
+namespace detail
+{
+
+template <typename MatrixVectorType>
+std::unique_ptr<MatrixVectorType> newZeroedInstance(
+    MathLib::MatrixSpecifications const& matrix_specification)
+{
+    auto result = MathLib::MatrixVectorTraits<MatrixVectorType>::newInstance(
+        matrix_specification);
+    result->setZero();
+    return result;
+}
+
+void calculateFluxCorrectedTransport(
+    [[maybe_unused]] const double t, [[maybe_unused]] const double dt,
+    [[maybe_unused]] std::vector<GlobalVector*> const& x,
+    [[maybe_unused]] std::vector<GlobalVector*> const& x_prev,
+    [[maybe_unused]] int const process_id,
+    [[maybe_unused]] const MathLib::MatrixSpecifications& matrix_specification,
+    [[maybe_unused]] GlobalMatrix& M, [[maybe_unused]] GlobalMatrix& K,
+    [[maybe_unused]] GlobalVector& b)
+{
+#ifndef USE_PETSC
+    auto D = newZeroedInstance<GlobalMatrix>(matrix_specification);
+    auto F = newZeroedInstance<GlobalMatrix>(matrix_specification);
+
+    // compute artificial diffusion operator D
+    using RawMatrixType = Eigen::SparseMatrix<double, Eigen::RowMajor>;
+    auto& K_raw = K.getRawMatrix();
+    K_raw *= -1;
+    for (int k = 0; k < K_raw.outerSize(); ++k)
+    {
+        for (RawMatrixType::InnerIterator it(K_raw, k); it; ++it)
+        {
+            if (it.row() != it.col())
+            {
+                double const kij = it.value();
+                double const kji = K.get(it.col(), it.row());
+                double const dij = std::max({-kij, 0., -kji});
+                D->setValue(it.row(), it.col(), dij);
+            }
+        }
+    }
+    auto& D_raw = D->getRawMatrix();
+    D_raw -= (D_raw * Eigen::VectorXd::Ones(D_raw.cols())).eval().asDiagonal();
+
+    // compute F
+    for (int k = 0; k < D_raw.outerSize(); ++k)
+    {
+        for (RawMatrixType::InnerIterator it(D_raw, k); it; ++it)
+        {
+            double const dij = it.value();
+            double const xj = x[process_id]->get(it.col());
+            double const xi = x[process_id]->get(it.row());
+            double const fij = -dij * (xj - xi);
+            F->setValue(it.row(), it.col(), fij);
+        }
+    }
+
+    auto& M_raw = M.getRawMatrix();
+    for (int k = 0; k < M_raw.outerSize(); ++k)
+    {
+        for (RawMatrixType::InnerIterator it(M_raw, k); it; ++it)
+        {
+            double const mij = it.value();
+            double const xdotj = (x[process_id]->get(it.col()) -
+                                  x_prev[process_id]->get(it.col())) /
+                                 dt;
+            double const xdoti = (x[process_id]->get(it.row()) -
+                                  x_prev[process_id]->get(it.row())) /
+                                 dt;
+            double const fij = -mij * (xdotj - xdoti);
+            F->add(it.row(), it.col(), fij);
+        }
+    }
+
+    auto P_plus = newZeroedInstance<GlobalVector>(matrix_specification);
+    auto P_minus = newZeroedInstance<GlobalVector>(matrix_specification);
+    auto Q_plus = newZeroedInstance<GlobalVector>(matrix_specification);
+    auto Q_minus = newZeroedInstance<GlobalVector>(matrix_specification);
+    auto R_plus = newZeroedInstance<GlobalVector>(matrix_specification);
+    auto R_minus = newZeroedInstance<GlobalVector>(matrix_specification);
+
+    auto& F_raw = F->getRawMatrix();
+    for (int k = 0; k < F_raw.outerSize(); ++k)
+    {
+        for (RawMatrixType::InnerIterator it(F_raw, k); it; ++it)
+        {
+            if (it.row() != it.col())
+            {
+                double const fij = it.value();
+                P_plus->add(it.row(), std::max(0., fij));
+                P_minus->add(it.row(), std::min(0., fij));
+
+                double const x_prev_i = x_prev[process_id]->get(it.row());
+                double const x_prev_j = x_prev[process_id]->get(it.col());
+
+                double const Q_plus_i = Q_plus->get(it.row());
+                double Q_plus_i_tmp =
+                    std::max({Q_plus_i, 0., x_prev_j - x_prev_i});
+                Q_plus->set(it.row(), Q_plus_i_tmp);
+
+                double const Q_minus_i = Q_minus->get(it.row());
+                double Q_minus_i_tmp =
+                    std::min({Q_minus_i, 0., x_prev_j - x_prev_i});
+                Q_minus->set(it.row(), Q_minus_i_tmp);
+            }
+        }
+    }
+
+    Eigen::VectorXd const M_L =
+        (M_raw * Eigen::VectorXd::Ones(M_raw.cols())).eval();
+    for (auto k = R_plus->getRangeBegin(); k < R_plus->getRangeEnd(); ++k)
+    {
+        double const mi = M_L(k);
+        double const Q_plus_i = Q_plus->get(k);
+        double const P_plus_i = P_plus->get(k);
+
+        double const tmp =
+            P_plus_i == 0. ? 0.0 : std::min(1.0, mi * Q_plus_i / dt / P_plus_i);
+
+        R_plus->set(k, tmp);
+    }
+
+    for (auto k = R_minus->getRangeBegin(); k < R_minus->getRangeEnd(); ++k)
+    {
+        double const mi = M_L(k);
+        double const Q_minus_i = Q_minus->get(k);
+        double const P_minus_i = P_minus->get(k);
+
+        double const tmp = P_minus_i == 0.
+                               ? 0.0
+                               : std::min(1.0, mi * Q_minus_i / dt / P_minus_i);
+        R_minus->set(k, tmp);
+    }
+
+    auto alpha = newZeroedInstance<GlobalMatrix>(matrix_specification);
+    for (int k = 0; k < F_raw.outerSize(); ++k)
+    {
+        for (RawMatrixType::InnerIterator it(F_raw, k); it; ++it)
+        {
+            double const fij = it.value();
+            if (fij > 0.)
+            {
+                double const R_plus_i = R_plus->get(it.row());
+                double const R_minus_j = R_minus->get(it.col());
+                double const alpha_ij = std::min(R_plus_i, R_minus_j);
+
+                alpha->setValue(it.row(), it.col(), alpha_ij);
+            }
+            else
+            {
+                double const R_minus_i = R_minus->get(it.row());
+                double const R_plus_j = R_plus->get(it.col());
+                double const alpha_ij = std::min(R_minus_i, R_plus_j);
+
+                alpha->setValue(it.row(), it.col(), alpha_ij);
+            }
+        }
+    }
+
+    // compute limited antidiffusive fluxes
+    for (int k = 0; k < F_raw.outerSize(); ++k)
+    {
+        for (RawMatrixType::InnerIterator it(F_raw, k); it; ++it)
+        {
+            if (it.row() != it.col())
+            {
+                double const fij = it.value();
+                double const alpha_ij = alpha->get(it.row(), it.col());
+
+                b.add(it.row(), alpha_ij * fij);
+            }
+        }
+    }
+
+    // compute low-order operator
+    K_raw += D_raw;
+    K_raw *= -1;
+
+    // overwrite with the lumped mass matrix
+    M.setZero();
+    for (int k = 0; k < M.getNumberOfRows(); ++k)
+    {
+        M.setValue(k, k, M_L(k));
+    }
+#endif  // end of ifndef USE_PETSC
+}
+}  // namespace detail
+
+inline void computeFluxCorrectedTransport(
+    NumericalStabilization const& stabilizer, const double t, const double dt,
+    std::vector<GlobalVector*> const& x,
+    std::vector<GlobalVector*> const& x_prev, int const process_id,
+    const MathLib::MatrixSpecifications& matrix_specification, GlobalMatrix& M,
+    GlobalMatrix& K, GlobalVector& b)
+{
+    // if needed, calling the Flux-Corrected-Transport function
+    return std::visit(
+        [&](auto&& stabilizer)
+        {
+            using Stabilizer = std::decay_t<decltype(stabilizer)>;
+            if constexpr (std::is_same_v<Stabilizer,
+                                         NumLib::FluxCorrectedTransport>)
+            {
+                return detail::calculateFluxCorrectedTransport(
+                    t, dt, x, x_prev, process_id, matrix_specification, M, K,
+                    b);
+            }
+        },
+        stabilizer);
+}
+}  // namespace NumLib

NumLib/NumericalStability/NumericalStabilization.h

@@ -15,6 +15,9 @@
 #include <variant>
 #include <vector>
 
+#include "BaseLib/Error.h"
+#include "NumLib/DOF/GlobalMatrixProviders.h"
+
 namespace NumLib
 {
 /** It defines stabilization method for solving the advection diffusion
@@ -170,6 +173,26 @@ class FullUpwind final
     double const cutoff_velocity_;
 };
 
+class FluxCorrectedTransport final
+{
+public:
+    explicit FluxCorrectedTransport()
+    {
+#ifdef USE_PETSC
+        OGS_FATAL(
+            "FluxCorrectedTransport scheme is not implemented to work with MPI "
+            "parallelization.");
+#endif
+    }
+
+    void calculateFluxCorrectedTransport(
+        const double t, const double dt, std::vector<GlobalVector*> const& x,
+        std::vector<GlobalVector*> const& x_prev, int const process_id,
+        const MathLib::MatrixSpecifications& matrix_specification,
+        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b);
+};
+
 using NumericalStabilization =
-    std::variant<NoStabilization, IsotropicDiffusionStabilization, FullUpwind>;
+    std::variant<NoStabilization, IsotropicDiffusionStabilization, FullUpwind,
+                 FluxCorrectedTransport>;
 }  // namespace NumLib

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -18,9 +18,12 @@
 #include "MathLib/LinAlg/Eigen/EigenTools.h"
 #include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
 #include "MathLib/LinAlg/FinalizeVectorAssembly.h"
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "NumLib/NumericalStability/FluxCorrectedTransport.h"
+#include "NumLib/NumericalStability/NumericalStabilization.h"
 #include "NumLib/ODESolver/NonlinearSystem.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 #include "ProcessLib/SurfaceFlux/SurfaceFlux.h"
@@ -220,6 +223,16 @@ void ComponentTransportProcess::assembleConcreteProcess(
     _asm_mat_cache.assemble(t, dt, x, x_prev, process_id, M, K, b, dof_tables,
                             _global_assembler, _local_assemblers,
                             pv.getActiveElementIDs());
+
+    if (process_id == _process_data.hydraulic_process_id)
+    {
+        return;
+    }
+    auto const matrix_specification = getMatrixSpecifications(process_id);
+
+    NumLib::computeFluxCorrectedTransport(_process_data.stabilizer, t, dt, x,
+                                          x_prev, process_id,
+                                          matrix_specification, M, K, b);
 }
 
 void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
@@ -515,6 +528,5 @@ void ComponentTransportProcess::preOutputConcreteProcess(
     INFO("[time] Computing residuum flow rates took {:g} s",
          time_residuum.elapsed());
 }
-
 }  // namespace ComponentTransport
 }  // namespace ProcessLib

ProcessLib/ComponentTransport/Tests.cmake

@@ -807,6 +807,10 @@ if (NOT OGS_USE_MPI)
     OgsTest(PROJECTFILE Parabolic/ComponentTransport/ThermalDiffusion/TemperatureField_transport.prj RUNTIME 27)
 endif()
 
+if(NOT OGS_USE_PETSC)
+    OgsTest(PROJECTFILE Parabolic/ComponentTransport/FCT_test/1d_step_func.prj RUNTIME 5)
+endif()
+
 if(NOT OGS_USE_PETSC)
     NotebookTest(
         NOTEBOOKFILE Parabolic/ComponentTransport/ReactiveTransport/DecayChain/GlobalImplicitApproach/performance_measurements.ipynb