Merge pull request #28732 from joshuahansel/function-time-integral

Added `Function::timeIntegral(t1, t2, p)`

framework/doc/content/syntax/Functions/index.md

@@ -15,13 +15,20 @@ Note that there are exceptions to the rule that `Function`s only depend on
 space and time; for example, [MooseParsedFunction.md] may depend on post-processor
 values (which may depend on the solution) and scalar variable values.
 
-Moose `Function`s should override the following member functions
+Moose `Function`s should override the following member functions:
 
-- `Real value(Real, Point)` - returning the value of the function at a point in space and time
+- `Real value(Real, Point)` - the value of the function at a point in space and time
 - `Real value(ADReal, ADPoint)` - the AD enabled version of the above function
-- `Real timeDerivative(Real, Point)` - retuning the derivative of the function with respect to the first argument (time)
+- `Real timeDerivative(Real, Point)` - the derivative of the function with respect to the first argument (time)
 - `RealVectorValue gradient(Real, Point)` -  the spatial derivative with respect to the second argument
 
+and may optionally override the following member functions, which is only needed
+for some particular functionality:
+
+- `Real timeIntegral(Real t1, Real t1, const Point & p)`, which computes the
+  time integral of the function at the spatial point `p` between the time values
+  `t1` and `t2`.
+
 For vector valued functions
 
 - `RealVectorValue vectorValue(Real, Point)` - returning a vector value at a point in space and time

framework/include/functions/ConstantFunction.h

@@ -28,6 +28,8 @@ class ConstantFunction : public Function
   virtual Real timeDerivative(Real t, const Point & p) const override;
   virtual RealVectorValue gradient(Real t, const Point & p) const override;
 
+  virtual Real timeIntegral(Real t1, Real t2, const Point & p) const override;
+
 protected:
   const Real & _value;
 };

framework/include/functions/Function.h

@@ -136,12 +136,21 @@ class Function : public MooseObject,
   auto timeDerivative(const U & t, const U & x, const U & y = 0, const U & z = 0) const;
   ///@}
 
-  // Not defined
+  /// Returns the integral of the function over its domain
   virtual Real integral() const;
 
-  // Not defined
+  /// Returns the average of the function over its domain
   virtual Real average() const;
 
+  /**
+   * Computes the time integral at a spatial point between two time values
+   *
+   * @param[in] t1  Beginning time value
+   * @param[in] t2  End time value
+   * @param[in] p   Spatial point
+   */
+  virtual Real timeIntegral(Real t1, Real t2, const Point & p) const;
+
   void timestepSetup() override;
   // We will only allow initialSetup() and timestepSetup() to be overriden
   void residualSetup() override final;

framework/include/functions/PiecewiseLinearBase.h

@@ -32,6 +32,7 @@ class PiecewiseLinearBase : public PiecewiseTabularBase
   virtual RealGradient gradient(Real, const Point & p) const override;
   virtual Real integral() const override;
   virtual Real average() const override;
+  virtual Real timeIntegral(Real t1, Real t2, const Point & p) const override;
   virtual void setData(const std::vector<Real> & x, const std::vector<Real> & y) override;
 
 protected:

framework/include/utils/LinearInterpolation.h

@@ -66,10 +66,18 @@ class LinearInterpolation
   unsigned int getSampleSize() const;
 
   /**
-   * This function returns the integral of the function
+   * This function returns the integral of the function over the whole domain
    */
   Real integrate();
 
+  /**
+   * Returns the integral of the function over a specified domain
+   *
+   * @param[in] x1  First point in integral domain
+   * @param[in] x2  Second point in integral domain
+   */
+  Real integratePartial(Real x1, Real x2) const;
+
   Real domain(int i) const;
   Real range(int i) const;
 

framework/src/functions/ConstantFunction.C

@@ -50,3 +50,9 @@ ConstantFunction::gradient(Real /*t*/, const Point & /*p*/) const
 {
   return RealVectorValue(0);
 }
+
+Real
+ConstantFunction::timeIntegral(Real t1, Real t2, const Point & /*p*/) const
+{
+  return _value * (t2 - t1);
+}

framework/src/functions/Function.C

@@ -85,6 +85,12 @@ Function::timeDerivative(Real /*t*/, const Point & /*p*/) const
   return 0;
 }
 
+Real
+Function::timeIntegral(Real /*t1*/, Real /*t2*/, const Point & /*p*/) const
+{
+  mooseError("timeIntegral() not implemented.");
+}
+
 RealVectorValue
 Function::vectorValue(Real /*t*/, const Point & /*p*/) const
 {

framework/src/functions/PiecewiseLinearBase.C

@@ -85,6 +85,18 @@ PiecewiseLinearBase::average() const
          (_linear_interp->domain(_linear_interp->getSampleSize() - 1) - _linear_interp->domain(0));
 }
 
