RhieChowMassFlux.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include "RhieChowFaceFluxProvider.h"
#include "CellCenteredMapFunctor.h"
#include "FaceCenteredMapFunctor.h"
#include "VectorComponentFunctor.h"
#include <unordered_map>
#include <set>
#include <unordered_set>
 
#include "libmesh/petsc_vector.h"
 
class MooseMesh;
class INSFVVelocityVariable;
class INSFVPressureVariable;
class LinearFVPressureCorrectionDiffusion;
class LinearFVGradientReader;
namespace libMesh
{
class Elem;
class MeshBase;
}
 
class RhieChowMassFlux : public RhieChowFaceFluxProvider, public NonADFunctorInterface
{
public:
  static InputParameters validParams();
  RhieChowMassFlux(const InputParameters & params);
 
  Real getMassFlux(const FaceInfo & fi) const;
 
  Real getVolumetricFaceFlux(const FaceInfo & fi) const;
 
  const LinearFVGradientReader & pressureGradientField() const;
 
  virtual Real getVolumetricFaceFlux(const Moose::FV::InterpMethod m,
                                     const FaceInfo & fi,
                                     const Moose::StateArg & time,
                                     const THREAD_ID tid,
                                     bool subtract_mesh_velocity) const override;
 
  void initFaceMassFlux();
  void initCouplingField();
  void computeFaceMassFlux();
  void computeCellVelocity();
 
  virtual void meshChanged() override;
  virtual void initialize() override;
  virtual void execute() override {}
  virtual void finalize() override {}
  virtual void initialSetup() override;
 
  void linkMomentumPressureSystems(const std::vector<LinearSystem *> & momentum_systems,
                                   LinearSystem & pressure_system,
                                   const std::vector<unsigned int> & momentum_system_numbers);
 
  void computeHbyA(const bool with_updated_pressure, const bool verbose);
 
protected:
  std::vector<std::unique_ptr<NumericVector<Number>>> &
  selectPressureGradient(const bool updated_pressure);
 
  const std::vector<std::unique_ptr<NumericVector<Number>>> & pressureGradientComponents() const;
 
  void setupMeshInformation();
 
  void
  populateCouplingFunctors(const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_hbya,
                           const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_Ainv);
 
  template <typename VarType>
  void checkBlocks(const VarType & var) const;
 
  virtual bool supportMeshVelocity() const override { return false; }
 
  const MooseMesh & _moose_mesh;
 
  const libMesh::MeshBase & _mesh;
 
  const unsigned int _dim;
 
  const MooseLinearVariableFVReal * const _p;
 
  std::vector<const MooseLinearVariableFVReal *> _vel;
 
  LinearFVPressureCorrectionDiffusion * _p_diffusion_kernel;
 
  FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> _HbyA_flux;
 
  std::vector<std::unique_ptr<NumericVector<Number>>> _HbyA_raw;
 
  FaceCenteredMapFunctor<RealVectorValue, std::unordered_map<dof_id_type, RealVectorValue>> _Ainv;
 
  std::vector<std::unique_ptr<NumericVector<Number>>> _Ainv_raw;
 
  std::unique_ptr<NumericVector<Number>> _A_avg;
 
  FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> & _face_mass_flux;
 
  std::vector<std::unique_ptr<NumericVector<Number>>> _grad_p_current;
 
  const Moose::Functor<Real> & _rho;
 
  std::vector<LinearSystem *> _momentum_systems;
 
  std::vector<unsigned int> _momentum_system_numbers;
 
  std::vector<unsigned int> _global_momentum_system_numbers;
 
  std::vector<libMesh::LinearImplicitSystem *> _momentum_implicit_systems;
 
  LinearSystem * _pressure_system;
 
  const LinearFVGradientReader * _pressure_gradient_field;
 
  unsigned int _global_pressure_system_number;
 
  std::unique_ptr<NumericVector<Number>> _cell_volumes;
 
  const MooseEnum _pressure_projection_method;
 
  const Moose::FV::InterpMethod _pressure_diffusion_interp_method;
 
private:
  std::vector<const FaceInfo *> _flow_face_info;
};
 
template <typename VarType>
void
RhieChowMassFlux::checkBlocks(const VarType & var) const
{
  const auto & var_blocks = var.blockIDs();
  const auto & uo_blocks = blockIDs();
 
  // Error if this UO has any blocks that the variable does not
  std::set<SubdomainID> uo_blocks_minus_var_blocks;
  std::set_difference(uo_blocks.begin(),
                      uo_blocks.end(),
                      var_blocks.begin(),
                      var_blocks.end(),
                      std::inserter(uo_blocks_minus_var_blocks, uo_blocks_minus_var_blocks.end()));
  if (uo_blocks_minus_var_blocks.size() > 0)
    mooseError("Block restriction of interpolator user object '",
               this->name(),
               "' (",
               Moose::stringify(blocks()),
               ") includes blocks not in the block restriction of variable '",
               var.name(),
               "' (",
               Moose::stringify(var.blocks()),
               ")");
 
  // Get the blocks in the variable but not this UO
  std::set<SubdomainID> var_blocks_minus_uo_blocks;
  std::set_difference(var_blocks.begin(),
                      var_blocks.end(),
                      uo_blocks.begin(),
                      uo_blocks.end(),
                      std::inserter(var_blocks_minus_uo_blocks, var_blocks_minus_uo_blocks.end()));
 
  // For each block in the variable but not this UO, error if there is connection
  // to any blocks on the UO.
  for (auto & block_id : var_blocks_minus_uo_blocks)
  {
    const auto connected_blocks = _moose_mesh.getBlockConnectedBlocks(block_id);
    std::set<SubdomainID> connected_blocks_on_uo;
    std::set_intersection(connected_blocks.begin(),
                          connected_blocks.end(),
                          uo_blocks.begin(),
                          uo_blocks.end(),
                          std::inserter(connected_blocks_on_uo, connected_blocks_on_uo.end()));
    if (connected_blocks_on_uo.size() > 0)
      mooseError("Block restriction of interpolator user object '",
                 this->name(),
                 "' (",
                 Moose::stringify(uo_blocks),
                 ") doesn't match the block restriction of variable '",
                 var.name(),
                 "' (",
                 Moose::stringify(var_blocks),
                 ")");
  }
}