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 "LBMFixedFirstOrderBC.h"
      10             : #include "LatticeBoltzmannProblem.h"
      11             : #include "LatticeBoltzmannStencilBase.h"
      12             : 
      13             : #include <cstdlib>
      14             : 
      15             : using namespace torch::indexing;
      16             : 
      17             : registerMooseObject("SwiftApp", LBMFixedFirstOrderBC);
      18             : 
      19             : InputParameters
      20           0 : LBMFixedFirstOrderBC::validParams()
      21             : {
      22           0 :   InputParameters params = LBMBoundaryCondition::validParams();
      23           0 :   params.addClassDescription("LBMFixedFirstOrderBC object");
      24           0 :   params.addRequiredParam<TensorInputBufferName>("f", "Input buffer distribution function");
      25           0 :   params.addRequiredParam<std::string>("value", "Fixed input velocity");
      26           0 :   params.addParam<bool>("perturb", false, "Whether to perturb first order moment at the boundary");
      27           0 :   return params;
      28           0 : }
      29             : 
      30           0 : LBMFixedFirstOrderBC::LBMFixedFirstOrderBC(const InputParameters & parameters)
      31             :   : LBMBoundaryCondition(parameters),
      32           0 :     _f(getInputBufferByName(getParam<TensorInputBufferName>("f"))),
      33           0 :     _grid_size(_lb_problem.getGridSize()),
      34           0 :     _value(_lb_problem.getConstant<Real>(getParam<std::string>("value"))),
      35           0 :     _perturb(getParam<bool>("perturb"))
      36             : {
      37           0 : }
      38             : 
      39             : void
      40           0 : LBMFixedFirstOrderBC::frontBoundary()
      41             : {
      42           0 :   if (_domain.getDim() == 2)
      43           0 :     mooseError("There is no front boundary in 2 dimensions.");
      44             :   else
      45           0 :     mooseError("Front boundary is not implemented, but it can be replaced by any other boundary by "
      46             :                "rotating the domain.");
      47             : }
      48             : 
      49             : void
      50           0 : LBMFixedFirstOrderBC::backBoundary()
      51             : {
      52           0 :   if (_domain.getDim() == 2)
      53           0 :     mooseError("There is no back boundary in 2 dimensions.");
      54             :   else
      55           0 :     mooseError("Back boundary is not implemented, but it can be replaced by any other boundary by "
      56             :                "rotating the domain.");
      57             : }
      58             : 
      59             : void
      60           0 : LBMFixedFirstOrderBC::leftBoundaryD2Q9()
      61             : {
      62             :   Real deltaU = 0.0;
      63           0 :   torch::Tensor u_x_perturbed = torch::zeros({_grid_size[1], 1}, MooseTensor::floatTensorOptions());
      64             : 
      65           0 :   if (_perturb)
      66             :   {
      67           0 :     deltaU = 1.0e-6 * _value;
      68             :     Real phi = 0.0; // static_cast<Real>(rand()) / static_cast<float>(RAND_MAX) * 2.0 * M_PI;
      69             :     torch::Tensor y_coords =
      70           0 :         torch::arange(0, _grid_size[1], MooseTensor::floatTensorOptions()).unsqueeze(1);
      71           0 :     u_x_perturbed = _value + deltaU * torch::sin(y_coords / _grid_size[1] * 2.0 * M_PI + phi);
      72             :   }
      73             :   else
      74           0 :     u_x_perturbed.fill_(_value);
      75             : 
      76             :   torch::Tensor density =
      77           0 :       1.0 / (1.0 - u_x_perturbed) *
      78           0 :       (_f.index({0, Slice(), Slice(), 0}) + _f.index({0, Slice(), Slice(), 2}) +
      79           0 :        _f.index({0, Slice(), Slice(), 4}) +
      80           0 :        2.0 * (_f.index({0, Slice(), Slice(), 3}) + _f.index({0, Slice(), Slice(), 6}) +
      81           0 :               _f.index({0, Slice(), Slice(), 7})));
      82             : 
      83             :   // axis aligned direction
      84           0 :   const auto & opposite_dir = _stencil._op[_stencil._left[0]];
      85           0 :   _u.index_put_({0, Slice(), Slice(), _stencil._left[0]},
      86           0 :                 _f.index({0, Slice(), Slice(), opposite_dir}) +
      87           0 :                     2.0 / 3.0 * density * u_x_perturbed);
      88             : 
      89             :   // other directions
      90           0 :   for (unsigned int i = 1; i < _stencil._left.size(0); i++)
      91             :   {
      92           0 :     const auto & opposite_dir = _stencil._op[_stencil._left[i]];
      93           0 :     _u.index_put_(
      94           0 :         {0, Slice(), Slice(), _stencil._left[i]},
      95           0 :         _f.index({0, Slice(), Slice(), opposite_dir}) -
      96           0 :             0.5 * _stencil._ey[_stencil._left[i]] *
      97           0 :                 (_f.index({0, Slice(), Slice(), 2}) - _f.index({0, Slice(), Slice(), 4})) +
      98           0 :             1.0 / 6.0 * density * u_x_perturbed);
      99             :   }
     100           0 : }
     101             : 
     102             : void
     103           0 : LBMFixedFirstOrderBC::leftBoundary()
     104             : {
     105           0 :   if (_stencil._q == 9)
     106           0 :     leftBoundaryD2Q9(); // higher order specialization for D2Q9
     107             :   else
     108             :   {
     109           0 :     torch::Tensor density = 1.0 / (1.0 - _value) *
     110           0 :                             (torch::sum(_f.index({0, Slice(), Slice(), -_stencil._neutral_x}), -1) +
     111           0 :                              2 * torch::sum(_f.index({0, Slice(), Slice(), _stencil._right}), -1));
     112             : 
     113           0 :     _u.index_put_({0, Slice(), Slice(), _stencil._left[0]},
     114           0 :                   _f.index({0, Slice(), Slice(), _stencil._right[0]}) +
     115           0 :                       2.0 * _stencil._weights[_stencil._left[0]] / _lb_problem._cs2 * _value *
     116             :                           density);
     117             : 
     118           0 :     for (unsigned int i = 1; i < _stencil._left.size(0); i++)
     119             :     {
     120           0 :       _u.index_put_({0, Slice(), Slice(), _stencil._left[i]},
     121           0 :                     _f.index({0, Slice(), Slice(), _stencil._right[i]}) +
     122           0 :                         2.0 * _stencil._weights[_stencil._left[i]] / _lb_problem._cs2 * _value *
     123             :                             density);
     124             :     }
     125             :   }
     126           0 : }
     127             : 
     128             : void
     129           0 : LBMFixedFirstOrderBC::rightBoundaryD2Q9()
     130             : {
     131           0 :   torch::Tensor density = 1.0 / (1.0 + _value) *
     132           0 :                           (_f.index({_grid_size[0] - 1, Slice(), Slice(), 0}) +
     133           0 :                            _f.index({_grid_size[0] - 1, Slice(), Slice(), 2}) +
     134           0 :                            _f.index({_grid_size[0] - 1, Slice(), Slice(), 4}) +
     135           0 :                            2 * (_f.index({_grid_size[0] - 1, Slice(), Slice(), 1}) +
     136           0 :                                 _f.index({_grid_size[0] - 1, Slice(), Slice(), 5}) +
     137           0 :                                 _f.index({_grid_size[0] - 1, Slice(), Slice(), 8})));
     138             : 
     139             :   // axis aligned direction
     140           0 :   const auto & opposite_dir = _stencil._op[_stencil._left[0]];
     141           0 :   _u.index_put_({_grid_size[0] - 1, Slice(), Slice(), opposite_dir},
     142           0 :                 _f.index({_grid_size[0] - 1, Slice(), Slice(), _stencil._left[0]}) -
     143           0 :                     2.0 / 3.0 * density * _value);
     144             : 
     145             :   // other directions
     146           0 :   for (unsigned int i = 1; i < _stencil._left.size(0); i++)
     147             :   {
     148           0 :     const auto & opposite_dir = _stencil._op[_stencil._left[i]];
     149           0 :     _u.index_put_({_grid_size[0] - 1, Slice(), Slice(), opposite_dir},
     150           0 :                   _f.index({_grid_size[0] - 1, Slice(), Slice(), _stencil._left[i]}) +
     151           0 :                       0.5 * _stencil._ey[opposite_dir] *
     152           0 :                           (_f.index({_grid_size[0] - 1, Slice(), Slice(), 4}) -
     153           0 :                            _f.index({_grid_size[0] - 1, Slice(), Slice(), 2})) -
     154           0 :                       1.0 / 6.0 * density * _value);
     155             :   }
     156           0 : }
     157             : 
     158             : void
     159           0 : LBMFixedFirstOrderBC::rightBoundary()
     160             : {
     161           0 :   if (_stencil._q == 9)
     162           0 :     rightBoundaryD2Q9(); // higher order specialization for D2Q9
     163             :   else
     164             :   {
     165             :     torch::Tensor density =
     166           0 :         1.0 / (1.0 + _value) *
     167           0 :         (torch::sum(_f.index({_grid_size[0] - 1, Slice(), Slice(), -_stencil._neutral_x}), -1) +
     168           0 :          2 * torch::sum(_f.index({_grid_size[0] - 1, Slice(), Slice(), _stencil._left}), -1));
     169             : 
     170           0 :     _u.index_put_({_grid_size[0] - 1, Slice(), Slice(), _stencil._right[0]},
     171           0 :                   _f.index({_grid_size[0] - 1, Slice(), Slice(), _stencil._left[0]}) -
     172           0 :                       2.0 * _stencil._weights[_stencil._right[0]] / _lb_problem._cs2 * _value *
     173             :                           density);
     174             : 
     175           0 :     for (unsigned int i = 1; i < _stencil._right.size(0); i++)
     176             :     {
     177           0 :       _u.index_put_({_grid_size[0] - 1, Slice(), Slice(), _stencil._right[i]},
     178           0 :                     _f.index({_grid_size[0] - 1, Slice(), Slice(), _stencil._left[i]}) -
     179           0 :                         2.0 * _stencil._weights[_stencil._right[i]] / _lb_problem._cs2 * _value *
     180             :                             density);
     181             :     }
     182             :   }
     183           0 : }
     184             : 
     185             : void
     186           0 : LBMFixedFirstOrderBC::bottomBoundaryD2Q9()
     187             : {
     188             :   torch::Tensor density =
     189           0 :       1.0 / (1.0 - _value) *
     190           0 :       (_f.index({Slice(), 0, Slice(), 0}) + _f.index({Slice(), 0, Slice(), 1}) +
     191           0 :        _f.index({Slice(), 0, Slice(), 3}) +
     192           0 :        2 * (_f.index({Slice(), 0, Slice(), 4}) + _f.index({Slice(), 0, Slice(), 7}) +
     193           0 :             _f.index({Slice(), 0, Slice(), 8})));
     194             : 
     195             :   // axis aligned direction
     196           0 :   const auto & opposite_dir = _stencil._op[_stencil._bottom[0]];
     197           0 :   _u.index_put_({Slice(), 0, Slice(), _stencil._bottom[0]},
     198           0 :                 _f.index({Slice(), 0, Slice(), opposite_dir}) + 2.0 / 3.0 * density * _value);
     199             : 
     200             :   // other directions
     201           0 :   for (unsigned int i = 1; i < _stencil._bottom.size(0); i++)
     202             :   {
     203           0 :     const auto & opposite_dir = _stencil._op[_stencil._bottom[i]];
     204           0 :     _u.index_put_(
     205           0 :         {Slice(), 0, Slice(), _stencil._bottom[i]},
     206           0 :         _f.index({Slice(), 0, Slice(), opposite_dir}) -
     207           0 :             0.5 * _stencil._ex[_stencil._bottom[i]] *
     208           0 :                 (_f.index({Slice(), 0, Slice(), 1}) - _f.index({Slice(), 0, Slice(), 3})) +
     209           0 :             1.0 / 6.0 * density * _value);
     210             :   }
     211           0 : }
     212             : 
     213             : void
     214           0 : LBMFixedFirstOrderBC::bottomBoundary()
     215             : {
     216           0 :   if (_stencil._q == 9)
     217           0 :     bottomBoundaryD2Q9();
     218             :   else
     219           0 :     mooseError("Bottom boundary is not implemented, but it can be replaced by another boundary by "
     220             :                "rotating the domain.");
     221           0 : }
     222             : 
     223             : void
     224           0 : LBMFixedFirstOrderBC::topBoundaryD2Q9()
     225             : {
     226           0 :   torch::Tensor density = 1.0 / (1.0 + _value) *
     227           0 :                           (_f.index({Slice(), _grid_size[1] - 1, Slice(), 0}) +
     228           0 :                            _f.index({Slice(), _grid_size[1] - 1, Slice(), 1}) +
     229           0 :                            _f.index({Slice(), _grid_size[1] - 1, Slice(), 3}) +
     230           0 :                            2 * (_f.index({Slice(), _grid_size[1] - 1, Slice(), 2}) +
     231           0 :                                 _f.index({Slice(), _grid_size[1] - 1, Slice(), 5}) +
     232           0 :                                 _f.index({Slice(), _grid_size[1] - 1, Slice(), 6})));
     233             : 
     234             :   // axis aligned direction
     235           0 :   const auto & opposite_dir = _stencil._op[_stencil._bottom[0]];
     236           0 :   _u.index_put_({Slice(), _grid_size[1] - 1, Slice(), opposite_dir},
     237           0 :                 _f.index({Slice(), _grid_size[1] - 1, Slice(), _stencil._bottom[0]}) -
     238           0 :                     2.0 / 3.0 * density * _value);
     239             : 
     240             :   // other directions
     241           0 :   for (unsigned int i = 1; i < _stencil._bottom.size(0); i++)
     242             :   {
     243           0 :     const auto & opposite_dir = _stencil._op[_stencil._bottom[i]];
     244           0 :     _u.index_put_({Slice(), _grid_size[1] - 1, Slice(), opposite_dir},
     245           0 :                   _f.index({Slice(), _grid_size[1] - 1, Slice(), _stencil._bottom[i]}) +
     246           0 :                       0.5 * _stencil._ex[opposite_dir] *
     247           0 :                           (_f.index({Slice(), _grid_size[1] - 1, Slice(), 3}) -
     248           0 :                            _f.index({Slice(), _grid_size[1] - 1, Slice(), 1})) -
     249           0 :                       1.0 / 6.0 * density * _value);
     250             :   }
     251           0 : }
     252             : 
     253             : void
     254           0 : LBMFixedFirstOrderBC::topBoundary()
     255             : {
     256           0 :   if (_stencil._q == 9)
     257           0 :     topBoundaryD2Q9();
     258             :   else
     259           0 :     mooseError("Top boundary is not implemented, but it can be replaced by another boundary by "
     260             :                "rotating the domain.");
     261           0 : }
     262             : 
     263             : void
     264           0 : LBMFixedFirstOrderBC::computeBuffer()
     265             : {
     266           0 :   _u = _u.clone();
     267           0 :   switch (_boundary)
     268             :   {
     269           0 :     case Boundary::top:
     270           0 :       topBoundary();
     271           0 :       break;
     272           0 :     case Boundary::bottom:
     273           0 :       bottomBoundary();
     274           0 :       break;
     275           0 :     case Boundary::left:
     276           0 :       leftBoundary();
     277           0 :       break;
     278           0 :     case Boundary::right:
     279           0 :       rightBoundary();
     280           0 :       break;
     281           0 :     case Boundary::front:
     282           0 :       frontBoundary();
     283           0 :       break;
     284           0 :     case Boundary::back:
     285           0 :       backBoundary();
     286           0 :       break;
     287           0 :     case Boundary::wall:
     288           0 :       wallBoundary();
     289           0 :       break;
     290           0 :     default:
     291           0 :       mooseError("Undefined boundary names");
     292             :   }
     293           0 :   _lb_problem.maskedFillSolids(_u, 0);
     294           0 : }