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);
  }
 
  void writeVolumeDisplacement(const int elem_id,
                               double * s,
                               const field::NekWriteEnum f,
                               const std::vector<double> * add = nullptr);
  void writeVolumeSolution(const int slot,
                           const int elem_id,
                           double * s,
                           const std::vector<double> * add = nullptr);
 
  void writeBoundarySolution(const int slot, const int elem_id, double * s);
 
  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 copyHostToDevice();
 
  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;
};