23 #include "libmesh/mesh_base.h" 24 #include "libmesh/elem_range.h" 25 #include "libmesh/petsc_matrix.h" 37 params.addClassDescription(
"Computes H/A and 1/A together with face mass fluxes for segregated " 38 "momentum-pressure equations using linear systems.");
40 params.addRequiredParam<VariableName>(
NS::pressure,
"The pressure variable.");
41 params.addRequiredParam<VariableName>(
"u",
"The x-component of velocity");
42 params.addParam<VariableName>(
"v",
"The y-component of velocity");
43 params.addParam<VariableName>(
"w",
"The z-component of velocity");
44 params.addRequiredParam<std::string>(
"p_diffusion_kernel",
45 "The diffusion kernel acting on the pressure.");
47 params.addRequiredParam<MooseFunctorName>(
NS::density,
"Density functor");
57 params.addParam<
MooseEnum>(
"pressure_projection_method",
58 MooseEnum(
"standard consistent",
"standard"),
59 "The method to use in the pressure projection for Ainv - " 60 "standard (SIMPLE) or consistent (SIMPLEC)");
70 _dim(blocksMaxDimension()),
74 _HbyA_flux(_moose_mesh, blockIDs(),
"HbyA_flux"),
75 _Ainv(_moose_mesh, blockIDs(),
"Ainv"),
78 "face_flux", _moose_mesh, blockIDs(),
"face_values")),
80 _pressure_projection_method(getParam<
MooseEnum>(
"pressure_projection_method"))
86 std::vector<std::string> vel_names = {
"u",
"v",
"w"};
93 paramError(vel_names[i],
"the velocity must be a MOOSELinearVariableFVReal.");
104 if (!dynamic_cast<SIMPLE *>(
getMooseApp().getExecutioner()) &&
105 !dynamic_cast<PIMPLE *>(
getMooseApp().getExecutioner()))
107 " should only be used with a linear segregated thermal-hydraulics solver!");
112 const std::vector<LinearSystem *> & momentum_systems,
114 const std::vector<unsigned int> & momentum_system_numbers)
145 std::vector<LinearFVFluxKernel *> flux_kernel;
148 .template condition<AttribThread>(
_tid)
149 .
template condition<AttribSysNum>(
_p->
sys().
number())
150 .
template condition<AttribSystem>(
"LinearFVFluxKernel")
151 .template condition<AttribName>(getParam<std::string>(
"p_diffusion_kernel"))
152 .queryInto(flux_kernel);
153 if (flux_kernel.size() != 1)
155 "p_diffusion_kernel",
156 "The kernel with the given name could not be found or multiple instances were identified.");
160 "The provided diffusion kernel should of type LinearFVAnisotropicDiffusion!");
172 if (
hasBlocks(elem_info->subdomain_id()))
175 _cell_volumes->set(elem_dof, elem_info->volume() * elem_info->coordFactor());
182 if (
hasBlocks(fi->elemPtr()->subdomain_id()) ||
183 (fi->neighborPtr() &&
hasBlocks(fi->neighborPtr()->subdomain_id())))
193 for (
const auto & pair :
_Ainv)
194 _Ainv[pair.first] = 0;
211 if (
_vel[0]->isInternalFace(*fi))
213 const auto & elem_info = *fi->elemInfo();
214 const auto & neighbor_info = *fi->neighborInfo();
221 density_times_velocity(dim_i),
222 _vel[dim_i]->getElemValue(elem_info, time_arg) * elem_rho,
223 _vel[dim_i]->getElemValue(neighbor_info, time_arg) * neighbor_rho,
230 const bool elem_is_fluid =
hasBlocks(fi->elemPtr()->subdomain_id());
231 const Elem *
const boundary_elem = elem_is_fluid ? fi->elemPtr() : fi->neighborPtr();
234 const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
238 const Real face_rho =
_rho(boundary_face, time_arg);
240 density_times_velocity(dim_i) = boundary_normal_multiplier * face_rho *
272 bool libmesh_dbg_var(subtract_mesh_velocity))
const 274 mooseAssert(!subtract_mesh_velocity,
"RhieChowMassFlux does not support moving meshes yet!");
277 mooseError(
"Interpolation methods other than Rhie-Chow are not supported!");
279 mooseError(
"Older interpolation times are not supported!");
305 Real p_grad_flux = 0.0;
308 const auto & elem_info = *fi->
elemInfo();
316 const auto p_elem_value = p_reader(elem_dof);
317 const auto p_neighbor_value = p_reader(neighbor_dof);
321 const auto neighbor_matrix_contribution =
323 const auto elem_rhs_contribution =
327 p_grad_flux = (p_neighbor_value * neighbor_matrix_contribution +
328 p_elem_value * elem_matrix_contribution) -
329 elem_rhs_contribution;
333 mooseAssert(fi->
boundaryIDs().size() == 1,
"We should only have one boundary on every face.");
335 bc_pointer->setupFaceData(
341 const auto matrix_contribution =
343 const auto rhs_contribution =
347 p_grad_flux = (p_elem_value * matrix_contribution - rhs_contribution);
363 auto working_vector =
_Ainv_raw[system_i]->clone();
364 working_vector->pointwise_mult(*working_vector, *pressure_gradient[system_i]);
365 working_vector->add(*
_HbyA_raw[system_i]);
366 working_vector->scale(-1.0);
398 std::vector<PetscVectorReader> hbya_reader;
400 hbya_reader.emplace_back(*raw_hbya[dim_i]);
402 std::vector<PetscVectorReader> ainv_reader;
404 ainv_reader.emplace_back(*raw_Ainv[dim_i]);
416 if (
_vel[0]->isInternalFace(*fi))
419 const auto & elem_info = *fi->
elemInfo();
435 hbya_reader[dim_i](elem_dof),
436 hbya_reader[dim_i](neighbor_dof),
441 elem_rho * ainv_reader[dim_i](elem_dof),
442 neighbor_rho * ainv_reader[dim_i](neighbor_dof),
452 const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
459 if (
_vel[0]->isDirichletBoundaryFace(*fi))
467 -boundary_normal_multiplier *
477 face_hbya(dim_i) = boundary_normal_multiplier * hbya_reader[dim_i](elem_dof);
483 Ainv(dim_i) = elem_rho * ainv_reader[dim_i](elem_dof);
495 _console <<
"************************************" << std::endl;
496 _console <<
"Computing HbyA" << std::endl;
497 _console <<
"************************************" << std::endl;
500 "The momentum system shall be linked before calling this function!");
515 "The matrices used in the segregated INSFVRhieChow objects need to be convertable " 520 _console <<
"Matrix in rc object" << std::endl;
532 _console <<
"Velocity solution in H(u)" << std::endl;
550 mooseAssert(working_vector_petsc,
551 "The vectors used in the RhieChowMassFlux objects need to be convertable " 556 HbyA.
add(-1.0, *working_vector_petsc);
570 HbyA.
add(-1.0, *working_vector_petsc);
574 _console <<
"total RHS" << std::endl;
576 _console <<
"pressure RHS" << std::endl;
577 pressure_gradient[system_i]->print();
578 _console <<
" H(u)-rhs-relaxation_source" << std::endl;
583 *working_vector_petsc = 1.0;
591 _console <<
" (H(u)-rhs)/A" << std::endl;
603 _console <<
"Performing SIMPLEC projection." << std::endl;
612 const auto local_size = mmat->
local_m();
614 for (
const auto row_i :
make_range(local_size))
617 const auto global_index = mmat->
row_start() + row_i;
618 std::vector<numeric_index_type> indices;
619 std::vector<Real> values;
620 mmat->
get_row(global_index, indices, values);
623 const Real row_sum = std::accumulate(values.cbegin(), values.cend(), 0.0);
626 sum_vector.add(global_index, row_sum);
632 auto row_sum = current_local_solution.
zero_clone();
633 get_row_sum(*row_sum);
636 auto Ainv_full = current_local_solution.
zero_clone();
637 *working_vector_petsc = 1.0;
639 const auto Ainv_full_old = Ainv_full->clone();
642 Ainv_full->add(-1.0, Ainv);
643 working_vector_petsc->
pointwise_mult(*Ainv_full, *pressure_gradient[system_i]);
645 HbyA.
add(-1.0, *working_vector_petsc);
648 Ainv = *Ainv_full_old;
659 _console <<
"************************************" << std::endl;
660 _console <<
"DONE Computing HbyA " << std::endl;
661 _console <<
"************************************" << std::endl;
665 std::vector<std::unique_ptr<NumericVector<Number>>> &
668 if (updated_pressure)
std::vector< LinearSystem * > _momentum_systems
Pointers to the linear system(s) in moose corresponding to the momentum equation(s) ...
const MooseEnum _pressure_projection_method
Enumerator for the method used for pressure projection.
virtual void initialize() override
void setupMeshInformation()
Compute the cell volumes on the mesh.
User object responsible for determining the face fluxes using the Rhie-Chow interpolation in a segreg...
A functor whose evaluation relies on querying a map where the keys are face info ids and the values c...
const std::set< BoundaryID > & boundaryIDs() const
virtual numeric_index_type local_m() const final
T & getMesh(MooseMesh &mesh)
function to cast mesh
void checkBlocks(const VarType &var) const
Check the block consistency between the passed in var and us.
const std::vector< const ElemInfo *> & elemInfoVector() const
std::vector< libMesh::LinearImplicitSystem * > _momentum_implicit_systems
Pointers to the momentum equation implicit system(s) from libmesh.
unsigned int number() const
const std::vector< std::unique_ptr< NumericVector< Number > > > & gradientContainer() const
Real getMassFlux(const FaceInfo &fi) const
Get the face velocity times density (used in advection terms)
void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const
virtual Real computeElemMatrixContribution() override
static const std::string component
const Elem & elem() const
virtual std::unique_ptr< NumericVector< T > > zero_clone() const=0
const ElemInfo * neighborInfo() const
const ExecFlagType EXEC_NONE
const Elem * elem() const
registerMooseObject("NavierStokesApp", RhieChowMassFlux)
void computeFaceMassFlux()
Update the values of the face velocities in the containers.
void addFunctor(const std::string &name, const Moose::FunctorBase< T > &functor, const THREAD_ID tid)
static const std::string density
virtual void setupFaceData(const FaceInfo *face_info)
NumericVector< Number > * rhs
void addAvailableFlags(const ExecFlagType &flag, Args... flags)
const ElemInfo * elemInfo() const
std::vector< std::unique_ptr< NumericVector< Number > > > _Ainv_raw
We hold on to the cell-based 1/A vectors so that we can easily reconstruct the cell velocities as wel...
const LinearSystem * _pressure_system
Pointer to the pressure system.
LinearFVBoundaryCondition * getBoundaryCondition(const BoundaryID bd_id) const
FaceCenteredMapFunctor< RealVectorValue, std::unordered_map< dof_id_type, RealVectorValue > > _Ainv
A map functor from faces to $(1/A)_f$.
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
std::vector< unsigned int > _momentum_system_numbers
Numbers of the momentum system(s)
virtual void initialSetup() override
virtual void meshChanged() override
virtual void pointwise_divide(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
virtual const std::string & name() const
virtual Real computeBoundaryMatrixContribution(const LinearFVBoundaryCondition &bc) override
MooseApp & getMooseApp() const
std::vector< unsigned int > _global_momentum_system_numbers
Global numbers of the momentum system(s)
static InputParameters validParams()
Moose::ElemArg makeElemArg(const Elem *elem, bool correct_skewnewss=false) const
bool isInternalFace(const FaceInfo &) const
void computeHbyA(const bool with_updated_pressure, const bool verbose)
Computes the inverse of the diagonal (1/A) of the system matrix plus the H/A components for the press...
const std::vector< const FaceInfo *> & faceInfo() const
ValueType evaluate(const FaceInfo *const fi) const
Evaluate the face functor using a FaceInfo argument.
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
FaceCenteredMapFunctor< Real, std::unordered_map< dof_id_type, Real > > _HbyA_flux
A map functor from faces to $HbyA_{ij} = (A_{offdiag}*{(predicted~velocity)} - {Source})_{ij}/A_{ij}$...
const Elem * neighborPtr() const
virtual const NumericVector< Number > *const & currentSolution() const override final
TheWarehouse & theWarehouse() const
std::vector< const MooseLinearVariableFVReal * > _vel
The thread 0 copy of the x-velocity variable.
const MooseLinearVariableFVReal *const _p
The thread 0 copy of the pressure variable.
unsigned int _global_pressure_system_number
Global number of the pressure system.
RhieChowMassFlux(const InputParameters ¶ms)
void linkMomentumPressureSystems(const std::vector< LinearSystem *> &momentum_systems, const LinearSystem &pressure_system, const std::vector< unsigned int > &momentum_system_numbers)
Update the momentum system-related information.
const Moose::Functor< Real > & _rho
Functor describing the density of the fluid.
static InputParameters validParams()
virtual Real computeElemRightHandSideContribution() override
std::unique_ptr< NumericVector< Number > > solution
void computeCellVelocity()
Update the cell values of the velocity variables.
virtual Real computeNeighborMatrixContribution() override
Real getVolumetricFaceFlux(const FaceInfo &fi) const
Get the volumetric face flux (used in advection terms)
void setCurrentFaceArea(const Real area)
virtual void print(std::ostream &os=libMesh::out) const
virtual Real computeBoundaryRHSContribution(const LinearFVBoundaryCondition &bc) override
const Point & normal() const
void paramError(const std::string ¶m, Args... args) const
unsigned int number() const
virtual const MooseVariableFieldBase & getVariable(const THREAD_ID tid, const std::string &var_name, Moose::VarKindType expected_var_type=Moose::VarKindType::VAR_ANY, Moose::VarFieldType expected_var_field_type=Moose::VarFieldType::VAR_FIELD_ANY) const=0
virtual void get_row(numeric_index_type i, std::vector< numeric_index_type > &indices, std::vector< T > &values) const override
std::vector< std::unique_ptr< NumericVector< Number > > > & selectPressureGradient(const bool updated_pressure)
Select the right pressure gradient field and return a reference to the container. ...
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2)=0
const Elem * elemPtr() const
const std::vector< std::vector< dof_id_type > > & dofIndices() const
virtual void get_diagonal(NumericVector< T > &dest) const override
FaceCenteredMapFunctor< Real, std::unordered_map< dof_id_type, Real > > & _face_mass_flux
A map functor from faces to mass fluxes which are used in the advection terms.
std::vector< std::unique_ptr< NumericVector< Number > > > _grad_p_current
for a PISO iteration we need to hold on to the original pressure gradient field.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
FEProblemBase & _fe_problem
SparseMatrix< Number > * matrix
Real getElemValue(const ElemInfo &elem_info, const StateArg &state) const
static const std::string pressure
IntRange< T > make_range(T beg, T end)
virtual MooseMesh & mesh() override
void mooseError(Args &&... args) const
virtual numeric_index_type row_start() const override
std::unique_ptr< NumericVector< Number > > current_local_solution
static InputParameters validParams()
const ConsoleStream _console
std::vector< const FaceInfo * > _flow_face_info
The subset of the FaceInfo objects that actually cover the subdomains which the flow field is defined...
bool hasBlocks(const SubdomainName &name) const
void populateCouplingFunctors(const std::vector< std::unique_ptr< NumericVector< Number >>> &raw_hbya, const std::vector< std::unique_ptr< NumericVector< Number >>> &raw_Ainv)
Populate the face values of the H/A and 1/A fields.
virtual void add(const numeric_index_type i, const T value)=0
void interpolate(InterpMethod m, T &result, const T2 &value1, const T3 &value2, const FaceInfo &fi, const bool one_is_elem)
LinearFVAnisotropicDiffusion * _p_diffusion_kernel
Pointer to the pressure diffusion term in the pressure Poisson equation.
auto index_range(const T &sizable)
void print(std::ostream &os=libMesh::out, const bool sparse=false) const
std::vector< std::unique_ptr< NumericVector< Number > > > _HbyA_raw
We hold on to the cell-based HbyA vectors so that we can easily reconstruct the cell velocities as we...
virtual void pointwise_divide(const NumericVector< T > &vec1, const NumericVector< T > &vec2)=0
void initCouplingField()
Initialize the coupling fields (HbyA and Ainv)
virtual System & system() override
void initFaceMassFlux()
Initialize the container for face velocities.
VarFaceNeighbors faceType(const std::pair< unsigned int, unsigned int > &var_sys) const
std::unique_ptr< NumericVector< Number > > _cell_volumes
We will hold a vector of cell volumes to make sure we can do volume corrections rapidly.
const unsigned int _dim
The dimension of the mesh, e.g. 3 for hexes and tets, 2 for quads and tris.