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 "NekPointValue.h"
      22             : #include "pointInterpolation.hpp"
      23             : 
      24             : registerMooseObject("CardinalApp", NekPointValue);
      25             : 
      26             : InputParameters
      27         441 : NekPointValue::validParams()
      28             : {
      29         441 :   InputParameters params = NekFieldPostprocessor::validParams();
      30         882 :   params.addRequiredParam<Point>("point", "The physical point where the field will be evaluated");
      31         441 :   params.addClassDescription("Uses NekRS's pointInterpolation to query the NekRS solution at a "
      32             :                              "point (does not need to be a grid point).");
      33         441 :   return params;
      34           0 : }
      35             : 
      36         147 : NekPointValue::NekPointValue(const InputParameters & parameters)
      37         294 :   : NekFieldPostprocessor(parameters), _point(getParam<Point>("point")), _value(0)
      38             : {
      39         147 : }
      40             : 
      41             : void
      42         147 : NekPointValue::execute()
      43             : {
      44         147 :   std::vector<dfloat> x = {_point(0) / nekrs::referenceLength()};
      45         147 :   std::vector<dfloat> y = {_point(1) / nekrs::referenceLength()};
      46         147 :   std::vector<dfloat> z = {_point(2) / nekrs::referenceLength()};
      47         147 :   int n = x.size();
      48             : 
      49         147 :   nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
      50             : 
      51             :   // set the points to be interpolated; we include this every time we call the
      52             :   // interpolator, in case the mesh is moving. TODO: auto-detect for efficiency
      53         147 :   auto interp = pointInterpolation_t(nrs);
      54         147 :   interp.setPoints(n, x.data(), y.data(), z.data());
      55             :   const auto verbosity = pointInterpolation_t::VerbosityLevel::Basic;
      56         147 :   interp.find(verbosity);
      57             : 
      58             :   // if this slot is not used for coupling, we are responsible for copying it to
      59             :   // device before calling the interp function. Otherwise, Cardinal handles copying
      60             :   // to device automatically.
      61         147 :   if (_field == field::usrwrk00 && !_nek_problem->isUsrWrkSlotReservedForCoupling(0))
      62           9 :     _nek_problem->copyIndividualScratchSlot(0);
      63         147 :   if (_field == field::usrwrk01 && !_nek_problem->isUsrWrkSlotReservedForCoupling(1))
      64           9 :     _nek_problem->copyIndividualScratchSlot(1);
      65         147 :   if (_field == field::usrwrk02 && !_nek_problem->isUsrWrkSlotReservedForCoupling(2))
      66           9 :     _nek_problem->copyIndividualScratchSlot(2);
      67             : 
      68             :   // interpolate the field onto those points
      69         147 :   occa::memory o_interpolated;
      70             :   int n_values = n;
      71         147 :   switch (_field)
      72             :   {
      73          64 :     case field::velocity_component:
      74             :     case field::velocity_x:
      75             :     case field::velocity_y:
      76             :     case field::velocity_z:
      77             :     case field::velocity:
      78             :     case field::velocity_x_squared:
      79             :     case field::velocity_y_squared:
      80             :     case field::velocity_z_squared:
      81          64 :       n_values = n * nrs->NVfields;
      82          64 :       o_interpolated = platform->device.malloc<dfloat>(n_values);
      83          64 :       interp.eval(n_values, nrs->fieldOffset, nrs->cds->o_U, n, o_interpolated);
      84             :       break;
      85           8 :     case field::pressure:
      86           8 :       o_interpolated = platform->device.malloc<dfloat>(n);
      87           8 :       interp.eval(1, nrs->fieldOffset, nrs->o_P, n, o_interpolated);
      88             :       break;
      89          32 :     case field::temperature:
      90             :     case field::scalar01:
      91             :     case field::scalar02:
      92             :     case field::scalar03:
      93          32 :       n_values = n * nrs->Nscalar;
      94          32 :       o_interpolated = platform->device.malloc<dfloat>(n_values);
      95          32 :       interp.eval(n_values, nekrs::scalarFieldOffset(), nrs->cds->o_S, n, o_interpolated);
      96             :       break;
      97             :     case field::unity:
      98             :       break;
      99          35 :     case field::usrwrk00:
     100             :     case field::usrwrk01:
     101             :     case field::usrwrk02:
     102          35 :       n_values = n * _nek_problem->nUsrWrkSlots();
     103          35 :       o_interpolated = platform->device.malloc<dfloat>(n_values);
     104          35 :       interp.eval(n_values, nekrs::fieldOffset(), nrs->o_usrwrk, n, o_interpolated);
     105             :       break;
     106           0 :     default:
     107           0 :       mooseError("Unhandled NekFieldEnum in NekPointValue!");
     108             :   }
     109             : 
     110             :   // the interpolation happens on device, so we need to copy it back to the host
     111         147 :   std::vector<dfloat> interpolated(n_values);
     112         147 :   o_interpolated.copyTo(interpolated.data(), n_values);
     113             : 
     114             :   // because of NekRS's way of storing solution, we need extra steps to actually
     115             :   // return what the user wants
     116         147 :   switch (_field)
     117             :   {
     118             :     case field::velocity_component:
     119           8 :       _value = interpolated[0] * _velocity_direction(0) + interpolated[1] * _velocity_direction(1) +
     120           8 :                interpolated[2] * _velocity_direction(2);
     121           8 :       break;
     122             :     case field::velocity_x:
     123           8 :       _value = interpolated[0];
     124           8 :       break;
     125             :     case field::velocity_y:
     126           8 :       _value = interpolated[1];
     127           8 :       break;
     128             :     case field::velocity_z:
     129           8 :       _value = interpolated[2];
     130           8 :       break;
     131             :     case field::velocity:
     132           8 :       _value = std::sqrt(interpolated[0] * interpolated[0] + interpolated[1] * interpolated[1] +
     133           8 :                          interpolated[2] * interpolated[2]);
     134           8 :       break;
     135             :     case field::velocity_x_squared:
     136           8 :       _value = interpolated[0] * interpolated[0];
     137           8 :       break;
     138             :     case field::velocity_y_squared:
     139           8 :       _value = interpolated[1] * interpolated[1];
     140           8 :       break;
     141             :     case field::velocity_z_squared:
     142           8 :       _value = interpolated[2] * interpolated[2];
     143           8 :       break;
     144             :     case field::pressure:
     145           8 :       _value = interpolated[0];
     146           8 :       break;
     147             :     case field::temperature:
     148           8 :       _value = interpolated[0];
     149           8 :       break;
     150             :     case field::scalar01:
     151           8 :       _value = interpolated[1];
     152           8 :       break;
     153             :     case field::scalar02:
     154           8 :       _value = interpolated[2];
     155           8 :       break;
     156             :     case field::scalar03:
     157           8 :       _value = interpolated[3];
     158           8 :       break;
     159           8 :     case field::unity:
     160           8 :       _value = 1;
     161           8 :       break;
     162             :     case field::usrwrk00:
     163           9 :       _value = interpolated[0];
     164           9 :       break;
     165             :     case field::usrwrk01:
     166          13 :       _value = interpolated[1];
     167          13 :       break;
     168             :     case field::usrwrk02:
     169          13 :       _value = interpolated[2];
     170          13 :       break;
     171           0 :     default:
     172           0 :       mooseError("Unhandled NekFieldEnum in NekPointValue!");
     173             :   }
     174             : 
     175         147 :   _value = _value * nekrs::nondimensionalDivisor(_field) + nekrs::nondimensionalAdditive(_field);
     176         288 : }
     177             : 
     178             : Real
     179         144 : NekPointValue::getValue() const
     180             : {
     181         144 :   return _value;
     182             : }
     183             : 
     184             : #endif