NekRSProblem.h

Go to the documentation of this file.

/********************************************************************/
/*                  SOFTWARE COPYRIGHT NOTIFICATION                 */
/*                             Cardinal                             */
/*                                                                  */
/*                  (c) 2021 UChicago Argonne, LLC                  */
/*                        ALL RIGHTS RESERVED                       */
/*                                                                  */
/*                 Prepared by UChicago Argonne, LLC                */
/*               Under Contract No. DE-AC02-06CH11357               */
/*                With the U. S. Department of Energy               */
/*                                                                  */
/*             Prepared by Battelle Energy Alliance, LLC            */
/*               Under Contract No. DE-AC07-05ID14517               */
/*                With the U. S. Department of Energy               */
/*                                                                  */
/*                 See LICENSE for full restrictions                */
/********************************************************************/

#pragma once

#include "CardinalProblem.h"
#include "NekInterface.h"
#include "NekTimeStepper.h"
#include "NekScalarValue.h"
#include "NekRSMesh.h"
#include "Transient.h"

class FieldTransferBase;

class NekRSProblem : public CardinalProblem
{
public:
  NekRSProblem(const InputParameters & params);

  static InputParameters validParams();

  ~NekRSProblem();

  int nPoints() const { return _n_points; }

  bool isDataTransferHappening(ExternalProblem::Direction direction);

  Transient * transientExecutioner() const { return _transient_executioner; }

  bool isUsrWrkSlotReservedForCoupling(const unsigned int & slot) const;

  void copyIndividualScratchSlot(const unsigned int & slot) const;

  void writeFieldFile(const Real & time, const int & step) const;

  virtual bool isOutputStep() const;

  virtual void initialSetup() override;

  virtual void externalSolve() override;

  virtual bool converged(unsigned int) override { return true; }

  virtual void addExternalVariables() override;

  virtual void syncSolutions(ExternalProblem::Direction direction) override;

  virtual bool nondimensional() const { return _nondimensional; }

  virtual const bool hasMovingNekMesh() const { return false; }

  virtual bool synchronizeIn();

  virtual bool synchronizeOut();

  unsigned int nUsrWrkSlots() const { return _n_usrwrk_slots; }

  template <typename T>
  void volumeSolution(const T & field, double * s)
  {
    mesh_t * mesh = nekrs::entireMesh();
    auto vc = _nek_mesh->volumeCoupling();

    double (*f)(int, int);
    f = nekrs::solutionPointer(field);

    int start_1d = mesh->Nq;
    int end_1d = _moose_Nq;
    int start_3d = start_1d * start_1d * start_1d;
    int end_3d = end_1d * end_1d * end_1d;
    int n_to_write = vc.n_elems * end_3d * _nek_mesh->nBuildPerVolumeElem();

    // allocate temporary space:
    // - Ttmp: results of the search for each process
    // - Telem: scratch space for volume interpolation to avoid reallocating a bunch (only used if
    // interpolating)
    double * Ttmp = (double *)calloc(n_to_write, sizeof(double));
    double * Telem = (double *)calloc(start_3d, sizeof(double));
    auto indices = _nek_mesh->cornerIndices();

    int c = 0;
    for (int k = 0; k < mesh->Nelements; ++k)
    {
      int offset = k * start_3d;

      for (int build = 0; build < _nek_mesh->nBuildPerVolumeElem(); ++build)
      {
        if (_needs_interpolation)
        {
          // get the solution on the element
          for (int v = 0; v < start_3d; ++v)
            Telem[v] = f(offset + v, 0 /* unused for volumes */);

          // and then interpolate it
          nekrs::interpolateVolumeHex3D(
              _interpolation_outgoing, Telem, start_1d, &(Ttmp[c]), end_1d);
          c += end_3d;
        }
        else
        {
          // get the solution on the element - no need to interpolate
          for (int v = 0; v < end_3d; ++v, ++c)
            Ttmp[c] = f(offset + indices[build][v], 0 /* unused for volumes */);
        }
      }
    }

    // dimensionalize the solution if needed
    for (int v = 0; v < n_to_write; ++v)
      Ttmp[v] =
          Ttmp[v] * nekrs::nondimensionalDivisor(field) + nekrs::nondimensionalAdditive(field);

    nekrs::allgatherv(vc.mirror_counts, Ttmp, s, end_3d);

    freePointer(Ttmp);
    freePointer(Telem);
  }

