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 "NekBinnedSideIntegral.h"
      22             : #include "NekInterface.h"
      23             : 
      24             : registerMooseObject("CardinalApp", NekBinnedSideIntegral);
      25             : 
      26             : InputParameters
      27         113 : NekBinnedSideIntegral::validParams()
      28             : {
      29         113 :   InputParameters params = NekSideSpatialBinUserObject::validParams();
      30         113 :   params.addClassDescription(
      31             :       "Compute the spatially-binned side integral of a field over a boundary of the NekRS mesh");
      32         113 :   return params;
      33           0 : }
      34             : 
      35          57 : NekBinnedSideIntegral::NekBinnedSideIntegral(const InputParameters & parameters)
      36          57 :   : NekSideSpatialBinUserObject(parameters)
      37             : {
      38          56 :   if (_fixed_mesh)
      39          56 :     computeBinVolumes();
      40          56 : }
      41             : 
      42             : void
      43          56 : NekBinnedSideIntegral::getBinVolumes()
      44             : {
      45          56 :   resetPartialStorage();
      46             : 
      47          56 :   mesh_t * mesh = nekrs::entireMesh();
      48          56 :   const auto & sgeo = nekrs::getSgeo();
      49             : 
      50       92688 :   for (int i = 0; i < mesh->Nelements; ++i)
      51             :   {
      52      648424 :     for (int j = 0; j < mesh->Nfaces; ++j)
      53             :     {
      54      555792 :       int face_id = mesh->EToB[i * mesh->Nfaces + j];
      55      555792 :       if (std::find(_boundary.begin(), _boundary.end(), face_id) != _boundary.end())
      56             :       {
      57       38472 :         int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
      58      625224 :         for (int v = 0; v < mesh->Nfp; ++v)
      59             :         {
      60      586752 :           Point p = nekPoint(i, j, v);
      61      586752 :           unsigned int b = bin(p);
      62      586752 :           _bin_partial_values[b] += sgeo[mesh->Nsgeo * (offset + v) + WSJID];
      63      586752 :           _bin_partial_counts[b]++;
      64             :         }
      65             :       }
      66             :     }
      67             :   }
      68             : 
      69             :   // sum across all processes
      70          56 :   MPI_Allreduce(
      71             :       _bin_partial_values, _bin_volumes, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
      72          56 :   MPI_Allreduce(
      73             :       _bin_partial_counts, _bin_counts, _n_bins, MPI_INT, MPI_SUM, platform->comm.mpiComm);
      74             : 
      75             :   // dimensionalize
      76         272 :   for (unsigned int i = 0; i < _n_bins; ++i)
      77         216 :     nekrs::dimensionalizeArea(_bin_volumes[i]);
      78          56 : }
      79             : 
      80             : void
      81          72 : NekBinnedSideIntegral::binnedSideIntegral(const field::NekFieldEnum & integrand,
      82             :                                           double * total_integral)
      83             : {
      84          72 :   resetPartialStorage();
      85             : 
      86          72 :   mesh_t * mesh = nekrs::entireMesh();
      87          72 :   double (*f)(int) = nekrs::solutionPointer(integrand);
      88          72 :   const auto & sgeo = nekrs::getSgeo();
      89             : 
      90      107936 :   for (int i = 0; i < mesh->Nelements; ++i)
      91             :   {
      92      755048 :     for (int j = 0; j < mesh->Nfaces; ++j)
      93             :     {
      94      647184 :       int face_id = mesh->EToB[i * mesh->Nfaces + j];
      95      647184 :       if (std::find(_boundary.begin(), _boundary.end(), face_id) != _boundary.end())
      96             :       {
      97       39240 :         int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
      98      638280 :         for (int v = 0; v < mesh->Nfp; ++v)
      99             :         {
     100      599040 :           Point p = nekPoint(i, j, v);
     101      599040 :           unsigned int b = bin(p);
     102      599040 :           _bin_partial_values[b] +=
     103      599040 :               f(mesh->vmapM[offset + v]) * sgeo[mesh->Nsgeo * (offset + v) + WSJID];
     104             :         }
     105             :       }
     106             :     }
     107             :   }
     108             : 
     109             :   // sum across all processes
     110          72 :   MPI_Allreduce(
     111             :       _bin_partial_values, total_integral, _n_bins, MPI_DOUBLE, MPI_SUM, platform->comm.mpiComm);
     112             : 
     113             :   // dimensionalize
     114         368 :   for (unsigned int i = 0; i < _n_bins; ++i)
     115         296 :     nekrs::dimensionalizeSideIntegral(integrand, _bin_volumes[i], total_integral[i]);
     116          72 : }
     117             : 
     118             : Real
     119           0 : NekBinnedSideIntegral::spatialValue(const Point & p, const unsigned int & component) const
     120             : {
     121           0 :   const auto & i = bin(p);
     122           0 :   return _bin_values[i] * _velocity_bin_directions[i](component);
     123             : }
     124             : 
     125             : void
     126          72 : NekBinnedSideIntegral::computeIntegral()
     127             : {
     128             :   // if the mesh is changing, re-compute the areas of the bins
     129          72 :   if (!_fixed_mesh)
     130           0 :     computeBinVolumes();
     131             : 
     132          72 :   if (_field == field::velocity_component)
     133             :   {
     134           0 :     binnedSideIntegral(field::velocity_x, _bin_values_x);
     135           0 :     binnedSideIntegral(field::velocity_y, _bin_values_y);
     136           0 :     binnedSideIntegral(field::velocity_z, _bin_values_z);
     137             : 
     138           0 :     for (unsigned int i = 0; i < num_bins(); ++i)
     139             :     {
     140           0 :       Point velocity(_bin_values_x[i], _bin_values_y[i], _bin_values_z[i]);
     141           0 :       _bin_values[i] = _velocity_bin_directions[i] * velocity;
     142             :     }
     143             :   }
     144             :   else
     145          72 :     binnedSideIntegral(_field, _bin_values);
     146          72 : }
     147             : 
     148             : void
     149          72 : NekBinnedSideIntegral::executeUserObject()
     150             : {
     151          72 :   computeIntegral();
     152          72 : }
     153             : 
     154             : #endif