+Real
+PiecewiseLinearBase::timeIntegral(Real t1, Real t2, const Point & p) const
+{
+  if (_has_axis)
+    // function is constant in time; evaluate at arbitrary point in time
+    // and multiply by time integral domain width
+    return value(t1, p) * (t2 - t1);
+  else
+    // function is piecewise linear in time
+    return _scale_factor * _linear_interp->integratePartial(t1, t2);
+}
+
 void
 PiecewiseLinearBase::setData(const std::vector<Real> & x, const std::vector<Real> & y)
 {

framework/src/utils/LinearInterpolation.C

@@ -8,6 +8,7 @@
 //* https://www.gnu.org/licenses/lgpl-2.1.html
 
 #include "LinearInterpolation.h"
+#include "MooseUtils.h"
 
 #include "ChainedReal.h"
 
@@ -122,6 +123,130 @@ LinearInterpolation::integrate()
   return answer;
 }
 
+Real
+LinearInterpolation::integratePartial(Real xA, Real xB) const
+{
+  // integral computation below will assume that x2 > x1; if this is not the
+  // case, compute as if it is and then use identity to convert
+  Real x1, x2;
+  bool switch_bounds;
+  if (MooseUtils::absoluteFuzzyEqual(xA, xB))
+    return 0.0;
+  else if (xB > xA)
+  {
+    x1 = xA;
+    x2 = xB;
+    switch_bounds = false;
+  }
+  else
+  {
+    x1 = xB;
+    x2 = xA;
+    switch_bounds = true;
+  }
+
+  // compute integral with knowledge that x2 > x1
+  Real integral = 0.0;
+  // find minimum i : x[i] > x; if x > x[n-1], i = n
+  auto n = _x.size();
+  const unsigned int i1 =
+      x1 <= _x[n - 1] ? std::distance(_x.begin(), std::upper_bound(_x.begin(), _x.end(), x1)) : n;
+  const unsigned int i2 =
+      x2 <= _x[n - 1] ? std::distance(_x.begin(), std::upper_bound(_x.begin(), _x.end(), x2)) : n;
+  unsigned int i = i1;
+  while (i <= i2)
+  {
+    if (i == 0)
+    {
+      // note i1 = i
+      Real integral1, integral2;
+      if (_extrap)
+      {
+        const Real dydx = (_y[1] - _y[0]) / (_x[1] - _x[0]);
+        const Real y1 = _y[0] + dydx * (x1 - _x[0]);
+        integral1 = 0.5 * (y1 + _y[0]) * (_x[0] - x1);
+        if (i2 == i)
+        {
+          const Real y2 = _y[0] + dydx * (x2 - _x[0]);
+          integral2 = 0.5 * (y2 + _y[0]) * (_x[0] - x2);
+        }
+        else
+          integral2 = 0.0;
+      }
+      else
+      {
+        integral1 = _y[0] * (_x[0] - x1);
+        if (i2 == i)
+          integral2 = _y[0] * (_x[0] - x2);
+        else
+          integral2 = 0.0;
+      }
+
+      integral += integral1 - integral2;
+    }
+    else if (i == n)
+    {
+      // note i2 = i
+      Real integral1, integral2;
+      if (_extrap)
+      {
+        const Real dydx = (_y[n - 1] - _y[n - 2]) / (_x[n - 1] - _x[n - 2]);
+        const Real y2 = _y[n - 1] + dydx * (x2 - _x[n - 1]);
+        integral2 = 0.5 * (y2 + _y[n - 1]) * (x2 - _x[n - 1]);
+        if (i1 == n)
+        {
+          const Real y1 = _y[n - 1] + dydx * (x1 - _x[n - 1]);
+          integral1 = 0.5 * (y1 + _y[n - 1]) * (x1 - _x[n - 1]);
+        }
+        else
+          integral1 = 0.0;
+      }
+      else
+      {
+        integral2 = _y[n - 1] * (x2 - _x[n - 1]);
+        if (i1 == n)
+          integral1 = _y[n - 1] * (x1 - _x[n - 1]);
+        else
+          integral1 = 0.0;
+      }
+
+      integral += integral2 - integral1;
+    }
+    else
+    {
+      Real integral1;
+      if (i == i1)
+      {
+        const Real dydx = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
+        const Real y1 = _y[i - 1] + dydx * (x1 - _x[i - 1]);
+        integral1 = 0.5 * (y1 + _y[i - 1]) * (x1 - _x[i - 1]);
+      }
+      else
+        integral1 = 0.0;
+
+      Real integral2;
+      if (i == i2)
+      {
+        const Real dydx = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
+        const Real y2 = _y[i - 1] + dydx * (x2 - _x[i - 1]);
+        integral2 = 0.5 * (y2 + _y[i - 1]) * (x2 - _x[i - 1]);
+      }
+      else
+        integral2 = 0.5 * (_y[i] + _y[i - 1]) * (_x[i] - _x[i - 1]);
+
+      integral += integral2 - integral1;
+    }
+
+    i++;
+  }
+
+  // apply identity if bounds were switched
+  if (switch_bounds)
+    return -1.0 * integral;
+  else
+    return integral;
+}
+
 Real
 LinearInterpolation::domain(int i) const
 {

modules/solid_properties/include/solidproperties/ThermalFunctionSolidProperties.h

@@ -37,6 +37,8 @@ class ThermalFunctionSolidProperties : public ThermalSolidProperties
 
   virtual void cp_from_T(const Real & T, Real & cp, Real & dcp_dT) const override;
 
+  virtual Real cp_integral(const Real & T) const override;
+
   virtual Real rho_from_T(const Real & T) const override;
 
   virtual void rho_from_T(const Real & T, Real & rho, Real & drho_dT) const override;

modules/solid_properties/src/solidproperties/ThermalFunctionSolidProperties.C

@@ -44,6 +44,12 @@ ThermalFunctionSolidProperties::cp_from_T(const Real & T, Real & cp, Real & dcp_
   dcp_dT = _cp_function.timeDerivative(T);
 }
 
+Real
+ThermalFunctionSolidProperties::cp_integral(const Real & T) const
+{
+  return _cp_function.timeIntegral(_T_zero_e, T, Point(0, 0, 0));
+}
+
 Real
 ThermalFunctionSolidProperties::k_from_T(const Real & T) const
 {

modules/solid_properties/unit/include/ThermalFunctionSolidPropertiesTest.h

@@ -30,9 +30,10 @@ class ThermalFunctionSolidPropertiesTest : public MooseObjectUnitTest
     Function & rho = _fe_problem->getFunction("rho");
     rho.initialSetup();
 
-    InputParameters cp_pars = _factory.getValidParams("ParsedFunction");
-    cp_pars.set<std::string>("value") = "1200.0";
-    _fe_problem->addFunction("ParsedFunction", "cp", cp_pars);
+    InputParameters cp_pars = _factory.getValidParams("PiecewiseLinear");
+    cp_pars.set<std::vector<Real>>("x") = {700.0, 900.0};
+    cp_pars.set<std::vector<Real>>("y") = {1000.0, 1400.0};
+    _fe_problem->addFunction("PiecewiseLinear", "cp", cp_pars);
     Function & cp = _fe_problem->getFunction("cp");
     cp.initialSetup();
 

modules/solid_properties/unit/src/ThermalFunctionSolidPropertiesTest.C

@@ -28,12 +28,23 @@ TEST_F(ThermalFunctionSolidPropertiesTest, k)
  */
 TEST_F(ThermalFunctionSolidPropertiesTest, cp)
 {
-  Real T = 10.0;
+  Real T = 800.0;
 
   REL_TEST(_sp->cp_from_T(T), 1200.0, REL_TOL_SAVED_VALUE);
   DERIV_TEST(_sp->cp_from_T, T, REL_TOL_DERIVATIVE);
 }
 
+/**
+ * Test that the specific internal energy and its derivatives are
+ * correctly computed.
+ */
+TEST_F(ThermalFunctionSolidPropertiesTest, e)
+{
+  const Real T = 800.0;
+  REL_TEST(_sp->e_from_T(T), 536850.0, REL_TOL_SAVED_VALUE);
+  SPECIFIC_INTERNAL_ENERGY_TESTS(_sp, T, 1e-6, 1e-6);
+}
+
 /**
  * Test that the density and its derivatives are
  * correctly computed.

unit/include/ConstantFunctionTest.h

@@ -0,0 +1,18 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "MooseObjectUnitTest.h"
+
+class ConstantFunctionTest : public MooseObjectUnitTest
+{
+public:
+  ConstantFunctionTest() : MooseObjectUnitTest("MooseUnitApp") {}
+};

unit/include/FunctionTestUtils.h

@@ -0,0 +1,20 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "Function.h"
+
+namespace FunctionTestUtils
+{
+
+void testTimeIntegral(
+    const Function & f, Real t1, Real t2, const Point & p, unsigned int n_intervals, Real tol);
+
+}

unit/include/PiecewiseLinearTest.h

@@ -0,0 +1,22 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "MooseObjectUnitTest.h"
+
+class PiecewiseLinearTest : public MooseObjectUnitTest
+{
+public:
+  PiecewiseLinearTest() : MooseObjectUnitTest("MooseUnitApp") {}
+
+protected:
+  InputParameters getDefaultParams();
+  const Function & buildFunction(InputParameters & params, const std::string & obj_name);
+};