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 "LBMBoundaryCondition.h"
      10             : #include "LatticeBoltzmannProblem.h"
      11             : #include "LatticeBoltzmannStencilBase.h"
      12             : 
      13             : InputParameters
      14           0 : LBMBoundaryCondition::validParams()
      15             : {
      16           0 :   InputParameters params = LatticeBoltzmannOperator::validParams();
      17           0 :   MooseEnum boundary("top bottom left right front back wall");
      18           0 :   params.addRequiredParam<MooseEnum>(
      19             :       "boundary", boundary, "Edges/Faces where boundary condition is applied.");
      20           0 :   params.addClassDescription("LBMBoundaryCondition object.");
      21           0 :   return params;
      22           0 : }
      23             : 
      24           0 : LBMBoundaryCondition::LBMBoundaryCondition(const InputParameters & parameters)
      25             :   : LatticeBoltzmannOperator(parameters),
      26           0 :     _grid_size(_lb_problem.getGridSize()),
      27           0 :     _boundary(getParam<MooseEnum>("boundary").getEnum<Boundary>())
      28             : {
      29             :   /**
      30             :    * Nodes that are adjacent to boundary will be set to 2, this will later be used in determining
      31             :    * the nodes for bounce-back
      32             :    * This will be achieved by shifting the mesh around in streaming directions and finding where
      33             :    * boundary hit happens
      34             :    */
      35             :   // if (_lb_problem.isBinaryMedia())
      36             :   // {
      37             :   //   auto new_mesh = _mesh.getBinaryMesh().clone();
      38             :   //   for (int64_t ic = 1; ic < _stencil._q; ic++)
      39             :   //   {
      40             :   //     int64_t ex = _stencil._ex[ic].item<int64_t>();
      41             :   //     int64_t ey = _stencil._ey[ic].item<int64_t>();
      42             :   //     int64_t ez = _stencil._ez[ic].item<int64_t>();
      43             :   //     torch::Tensor shifted_mesh = torch::roll(new_mesh, {ex, ey, ez}, {0, 1, 2});
      44             :   //     torch::Tensor adjacent_to_boundary = (shifted_mesh == 0) & (new_mesh == 1);
      45             :   //     new_mesh.masked_fill_(adjacent_to_boundary, 2);
      46             :   //   }
      47             :   //   // Deep copy new mesh
      48             :   //   // MooseTensor::printField(new_mesh, 1, 0);
      49             :   //   _mesh.setBinaryMesh(new_mesh);
      50             :   // }
      51             : 
      52             :   // if (_lb_problem.isBinaryMedia() && _domain.getDim() != 3)
      53             :   //   LBMBoundaryCondition::buildBoundaryIndices();
      54           0 : }
      55             : 
      56             : #if 0
      57             : int64_t
      58             : LBMBoundaryCondition::countNumberofBoundaries()
      59             : {
      60             :   /**
      61             :    * For efficiency, we count the number of boundaries first
      62             :    */
      63             : 
      64             :   LBMBoundaryCondition::determineBoundaryTypes();
      65             : 
      66             :   int64_t k = 0;
      67             :   int64_t num_of_boundaries = 0;
      68             :   for (int64_t i = 0; i < _shape_q[0]; i++)
      69             :     for (int64_t j = 0; j < _shape_q[1]; j++)
      70             :       for (int64_t ic = 0; ic < _shape_q[3]; ic++)
      71             :       {
      72             :         // Avoid calling item() repeatedly
      73             :         int64_t boundary_type = _boundary_types[i][j][k].item<int64_t>();
      74             : 
      75             :         if (boundary_type != -1)
      76             :         {
      77             :           int64_t if_stream_index = _all_boundary_types.size(0) * ic + boundary_type;
      78             : 
      79             :           // when streaming along ic is NOT possible at i, j, k i.e zero
      80             :           if (_if_stream[if_stream_index].item<int64_t>() == 0)
      81             :           {
      82             :             num_of_boundaries++;
      83             :           }
      84             :         }
      85             :       }
      86             : 
      87             :   return num_of_boundaries;
      88             : }
      89             : 
      90             : void
      91             : LBMBoundaryCondition::buildBoundaryIndices()
      92             : {
      93             :   /**
      94             :    * Building boundary indices
      95             :    */
      96             :   // const torch::Tensor & mesh_expanded =
      97             :   // _mesh.getBinaryMesh().unsqueeze(3).expand(expected_shape); auto mask = (mesh_expanded == 2) &
      98             :   // (_u == 0); _boundary_indices = torch::nonzero(mask); _boundary_indices =
      99             :   // _boundary_indices.to(MooseTensor::intTensorOptions());
     100             : 
     101             :   int64_t num_of_boundaries = countNumberofBoundaries();
     102             : 
     103             :   // initialize boundary indices
     104             :   _boundary_indices = torch::zeros({num_of_boundaries, 4}, MooseTensor::intTensorOptions());
     105             : 
     106             :   int64_t row_index = 0;
     107             :   int64_t k = 0;
     108             :   for (int64_t i = 0; i < _shape_q[0]; i++)
     109             :     for (int64_t j = 0; j < _shape_q[1]; j++)
     110             :       for (int64_t ic = 0; ic < _shape_q[3]; ic++)
     111             :       {
     112             :         // Avoid calling item() repeatedly
     113             :         int64_t boundary_type = _boundary_types[i][j][k].item<int64_t>();
     114             : 
     115             :         if (boundary_type != -1)
     116             :         {
     117             :           int64_t if_stream_index = _all_boundary_types.size(0) * ic + boundary_type;
     118             : 
     119             :           // when streaming along ic is NOT possible at i, j, k i.e zero
     120             :           if (_if_stream[if_stream_index].item<int64_t>() == 0)
     121             :           {
     122             :             _boundary_indices[row_index] = torch::tensor(
     123             :                 {i, j, k, _stencil._op[ic].item<int64_t>()}, MooseTensor::intTensorOptions());
     124             :             row_index++;
     125             :           }
     126             :         }
     127             :       }
     128             : }
     129             : 
     130             : void
     131             : LBMBoundaryCondition::determineBoundaryTypes()
     132             : {
     133             :   /**
     134             :    * Scan the binary domain in a 3x3 window to determine boundary types
     135             :    * D2Q9 only
     136             :    */
     137             : 
     138             :   const torch::Tensor & binary_mesh = _lb_problem.getBinaryMedia();
     139             :   _boundary_types =
     140             :       torch::zeros({_shape[0], _shape[1], _shape[2]}, MooseTensor::intTensorOptions());
     141             : 
     142             :   _boundary_types.fill_(-1);
     143             : 
     144             :   // changes the order of D2Q9
     145             :   std::vector<int64_t> d2q9_order_to_new_order{7, 3, 6, 4, 0, 2, 8, 1, 5};
     146             : 
     147             :   int64_t k = 0;
     148             : 
     149             :   for (int64_t i = 0; i < _shape[0]; i++)
     150             :     for (int64_t j = 0; j < _shape[1]; j++)
     151             :     {
     152             :       if (binary_mesh[i][j][k].item<int64_t>() != 0)
     153             :       {
     154             :         // assembling binary string
     155             :         std::string string_of_binary_digits;
     156             :         // search in all directions
     157             :         for (int64_t ic = 0; ic < _stencil._q; ic++)
     158             :         {
     159             :           int64_t i_prime = i + _stencil._ex[d2q9_order_to_new_order[ic]].item<int64_t>();
     160             :           int64_t j_prime = j + _stencil._ey[d2q9_order_to_new_order[ic]].item<int64_t>();
     161             : 
     162             :           // ensuring periodicity
     163             :           i_prime = (i_prime < 0) ? i_prime + _shape[0] : i_prime;
     164             :           i_prime = (i_prime >= _shape[0]) ? i_prime - _shape[0] : i_prime;
     165             : 
     166             :           j_prime = (j_prime < 0) ? j_prime + _shape[1] : j_prime;
     167             :           j_prime = (j_prime >= _shape[1]) ? j_prime - _shape[1] : j_prime;
     168             : 
     169             :           //
     170             :           if (binary_mesh[i_prime][j_prime][k].item<int64_t>() == 0)
     171             :             string_of_binary_digits += '0';
     172             :           else
     173             :             string_of_binary_digits += '1';
     174             :         }
     175             : 
     176             :         // convert binary to decimal
     177             :         std::bitset<32> binary(string_of_binary_digits);
     178             :         int64_t decimal_number = binary.to_ulong();
     179             :         if (decimal_number != 511)
     180             :         {
     181             :           auto comparison_result = (_all_boundary_types == decimal_number);
     182             :           bool is_decimal_number_in_boundary_types = torch::any(comparison_result).item<bool>();
     183             : 
     184             :           if (!is_decimal_number_in_boundary_types)
     185             :           {
     186             :             std::string error_message =
     187             :                 "Boundary type " + std::to_string(decimal_number) + " is not found";
     188             :             mooseError(error_message);
     189             :           }
     190             :           else
     191             :           {
     192             :             // the index where boundary type matches decimal_number
     193             :             int64_t index = torch::nonzero(comparison_result).item<int64_t>();
     194             : 
     195             :             _boundary_types[i][j][k] = index;
     196             :           }
     197             :         }
     198             :       }
     199             :     }
     200             : }
     201             : #endif
     202             : 
     203             : void
     204           0 : LBMBoundaryCondition::computeBuffer()
     205             : {
     206           0 :   switch (_boundary)
     207             :   {
     208           0 :     case Boundary::top:
     209           0 :       topBoundary();
     210           0 :       break;
     211           0 :     case Boundary::bottom:
     212           0 :       bottomBoundary();
     213           0 :       break;
     214           0 :     case Boundary::left:
     215           0 :       leftBoundary();
     216           0 :       break;
     217           0 :     case Boundary::right:
     218           0 :       rightBoundary();
     219           0 :       break;
     220           0 :     case Boundary::front:
     221           0 :       frontBoundary();
     222           0 :       break;
     223           0 :     case Boundary::back:
     224           0 :       backBoundary();
     225           0 :       break;
     226           0 :     case Boundary::wall:
     227           0 :       wallBoundary();
     228           0 :       break;
     229           0 :     default:
     230           0 :       mooseError("Undefined boundary names");
     231             :   }
     232           0 :   _lb_problem.maskedFillSolids(_u, 0);
     233           0 : }