Merge branch 'ShapeMatrixCacheForLF_RemoveTemplateParameter' into 'ma…

…ster' HydroThermal and ComponentTransport: Use shape matrix cache See merge request ogs/ogs!4967
ufz · Apr 8, 2024 · 1a66142 · 1a66142
2 parents fc04742 + 64cf5e7
commit 1a66142
Show file tree

Hide file tree

Showing 17 changed files with 215 additions and 100 deletions.
diff --git a/NumLib/NumericalStability/AdvectionMatrixAssembler.h b/NumLib/NumericalStability/AdvectionMatrixAssembler.h
@@ -14,12 +14,35 @@
 #include <memory>
 #include <vector>
 
+#include "NumLib/Fem/ShapeMatrixCache.h"
 #include "NumericalStabilization.h"
 
 namespace NumLib
 {
 namespace detail
 {
+template <typename MeshElementType,
+          typename IPData,
+          typename FluxVectorType,
+          typename Derived>
+void assembleAdvectionMatrix(IPData const& ip_data_vector,
+                             NumLib::ShapeMatrixCache const& shape_matrix_cache,
+                             std::vector<FluxVectorType> const& ip_flux_vector,
+                             Eigen::MatrixBase<Derived>& laplacian_matrix)
+{
+    auto const& Ns = shape_matrix_cache.NsHigherOrder<MeshElementType>();
+
+    for (std::size_t ip = 0; ip < ip_flux_vector.size(); ++ip)
+    {
+        auto const& ip_data = ip_data_vector[ip];
+        auto const w = ip_data.integration_weight;
+        auto const& dNdx = ip_data.dNdx;
+        auto const& N = Ns[ip];
+        laplacian_matrix.noalias() +=
+            N.transpose() * ip_flux_vector[ip].transpose() * dNdx * w;
+    }
+}
+
 template <typename IPData, typename FluxVectorType, typename Derived>
 void assembleAdvectionMatrix(IPData const& ip_data_vector,
                              std::vector<FluxVectorType> const& ip_flux_vector,
@@ -79,9 +102,13 @@ void applyFullUpwind(IPData const& ip_data_vector,
 }
 }  // namespace detail
 
-template <typename IPData, typename FluxVectorType, typename Derived>
+template <typename MeshElementType,
+          typename IPData,
+          typename FluxVectorType,
+          typename Derived>
 void assembleAdvectionMatrix(NumericalStabilization const& stabilizer,
                              IPData const& ip_data_vector,
+                             NumLib::ShapeMatrixCache const& shape_matrix_cache,
                              std::vector<FluxVectorType> const& ip_flux_vector,
                              double const average_velocity,
                              Eigen::MatrixBase<Derived>& laplacian_matrix)
@@ -100,8 +127,36 @@ void assembleAdvectionMatrix(NumericalStabilization const& stabilizer,
                 }
             }
 
-            detail::assembleAdvectionMatrix(ip_data_vector, ip_flux_vector,
+            detail::assembleAdvectionMatrix<MeshElementType>(
+                ip_data_vector, shape_matrix_cache, ip_flux_vector,
+                laplacian_matrix);
+        },
+        stabilizer);
+}
+
+template <typename IPData, typename FluxVectorType, typename Derived>
+void assembleAdvectionMatrix(NumericalStabilization const& stabilizer,
+                             IPData const& ip_data_vector,
+                             std::vector<FluxVectorType> const& ip_flux_vector,
+                             double const average_velocity,
+                             Eigen::MatrixBase<Derived>& laplacian_matrix)
+{
+    std::visit(
+        [&](auto&& stabilizer)
+        {
+            using Stabilizer = std::decay_t<decltype(stabilizer)>;
+            if constexpr (std::is_same_v<Stabilizer, FullUpwind>)
+            {
+                if (average_velocity > stabilizer.getCutoffVelocity())
+                {
+                    detail::applyFullUpwind(ip_data_vector, ip_flux_vector,
                                             laplacian_matrix);
+                    return;
+                }
+            }
+
+            detail::assembleAdvectionMatrix(
+                ip_data_vector, ip_flux_vector, laplacian_matrix);
         },
         stabilizer);
 }