From 2e900b92538be7094fcd3f1930f4e957dfda2ee4 Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@uio.no>
Date: Tue, 4 Jun 2024 10:44:17 +0200
Subject: [PATCH] [math] Add smoothMin/Max to common/math.hh
---
dumux/common/math.hh | 29 +++++++++++++++++++++++++++++
test/common/math/test_math.cc | 20 ++++++++++++++++++++
2 files changed, 49 insertions(+)
diff --git a/dumux/common/math.hh b/dumux/common/math.hh
index 314bc2d2f2..623dfb1dab 100644
--- a/dumux/common/math.hh
+++ b/dumux/common/math.hh
@@ -111,6 +111,35 @@ void harmonicMeanMatrix(Dune::FieldMatrix<Scalar, m, n> &K,
}
}
+/*!
+ * \ingroup Core
+ * \brief A smoothed minimum function (using cubic interpolation)
+ * \note The returned value is guaranteed to be smaller or equal to the minimum of the two input values.
+ * \note The smoothing kicks in if the difference between the two input values is smaller than the smoothing parameter.
+ * \param a The first input value
+ * \param b The second input value
+ * \param k The smoothing parameter (k > 0)
+ */
+template<class Scalar>
+Scalar smoothMin(const Scalar a, const Scalar b, const Scalar k)
+{
+ using std::max; using std::min; using std::abs;
+ const auto h = max(k-abs(a-b), 0.0 )/k;
+ return min(a, b) - h*h*h*k*(1.0/6.0);
+}
+
+/*!
+ * \ingroup Core
+ * \brief A smoothed maximum function (using cubic interpolation)
+ *
+ * \param a The first input value
+ * \param b The second input value
+ * \param k The smoothing parameter (k > 0)
+ */
+template<class Scalar>
+Scalar smoothMax(const Scalar a, const Scalar b, const Scalar k)
+{ return -smoothMin(-a, -b, k); }
+
/*!
* \ingroup Core
* \brief Invert a linear polynomial analytically
diff --git a/test/common/math/test_math.cc b/test/common/math/test_math.cc
index 0f06b0a035..e7404d5521 100644
--- a/test/common/math/test_math.cc
+++ b/test/common/math/test_math.cc
@@ -377,5 +377,25 @@ int main()
if (!Dune::FloatCmp::eq(roots[0], -3.863448718244389e-02, 1e-14))
DUNE_THROW(Dune::Exception, "Root of cubic equation does not match reference");
+ // Test smoothMin and smoothMax
+ {
+ const double a = 1.0;
+ const double b = 2.0;
+ const double k = 0.5;
+ if (!Dune::FloatCmp::eq(Dumux::smoothMin(a, b, k), a))
+ DUNE_THROW(Dune::Exception, "smoothMin does not work as expected");
+ if (!Dune::FloatCmp::eq(Dumux::smoothMax(a, b, k), b))
+ DUNE_THROW(Dune::Exception, "smoothMax does not work as expected");
+ }
+ {
+ const double a = 1.0;
+ const double b = 2.0;
+ const double k = 2.0;
+ if (Dumux::smoothMin(a, b, k) >= a)
+ DUNE_THROW(Dune::Exception, "smoothMin does not work as expected");
+ if (Dumux::smoothMax(a, b, k) <= b)
+ DUNE_THROW(Dune::Exception, "smoothMax does not work as expected");
+ }
+
return 0;
}
--
GitLab