  template <typename T>
  void boundarySolution(const T & field, double * s)
  {
    mesh_t * mesh = nekrs::entireMesh();
    auto bc = _nek_mesh->boundaryCoupling();

    double (*f)(int, int);
    f = nekrs::solutionPointer(field);

    int start_1d = mesh->Nq;
    int end_1d = _moose_Nq;
    int start_2d = start_1d * start_1d;
    int end_2d = end_1d * end_1d;

    int n_to_write = bc.n_faces * end_2d * _nek_mesh->nBuildPerSurfaceElem();

    // allocate temporary space:
    // - Ttmp: results of the search for each process
    // - Tface: scratch space for face solution to avoid reallocating a bunch (only used if
    // interpolating)
    // - scratch: scratch for the interpolatino process to avoid reallocating a bunch (only used if
    // interpolating0
    double * Ttmp = (double *)calloc(n_to_write, sizeof(double));
    double * Tface = (double *)calloc(start_2d, sizeof(double));
    double * scratch = (double *)calloc(start_1d * end_1d, sizeof(double));

    auto indices = _nek_mesh->cornerIndices();

    int c = 0;
    for (int k = 0; k < bc.total_n_faces; ++k)
    {
      if (bc.process[k] == nekrs::commRank())
      {
        int i = bc.element[k];
        int j = bc.face[k];
        int offset = i * mesh->Nfaces * start_2d + j * start_2d;

        for (int build = 0; build < _nek_mesh->nBuildPerSurfaceElem(); ++build)
        {
          if (_needs_interpolation)
          {
            // get the solution on the face
            for (int v = 0; v < start_2d; ++v)
            {
              int id = mesh->vmapM[offset + v];
              Tface[v] = f(id, mesh->Nsgeo * (offset + v));
            }

            // and then interpolate it
            nekrs::interpolateSurfaceFaceHex3D(
                scratch, _interpolation_outgoing, Tface, start_1d, &(Ttmp[c]), end_1d);
            c += end_2d;
          }
          else
          {
            // get the solution on the face - no need to interpolate
            for (int v = 0; v < end_2d; ++v, ++c)
            {
              int id = mesh->vmapM[offset + indices[build][v]];
              Ttmp[c] = f(id, mesh->Nsgeo * (offset + v));
            }
          }
        }
      }
    }

    // dimensionalize the solution if needed
    for (int v = 0; v < n_to_write; ++v)
      Ttmp[v] =
          Ttmp[v] * nekrs::nondimensionalDivisor(field) + nekrs::nondimensionalAdditive(field);

    nekrs::allgatherv(bc.mirror_counts, Ttmp, s, end_2d);

    freePointer(Ttmp);
    freePointer(Tface);
    freePointer(scratch);
  }

  template <typename T>
  void writeVolumeSolution(const int elem_id,
                           const T & field,
                           double * s,
                           const std::vector<double> * add = nullptr)
  {
    mesh_t * mesh = nekrs::entireMesh();
    void (*write_solution)(int, dfloat);
    write_solution = nekrs::solutionWritePointer(field);

    auto vc = _nek_mesh->volumeCoupling();
    int id = vc.element[elem_id] * mesh->Np;

    if (_nek_mesh->exactMirror())
    {
      // can write directly into the NekRS solution
      for (int v = 0; v < mesh->Np; ++v)
      {
        double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
        write_solution(id + v, s[v] + extra);
      }
    }
    else
    {
      // need to interpolate onto the higher-order Nek mesh
      double * tmp = (double *)calloc(mesh->Np, sizeof(double));

      interpolateVolumeSolutionToNek(elem_id, s, tmp);

      for (int v = 0; v < mesh->Np; ++v)
      {
        double extra = (add == nullptr) ? 0.0 : (*add)[id + v];
        write_solution(id + v, tmp[v] + extra);
      }

      freePointer(tmp);
    }
  }

