Line data    Source code

       1             : /**********************************************************************/
       2             : /*                    DO NOT MODIFY THIS HEADER                       */
       3             : /*             Swift, a Fourier spectral solver for MOOSE             */
       4             : /*                                                                    */
       5             : /*            Copyright 2024 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : 
       9             : #include "SmoothRectangleCompute.h"
      10             : // #include "SwiftUtils.h"
      11             : #include "TensorProblem.h"
      12             : #include <torch/torch.h>
      13             : 
      14             : registerMooseObject("SwiftApp", SmoothRectangleCompute);
      15             : 
      16             : InputParameters
      17           0 : SmoothRectangleCompute::validParams()
      18             : {
      19           0 :   InputParameters params = TensorOperator::validParams();
      20           0 :   params.addClassDescription(
      21             :       "Interpolate a value between the inside and outside of a rectangle smoothly.");
      22           0 :   params.addRequiredParam<Real>("x1", "The x coordinate of the lower left-hand corner of the box.");
      23           0 :   params.addRequiredParam<Real>("x2",
      24             :                                 "The x coordinate of the upper right-hand corner of the box.");
      25           0 :   params.addRequiredParam<Real>("y1", "The x coordinate of the lower left-hand corner of the box.");
      26           0 :   params.addRequiredParam<Real>("y2",
      27             :                                 "The x coordinate of the upper right-hand corner of the box.");
      28           0 :   params.addParam<Real>("z1", 0, "The z coordinate of the lower left-hand corner of the box.");
      29           0 :   params.addParam<Real>("z2", 0, "The z coordinate of the upper right-hand corner of the box.");
      30           0 :   MooseEnum interpolationFunction("COS TANH");
      31           0 :   params.addParam<MooseEnum>(
      32             :       "profile", interpolationFunction, "Functional dependence for the interface profile");
      33           0 :   params.addParam<Real>(
      34           0 :       "int_width", 0, "The width of the diffuse interface. Set to 0 for sharp interface.");
      35           0 :   params.addParam<Real>("inside", 1, "The value inside the rectangle.");
      36           0 :   params.addParam<Real>("outside", 0, "The value outside the rectangle.");
      37           0 :   return params;
      38           0 : }
      39             : 
      40           0 : SmoothRectangleCompute::SmoothRectangleCompute(const InputParameters & parameters)
      41             :   : TensorOperator(parameters),
      42           0 :     _x1(getParam<Real>("x1")),
      43           0 :     _x2(getParam<Real>("x2")),
      44           0 :     _y1(getParam<Real>("y1")),
      45           0 :     _y2(getParam<Real>("y2")),
      46           0 :     _z1(getParam<Real>("z1")),
      47           0 :     _z2(getParam<Real>("z2")),
      48           0 :     _interp_func(getParam<MooseEnum>("profile").getEnum<interpolationFunction>()),
      49           0 :     _int_width(getParam<Real>("int_width")),
      50           0 :     _inside(getParam<Real>("inside")),
      51           0 :     _outside(getParam<Real>("outside"))
      52             : {
      53           0 :   if (_int_width < 0.0)
      54           0 :     mooseError("Interface width must be a non-negative real number.");
      55           0 : }
      56             : 
      57             : void
      58           0 : SmoothRectangleCompute::computeBuffer()
      59             : {
      60           0 :   auto dim = _domain.getDim();
      61           0 :   auto h_box = torch::zeros(_tensor_problem.getShape(), torch::kDouble);
      62             : 
      63             :   // sharp interpolation
      64           0 :   if (_int_width <= 0.0)
      65             :   {
      66           0 :     auto h_x = (_x >= _x1) & (_x <= _x2);
      67           0 :     auto h_y = (dim >= 2) ? ((_y >= _y1) & (_y <= _y2)) : torch::ones_like(_y, torch::kBool);
      68           0 :     auto h_z = (dim == 3) ? ((_z >= _z1) & (_z <= _z2)) : torch::ones_like(_z, torch::kBool);
      69             :     // Reshape to the same dimensions for logical_and operation
      70           0 :     h_x = h_x.reshape({-1, 1, 1});
      71           0 :     h_y = h_y.reshape({1, -1, 1});
      72           0 :     h_z = h_z.reshape({1, 1, -1});
      73             : 
      74             :     // Combine the conditions
      75           0 :     auto combined_conditions = torch::logical_and(h_x, torch::logical_and(h_y, h_z)).squeeze();
      76             :     // Apply the combined conditions to h_box
      77           0 :     h_box.index_put_({combined_conditions}, 1.0);
      78             :   }
      79             :   else
      80             :   {
      81             :     // generate distances based on the right problem dimension
      82             : 
      83           0 :     switch (_interp_func)
      84             :     {
      85           0 :       case interpolationFunction::COS:
      86             :       {
      87           0 :         auto min_x = torch::minimum(_x - _x1, _x2 - _x).clamp(-_int_width / 2.0, _int_width / 2.0);
      88             :         auto min_y =
      89             :             (dim >= 2)
      90           0 :                 ? torch::minimum(_y - _y1, _y2 - _y).clamp(-_int_width / 2.0, _int_width / 2.0)
      91           0 :                 : torch::full_like(_y, _int_width / 2.0);
      92             :         auto min_z =
      93             :             (dim == 3)
      94           0 :                 ? torch::minimum(_z - _z1, _z2 - _z).clamp(-_int_width / 2.0, _int_width / 2.0)
      95           0 :                 : torch::full_like(_z, _int_width / 2.0);
      96           0 :         auto h_x = 0.5 + 0.5 * torch::sin(pi * min_x / _int_width);
      97           0 :         auto h_y = 0.5 + 0.5 * torch::sin(pi * min_y / _int_width);
      98           0 :         auto h_z = 0.5 + 0.5 * torch::sin(pi * min_z / _int_width);
      99           0 :         h_box = h_x.reshape({-1, 1, 1}) * h_y.reshape({1, -1, 1}) * h_z.reshape({1, 1, -1});
     100             :         break;
     101             :       }
     102           0 :       case interpolationFunction::TANH:
     103             :       {
     104           0 :         auto min_x = torch::minimum(_x - _x1, _x2 - _x);
     105             :         auto min_y =
     106           0 :             (dim >= 2) ? torch::minimum(_y - _y1, _y2 - _y) : torch::full_like(_y, 10 * _int_width);
     107           0 :         auto min_z = (dim == 3) ? torch::minimum(_z - _z1, _z2 - _z)
     108           0 :                                 : torch::full_like(_z, 10 * _int_width / 2.0);
     109           0 :         auto h_x = 0.5 + 0.5 * torch::tanh(4 * min_x / _int_width);
     110           0 :         auto h_y = 0.5 + 0.5 * torch::tanh(4 * min_y / _int_width);
     111           0 :         auto h_z = 0.5 + 0.5 * torch::tanh(4 * min_z / _int_width);
     112           0 :         h_box = h_x.reshape({-1, 1, 1}) * h_y.reshape({1, -1, 1}) * h_z.reshape({1, 1, -1});
     113             :         break;
     114             :       }
     115             :     }
     116             :   }
     117           0 :   _u = (h_box * _inside + (1 - h_box) * _outside).squeeze();
     118           0 : }