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          50 : NekBinnedPlaneIntegral::NekBinnedPlaneIntegral(const InputParameters & parameters)
      36          50 :   : NekPlaneSpatialBinUserObject(parameters)
      37             : {
      38          49 :   if (_fixed_mesh)
      39          49 :     computeBinVolumes();
      40          48 : }
      41             : 
      42             : void
      43          49 : NekBinnedPlaneIntegral::getBinVolumes()
      44             : {
      45          49 :   resetPartialStorage();
      46          49 :   mesh_t * mesh = nekrs::entireMesh();
      47          49 :   const auto & vgeo = nekrs::getVgeo();
      48             : 
      49       66613 :   for (int k = 0; k < mesh->Nelements; ++k)
      50             :   {
      51       66564 :     int offset = k * mesh->Np;
      52    48757860 :     for (int v = 0; v < mesh->Np; ++v)
      53             :     {
      54    48691296 :       Point p = nekPoint(k, v);
      55             :       unsigned int gap_bin;
      56             :       double distance;
      57    48691296 :       gapIndexAndDistance(p, gap_bin, distance);
      58             : 
      59    48691296 :       if (distance < _gap_thickness / 2.0)
      60             :       {
      61    13949340 :         unsigned int b = bin(p);
      62    13949340 :         _bin_partial_values[b] += vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
      63    13949340 :         _bin_partial_counts[b]++;
      64             :       }
      65             :     }
      66             :   }
      67             : 
      68             :   // sum across all processes
      69          49 :   MPI_Allreduce(
      70             :       _bin_partial_values, _bin_volumes, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
      71          49 :   MPI_Allreduce(
      72             :       _bin_partial_counts, _bin_counts, _n_bins, MPI_INT, MPI_SUM, platform->comm.mpiComm);
      73             : 
      74       30089 :   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       30040 :     const auto local_bins = unrolledBin(i);
      78       30040 :     _bin_volumes[i] *= _side_bin->adjustBinValue(local_bins[_side_index]);
      79             : 
      80             :     // dimensionalize
      81       30040 :     nekrs::dimensionalizeVolume(_bin_volumes[i]);
      82             :   }
      83          49 : }
      84             : 
      85             : void
      86          88 : NekBinnedPlaneIntegral::binnedPlaneIntegral(const field::NekFieldEnum & integrand,
      87             :                                             double * total_integral)
      88             : {
      89          88 :   resetPartialStorage();
      90             : 
      91          88 :   mesh_t * mesh = nekrs::entireMesh();
      92          88 :   double (*f)(int) = nekrs::solutionPointer(integrand);
      93          88 :   const auto & vgeo = nekrs::getVgeo();
      94             : 
      95      119032 :   for (int k = 0; k < mesh->Nelements; ++k)
      96             :   {
      97      118944 :     int offset = k * mesh->Np;
      98    72472032 :     for (int v = 0; v < mesh->Np; ++v)
      99             :     {
     100    72353088 :       Point p = nekPoint(k, v);
     101             : 
     102             :       unsigned int gap_bin;
     103             :       double distance;
     104    72353088 :       gapIndexAndDistance(p, gap_bin, distance);
     105             : 
     106    72353088 :       if (distance < _gap_thickness / 2.0)
     107             :       {
     108    17020704 :         unsigned int b = bin(p);
     109    17020704 :         _bin_partial_values[b] += f(offset + v) * vgeo[mesh->Nvgeo * offset + v + mesh->Np * JWID];
     110             :       }
     111             :     }
     112             :   }
     113             : 
     114             :   // sum across all processes
     115          88 :   MPI_Allreduce(
     116             :       _bin_partial_values, total_integral, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
     117             : 
     118       11168 :   for (unsigned int i = 0; i < _n_bins; ++i)
     119             :   {
     120             :     // some bins require dividing by a different value, depending on the bin type
     121       11080 :     const auto local_bins = unrolledBin(i);
     122       11080 :     total_integral[i] *= _side_bin->adjustBinValue(local_bins[_side_index]);
     123             : 
     124       11080 :     nekrs::dimensionalizeVolumeIntegral(integrand, _bin_volumes[i], total_integral[i]);
     125             :   }
     126          88 : }
     127             : 
     128             : Real
     129       83808 : NekBinnedPlaneIntegral::spatialValue(const Point & p, const unsigned int & component) const
     130             : {
     131             :   // total bin index
     132       83808 :   const auto & i = bin(p);
     133             : 
     134             :   // get the index of the gap
     135       83808 :   auto local_indices = unrolledBin(i);
     136       83808 :   auto gap_index = local_indices[_side_index];
     137             : 
     138       83808 :   return _bin_values[i] * _velocity_bin_directions[gap_index](component);
     139             : }
     140             : 
     141             : void
     142          48 : NekBinnedPlaneIntegral::computeIntegral()
     143             : {
     144             :   // if the mesh is changing, re-compute the areas of the bins
     145          48 :   if (!_fixed_mesh)
     146           0 :     computeBinVolumes();
     147             : 
     148          48 :   if (_field == field::velocity_component)
     149             :   {
     150          20 :     binnedPlaneIntegral(field::velocity_x, _bin_values_x);
     151          20 :     binnedPlaneIntegral(field::velocity_y, _bin_values_y);
     152          20 :     binnedPlaneIntegral(field::velocity_z, _bin_values_z);
     153             : 
     154        2540 :     for (unsigned int i = 0; i < num_bins(); ++i)
     155             :     {
     156        2520 :       auto local_indices = unrolledBin(i);
     157        2520 :       auto gap_index = local_indices[_side_index];
     158             : 
     159        2520 :       Point velocity(_bin_values_x[i], _bin_values_y[i], _bin_values_z[i]);
     160        2520 :       _bin_values[i] = _velocity_bin_directions[gap_index] * velocity;
     161             :     }
     162             :   }
     163             :   else
     164          28 :     binnedPlaneIntegral(_field, _bin_values);
     165          48 : }
     166             : 
     167             : void
     168           8 : NekBinnedPlaneIntegral::executeUserObject()
     169             : {
     170           8 :   computeIntegral();
     171             : 
     172             :   // correct the values to areas
     173        1088 :   for (unsigned int i = 0; i < _n_bins; ++i)
     174        1080 :     _bin_values[i] /= _gap_thickness;
     175           8 : }
     176             : 
     177             : #endif