  template <typename T>
  void writeBoundarySolution(const int elem_id, const T & field, double * s)
  {
    mesh_t * mesh = nekrs::temperatureMesh();
    void (*write_solution)(int, dfloat);
    write_solution = nekrs::solutionWritePointer(field);

    const auto & bc = _nek_mesh->boundaryCoupling();
    int offset = bc.element[elem_id] * mesh->Nfaces * mesh->Nfp + bc.face[elem_id] * mesh->Nfp;

    if (_nek_mesh->exactMirror())
    {
      // can write directly into the NekRS solution
      for (int i = 0; i < mesh->Nfp; ++i)
        write_solution(mesh->vmapM[offset + i], s[i]);
    }
    else
    {
      // need to interpolate onto the higher-order Nek mesh
      double * tmp = (double *)calloc(mesh->Nfp, sizeof(double));
      interpolateBoundarySolutionToNek(s, tmp);

      for (int i = 0; i < mesh->Nfp; ++i)
        write_solution(mesh->vmapM[offset + i], tmp[i]);

      freePointer(tmp);
    }
  }

  std::string casename() const { return _casename; }

  void mapFaceDataToNekFace(const unsigned int & e,
                            const unsigned int & var_num,
                            const Real & divisor,
                            const Real & additive,
                            double ** outgoing_data);

  void mapVolumeDataToNekVolume(const unsigned int & e,
                                const unsigned int & var_num,
                                const Real & divisor,
                                const Real & additive,
                                double ** outgoing_data);

  void mapFaceDataToNekVolume(const unsigned int & e,
                              const unsigned int & var_num,
                              const Real & divisor,
                              const Real & additive,
                              double ** outgoing_data);

protected:
  void copyScratchToDevice();

  void interpolateBoundarySolutionToNek(double * incoming_moose_value, double * outgoing_nek_value);

  void interpolateVolumeSolutionToNek(const int elem_id,
                                      double * incoming_moose_value,
                                      double * outgoing_nek_value);

  void initializeInterpolationMatrices();

  std::unique_ptr<NumericVector<Number>> _serialized_solution;

  std::string fieldFilePrefix(const int & number) const;

  const std::string & _casename;

  const bool & _write_fld_files;

  const bool & _disable_fld_file_output;

  bool _nondimensional;

  unsigned int _n_usrwrk_slots;

  const unsigned int & _constant_interval;

  const bool & _skip_final_field_file;

  int _n_surface_elems;

  int _n_vertices_per_surface;

  int _n_volume_elems;

  int _n_vertices_per_volume;

  double _start_time;

  bool _is_output_step;

  NekRSMesh * _nek_mesh;

  NekTimeStepper * _timestepper = nullptr;

  Transient * _transient_executioner = nullptr;

  bool _needs_interpolation;

  int _n_points;

  const PostprocessorValue * _transfer_in = nullptr;

  double * _interpolation_outgoing = nullptr;

  double * _interpolation_incoming = nullptr;

  int _moose_Nq;

  const std::vector<unsigned int> * _usrwrk_output = nullptr;

  const std::vector<std::string> * _usrwrk_output_prefix = nullptr;

  double _elapsedStepSum;

  double _elapsedTime;

  double _tSolveStepMin;

  double _tSolveStepMax;

  synchronization::SynchronizationEnum _sync_interval;

  static bool _first;

  std::vector<FieldTransferBase *> _field_transfers;

  std::vector<ScalarTransferBase *> _scalar_transfers;

  std::set<unsigned int> _usrwrk_slots;
};