Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : #include "KokkosPiecewiseTabularBase.h"
      13             : #include "KokkosUtils.h"
      14             : 
      15             : /**
      16             :  * Function which provides a piecewise constant interpolation of a provided (x,y) point data set.
      17             :  */
      18             : class KokkosPiecewiseConstant : public KokkosPiecewiseTabularBase
      19             : {
      20             : public:
      21             :   static InputParameters validParams();
      22             : 
      23             :   KokkosPiecewiseConstant(const InputParameters & parameters);
      24             : 
      25             :   using Real3 = Moose::Kokkos::Real3;
      26             : 
      27             :   KOKKOS_FUNCTION Real value(Real t, Real3 p) const;
      28             :   KOKKOS_FUNCTION Real integral() const;
      29             :   KOKKOS_FUNCTION Real average() const;
      30             : 
      31             : private:
      32             :   /// Enum for which direction to apply values
      33             :   const enum class Direction { LEFT, RIGHT, LEFT_INCLUSIVE, RIGHT_INCLUSIVE } _direction;
      34             : };
      35             : 
      36             : KOKKOS_FUNCTION inline Real
      37     8611736 : KokkosPiecewiseConstant::value(Real t, Real3 p) const
      38             : {
      39             :   using Moose::Kokkos::Utils::sign;
      40             : 
      41     8611736 :   const Real x = _has_axis ? p(_axis) : t;
      42             : 
      43     8611736 :   const auto len = functionSize();
      44     8611736 :   constexpr Real tolerance = 1.0e-14;
      45             : 
      46             :   // endpoint cases
      47     8611736 :   if ((_direction == Direction::LEFT && x < (1 + tolerance * sign(domain(0))) * domain(0)) ||
      48     8611736 :       (_direction == Direction::RIGHT && x < (1 - tolerance * sign(domain(0))) * domain(0)) ||
      49     8611736 :       (_direction == Direction::LEFT_INCLUSIVE &&
      50    17223472 :        x < (1 - tolerance * sign(domain(0))) * domain(0)) ||
      51     8611736 :       (_direction == Direction::RIGHT_INCLUSIVE &&
      52           0 :        x < (1 + tolerance * sign(domain(0))) * domain(0)))
      53           0 :     return _scale_factor * range(0);
      54    25835208 :   else if ((_direction == Direction::LEFT &&
      55     8611736 :             x > (1 + tolerance * sign(domain(len - 1))) * domain(len - 1)) ||
      56     4376456 :            (_direction == Direction::RIGHT &&
      57           0 :             x > (1 - tolerance * sign(domain(len - 1))) * domain(len - 1)) ||
      58     4376456 :            (_direction == Direction::LEFT_INCLUSIVE &&
      59    17223472 :             x > (1 - tolerance * sign(domain(len - 1))) * domain(len - 1)) ||
      60     4376456 :            (_direction == Direction::RIGHT_INCLUSIVE &&
      61           0 :             x > (1 + tolerance * sign(domain(len - 1))) * domain(len - 1)))
      62     4235280 :     return _scale_factor * range(len - 1);
      63             : 
      64     4376456 :   for (unsigned int i = 1; i < len; ++i)
      65             :   {
      66     4376456 :     if (_direction == Direction::LEFT && x < (1 + tolerance * sign(domain(i))) * domain(i))
      67     4376456 :       return _scale_factor * range(i - 1);
      68           0 :     else if (_direction == Direction::LEFT_INCLUSIVE &&
      69           0 :              x < (1 - tolerance * sign(domain(i))) * domain(i))
      70           0 :       return _scale_factor * range(i - 1);
      71           0 :     else if ((_direction == Direction::RIGHT && x < (1 - tolerance * sign(domain(i))) * domain(i)))
      72           0 :       return _scale_factor * range(i);
      73           0 :     else if ((_direction == Direction::RIGHT_INCLUSIVE &&
      74           0 :               x < (1 + tolerance * sign(domain(i))) * domain(i)))
      75           0 :       return _scale_factor * range(i);
      76             :   }
      77             : 
      78           0 :   return 0.0;
      79             : }
      80             : 
      81             : KOKKOS_FUNCTION inline Real
      82           0 : KokkosPiecewiseConstant::integral() const
      83             : {
      84           0 :   const auto len = functionSize();
      85             : 
      86           0 :   unsigned int offset = 0;
      87             : 
      88           0 :   if (_direction == Direction::RIGHT || _direction == Direction::RIGHT_INCLUSIVE)
      89           0 :     offset = 1;
      90             : 
      91           0 :   Real sum = 0;
      92             : 
      93           0 :   for (unsigned int i = 0; i < len - 1; ++i)
      94           0 :     sum += range(i + offset) * (domain(i + 1) - domain(i));
      95             : 
      96           0 :   return _scale_factor * sum;
      97             : }
      98             : 
      99             : KOKKOS_FUNCTION inline Real
     100           0 : KokkosPiecewiseConstant::average() const
     101             : {
     102           0 :   return integral() / (domain(functionSize() - 1) - domain(0));
     103             : }