Line data    Source code

       1             : /********************************************************************/
       2             : /*                  SOFTWARE COPYRIGHT NOTIFICATION                 */
       3             : /*                             Cardinal                             */
       4             : /*                                                                  */
       5             : /*                  (c) 2021 UChicago Argonne, LLC                  */
       6             : /*                        ALL RIGHTS RESERVED                       */
       7             : /*                                                                  */
       8             : /*                 Prepared by UChicago Argonne, LLC                */
       9             : /*               Under Contract No. DE-AC02-06CH11357               */
      10             : /*                With the U. S. Department of Energy               */
      11             : /*                                                                  */
      12             : /*             Prepared by Battelle Energy Alliance, LLC            */
      13             : /*               Under Contract No. DE-AC07-05ID14517               */
      14             : /*                With the U. S. Department of Energy               */
      15             : /*                                                                  */
      16             : /*                 See LICENSE for full restrictions                */
      17             : /********************************************************************/
      18             : 
      19             : #ifdef ENABLE_NEK_COUPLING
      20             : 
      21             : #include "NekBinnedPlaneIntegral.h"
      22             : #include "NekInterface.h"
      23             : 
      24             : registerMooseObject("CardinalApp", NekBinnedPlaneIntegral);
      25             : 
      26             : InputParameters
      27          18 : NekBinnedPlaneIntegral::validParams()
      28             : {
      29          18 :   InputParameters params = NekPlaneSpatialBinUserObject::validParams();
      30          18 :   params.addClassDescription(
      31             :       "Compute the spatially-binned side integral of a field over the NekRS mesh");
      32          18 :   return params;
      33           0 : }
      34             : 
      35          34 : NekBinnedPlaneIntegral::NekBinnedPlaneIntegral(const InputParameters & parameters)
      36          34 :   : NekPlaneSpatialBinUserObject(parameters)
      37             : {
      38          33 :   if (_fixed_mesh)
      39          33 :     computeBinVolumes();
      40          32 : }
      41             : 
      42             : void
      43          33 : NekBinnedPlaneIntegral::getBinVolumes()
      44             : {
      45          33 :   resetPartialStorage();
      46          33 :   mesh_t * mesh = nekrs::entireMesh();
      47          33 :   const auto & vgeo = nekrs::getVgeo();
      48             : 
      49       52629 :   for (int k = 0; k < mesh->Nelements; ++k)
      50             :   {
      51       52596 :     int offset = k * mesh->Np;
      52    41592276 :     for (int v = 0; v < mesh->Np; ++v)
      53             :     {
      54    41539680 :       Point p = nekPoint(k, v);
      55             :       unsigned int gap_bin;
      56             :       double distance;
      57    41539680 :       gapIndexAndDistance(p, gap_bin, distance);
      58             : 
      59    41539680 :       if (distance < _gap_thickness / 2.0)
      60             :       {
      61    13583196 :         unsigned int b = bin(p);
      62    13583196 :         _bin_partial_values[b] += vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
      63    13583196 :         _bin_partial_counts[b]++;
      64             :       }
      65             :     }
      66             :   }
      67             : 
      68             :   // sum across all processes
      69          33 :   MPI_Allreduce(
      70             :       _bin_partial_values, _bin_volumes, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
      71          33 :   MPI_Allreduce(
      72             :       _bin_partial_counts, _bin_counts, _n_bins, MPI_INT, MPI_SUM, platform->comm.mpiComm);
      73             : 
      74       27385 :   for (unsigned int i = 0; i < _n_bins; ++i)
      75             :   {
      76             :     // some bins require dividing by a different value, depending on the bin type
      77       27352 :     const auto local_bins = unrolledBin(i);
      78       27352 :     _bin_volumes[i] *= _side_bin->adjustBinValue(local_bins[_side_index]);
      79             : 
      80             :     // dimensionalize
      81       27352 :     nekrs::dimensionalizeVolume(_bin_volumes[i]);
      82       27352 :   }
      83          33 : }
      84             : 
      85             : void
      86          56 : NekBinnedPlaneIntegral::binnedPlaneIntegral(const field::NekFieldEnum & integrand,
      87             :                                             double * total_integral)
      88             : {
      89          56 :   resetPartialStorage();
      90             : 
      91          56 :   mesh_t * mesh = nekrs::entireMesh();
      92          56 :   double (*f)(int, int) = nekrs::solutionPointer(integrand);
      93          56 :   const auto & vgeo = nekrs::getVgeo();
      94             : 
      95       91064 :   for (int k = 0; k < mesh->Nelements; ++k)
      96             :   {
      97       91008 :     int offset = k * mesh->Np;
      98    58140864 :     for (int v = 0; v < mesh->Np; ++v)
      99             :     {
     100    58049856 :       Point p = nekPoint(k, v);
     101             : 
     102             :       unsigned int gap_bin;
     103             :       double distance;
     104    58049856 :       gapIndexAndDistance(p, gap_bin, distance);
     105             : 
     106    58049856 :       if (distance < _gap_thickness / 2.0)
     107             :       {
     108    16288416 :         unsigned int b = bin(p);
     109    16288416 :         _bin_partial_values[b] +=
     110    16288416 :             f(offset + v, 0) * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
     111             :       }
     112             :     }
     113             :   }
     114             : 
     115             :   // sum across all processes
     116          56 :   MPI_Allreduce(
     117             :       _bin_partial_values, total_integral, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
     118             : 
     119        5760 :   for (unsigned int i = 0; i < _n_bins; ++i)
     120             :   {
     121             :     // some bins require dividing by a different value, depending on the bin type
     122        5704 :     const auto local_bins = unrolledBin(i);
     123        5704 :     total_integral[i] *= _side_bin->adjustBinValue(local_bins[_side_index]);
     124             : 
     125        5704 :     nekrs::dimensionalizeVolumeIntegral(integrand, _bin_volumes[i], total_integral[i]);
     126        5704 :   }
     127          56 : }
     128             : 
     129             : Real
     130       62856 : NekBinnedPlaneIntegral::spatialValue(const Point & p, const unsigned int & component) const
     131             : {
     132             :   // total bin index
     133       62856 :   const auto & i = bin(p);
     134             : 
     135             :   // get the index of the gap
     136       62856 :   auto local_indices = unrolledBin(i);
     137       62856 :   auto gap_index = local_indices[_side_index];
     138             : 
     139       62856 :   return _bin_values[i] * _velocity_bin_directions[gap_index](component);
     140       62856 : }
     141             : 
     142             : void
     143          32 : NekBinnedPlaneIntegral::computeIntegral()
     144             : {
     145             :   // if the mesh is changing, re-compute the areas of the bins
     146          32 :   if (!_fixed_mesh)
     147           0 :     computeBinVolumes();
     148             : 
     149          32 :   if (_field == field::velocity_component)
     150             :   {
     151          12 :     binnedPlaneIntegral(field::velocity_x, _bin_values_x);
     152          12 :     binnedPlaneIntegral(field::velocity_y, _bin_values_y);
     153          12 :     binnedPlaneIntegral(field::velocity_z, _bin_values_z);
     154             : 
     155        1188 :     for (unsigned int i = 0; i < num_bins(); ++i)
     156             :     {
     157        1176 :       auto local_indices = unrolledBin(i);
     158        1176 :       auto gap_index = local_indices[_side_index];
     159             : 
     160        1176 :       Point velocity(_bin_values_x[i], _bin_values_y[i], _bin_values_z[i]);
     161        1176 :       _bin_values[i] = _velocity_bin_directions[gap_index] * velocity;
     162        1176 :     }
     163             :   }
     164             :   else
     165          20 :     binnedPlaneIntegral(_field, _bin_values);
     166          32 : }
     167             : 
     168             : void
     169           8 : NekBinnedPlaneIntegral::executeUserObject()
     170             : {
     171           8 :   computeIntegral();
     172             : 
     173             :   // correct the values to areas
     174        1088 :   for (unsigned int i = 0; i < _n_bins; ++i)
     175        1080 :     _bin_values[i] /= _gap_thickness;
     176           8 : }
     177             : 
     178             : #endif