https://mooseframework.inl.gov
RhieChowMassFlux.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 // MOOSE includes
11 #include "RhieChowMassFlux.h"
12 #include "SubProblem.h"
13 #include "MooseMesh.h"
14 #include "NS.h"
15 #include "VectorCompositeFunctor.h"
16 #include "PIMPLE.h"
17 #include "SIMPLE.h"
18 #include "PetscVectorReader.h"
19 #include "LinearSystem.h"
21 
22 // libMesh includes
23 #include "libmesh/mesh_base.h"
24 #include "libmesh/elem_range.h"
25 #include "libmesh/petsc_matrix.h"
26 
27 using namespace libMesh;
28 
29 registerMooseObject("NavierStokesApp", RhieChowMassFlux);
30 
33 {
36 
37  params.addClassDescription("Computes H/A and 1/A together with face mass fluxes for segregated "
38  "momentum-pressure equations using linear systems.");
39 
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.");
46 
47  params.addRequiredParam<MooseFunctorName>(NS::density, "Density functor");
48 
49  // We disable the execution of this, should only provide functions
50  // for the SIMPLE executioner
51  ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
52  exec_enum.addAvailableFlags(EXEC_NONE);
53  exec_enum = {EXEC_NONE};
54  params.suppressParameter<ExecFlagEnum>("execute_on");
55 
56  // Pressure projection
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)");
61 
62  return params;
63 }
64 
66  : RhieChowFaceFluxProvider(params),
68  _moose_mesh(UserObject::_subproblem.mesh()),
69  _mesh(_moose_mesh.getMesh()),
70  _dim(blocksMaxDimension()),
71  _p(dynamic_cast<MooseLinearVariableFVReal *>(
72  &UserObject::_subproblem.getVariable(0, getParam<VariableName>(NS::pressure)))),
73  _vel(_dim, nullptr),
74  _HbyA_flux(_moose_mesh, blockIDs(), "HbyA_flux"),
75  _Ainv(_moose_mesh, blockIDs(), "Ainv"),
76  _face_mass_flux(
77  declareRestartableData<FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>>(
78  "face_flux", _moose_mesh, blockIDs(), "face_values")),
79  _rho(getFunctor<Real>(NS::density)),
80  _pressure_projection_method(getParam<MooseEnum>("pressure_projection_method"))
81 {
82  if (!_p)
83  paramError(NS::pressure, "the pressure must be a MooseLinearVariableFVReal.");
84  checkBlocks(*_p);
85 
86  std::vector<std::string> vel_names = {"u", "v", "w"};
87  for (const auto i : index_range(_vel))
88  {
89  _vel[i] = dynamic_cast<MooseLinearVariableFVReal *>(
90  &UserObject::_subproblem.getVariable(0, getParam<VariableName>(vel_names[i])));
91 
92  if (!_vel[i])
93  paramError(vel_names[i], "the velocity must be a MOOSELinearVariableFVReal.");
94  checkBlocks(*_vel[i]);
95  }
96 
97  // Register the elemental/face functors which will be queried in the pressure equation
98  for (const auto tid : make_range(libMesh::n_threads()))
99  {
102  }
103 
104  if (!dynamic_cast<SIMPLE *>(getMooseApp().getExecutioner()) &&
105  !dynamic_cast<PIMPLE *>(getMooseApp().getExecutioner()))
106  mooseError(this->name(),
107  " should only be used with a linear segregated thermal-hydraulics solver!");
108 }
109 
110 void
112  const std::vector<LinearSystem *> & momentum_systems,
113  const LinearSystem & pressure_system,
114  const std::vector<unsigned int> & momentum_system_numbers)
115 {
116  _momentum_systems = momentum_systems;
117  _momentum_system_numbers = momentum_system_numbers;
118  _pressure_system = &pressure_system;
120 
122  for (auto & system : _momentum_systems)
123  {
124  _global_momentum_system_numbers.push_back(system->number());
125  _momentum_implicit_systems.push_back(dynamic_cast<LinearImplicitSystem *>(&system->system()));
126  }
127 
129 }
130 
131 void
133 {
134  _HbyA_flux.clear();
135  _Ainv.clear();
136  _face_mass_flux.clear();
138 }
139 
140 void
142 {
143  // We fetch the pressure diffusion kernel to ensure that the face flux correction
144  // is consistent with the pressure discretization in the Poisson equation.
145  std::vector<LinearFVFluxKernel *> flux_kernel;
146  auto base_query = _fe_problem.theWarehouse()
147  .query()
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)
154  paramError(
155  "p_diffusion_kernel",
156  "The kernel with the given name could not be found or multiple instances were identified.");
157  _p_diffusion_kernel = dynamic_cast<LinearFVAnisotropicDiffusion *>(flux_kernel[0]);
158  if (!_p_diffusion_kernel)
159  paramError("p_diffusion_kernel",
160  "The provided diffusion kernel should of type LinearFVAnisotropicDiffusion!");
161 }
162 
163 void
165 {
166  // We cache the cell volumes into a petsc vector for corrections here so we can use
167  // the optimized petsc operations for the normalization
169  for (const auto & elem_info : _fe_problem.mesh().elemInfoVector())
170  // We have to check this because the variable might not be defined on the given
171  // block
172  if (hasBlocks(elem_info->subdomain_id()))
173  {
174  const auto elem_dof = elem_info->dofIndices()[_global_pressure_system_number][0];
175  _cell_volumes->set(elem_dof, elem_info->volume() * elem_info->coordFactor());
176  }
177 
178  _cell_volumes->close();
179 
180  _flow_face_info.clear();
181  for (auto & fi : _fe_problem.mesh().faceInfo())
182  if (hasBlocks(fi->elemPtr()->subdomain_id()) ||
183  (fi->neighborPtr() && hasBlocks(fi->neighborPtr()->subdomain_id())))
184  _flow_face_info.push_back(fi);
185 }
186 
187 void
189 {
190  for (const auto & pair : _HbyA_flux)
191  _HbyA_flux[pair.first] = 0;
192 
193  for (const auto & pair : _Ainv)
194  _Ainv[pair.first] = 0;
195 }
196 
197 void
199 {
200  using namespace Moose::FV;
201 
202  const auto time_arg = Moose::currentState();
203 
204  // We loop through the faces and compute the resulting face fluxes from the
205  // initial conditions for velocity
206  for (auto & fi : _flow_face_info)
207  {
208  RealVectorValue density_times_velocity;
209 
210  // On internal face we do a regular interpolation with geometric weights
211  if (_vel[0]->isInternalFace(*fi))
212  {
213  const auto & elem_info = *fi->elemInfo();
214  const auto & neighbor_info = *fi->neighborInfo();
215 
216  Real elem_rho = _rho(makeElemArg(fi->elemPtr()), time_arg);
217  Real neighbor_rho = _rho(makeElemArg(fi->neighborPtr()), time_arg);
218 
219  for (const auto dim_i : index_range(_vel))
220  interpolate(InterpMethod::Average,
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,
224  *fi,
225  true);
226  }
227  // On the boundary, we just take the boundary values
228  else
229  {
230  const bool elem_is_fluid = hasBlocks(fi->elemPtr()->subdomain_id());
231  const Elem * const boundary_elem = elem_is_fluid ? fi->elemPtr() : fi->neighborPtr();
232 
233  // We need this multiplier in case the face is an internal face and
234  const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
235  const Moose::FaceArg boundary_face{
236  fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr};
237 
238  const Real face_rho = _rho(boundary_face, time_arg);
239  for (const auto dim_i : index_range(_vel))
240  density_times_velocity(dim_i) = boundary_normal_multiplier * face_rho *
241  raw_value((*_vel[dim_i])(boundary_face, time_arg));
242  }
243 
244  _face_mass_flux[fi->id()] = density_times_velocity * fi->normal();
245  }
246 }
247 
248 Real
250 {
251  return _face_mass_flux.evaluate(&fi);
252 }
253 
254 Real
256 {
257  const Moose::FaceArg face_arg{&fi,
259  /*elem_is_upwind=*/true,
260  /*correct_skewness=*/false,
261  &fi.elem(),
262  /*state_limiter*/ nullptr};
263  const Real face_rho = _rho(face_arg, Moose::currentState());
264  return libmesh_map_find(_face_mass_flux, fi.id()) / face_rho;
265 }
266 
267 Real
269  const FaceInfo & fi,
270  const Moose::StateArg & time,
271  const THREAD_ID /*tid*/,
272  bool libmesh_dbg_var(subtract_mesh_velocity)) const
273 {
274  mooseAssert(!subtract_mesh_velocity, "RhieChowMassFlux does not support moving meshes yet!");
275 
277  mooseError("Interpolation methods other than Rhie-Chow are not supported!");
278  if (time.state != Moose::currentState().state)
279  mooseError("Older interpolation times are not supported!");
280 
281  return getVolumetricFaceFlux(fi);
282 }
283 
284 void
286 {
287  using namespace Moose::FV;
288 
289  const auto time_arg = Moose::currentState();
290 
291  // Petsc vector reader to make the repeated reading from the vector faster
293 
294  // We loop through the faces and compute the face fluxes using the pressure gradient
295  // and the momentum matrix/right hand side
296  for (auto & fi : _flow_face_info)
297  {
298  // Making sure the kernel knows which face we are on
300 
301  // We are setting this to 1.0 because we don't want to multiply the kernel contributions
302  // with the surface area yet. The surface area will be factored in in the advection kernels.
304 
305  Real p_grad_flux = 0.0;
306  if (_p->isInternalFace(*fi))
307  {
308  const auto & elem_info = *fi->elemInfo();
309  const auto & neighbor_info = *fi->neighborInfo();
310 
311  // Fetching the dof indices for the pressure variable
312  const auto elem_dof = elem_info.dofIndices()[_global_pressure_system_number][0];
313  const auto neighbor_dof = neighbor_info.dofIndices()[_global_pressure_system_number][0];
314 
315  // Fetching the values of the pressure for the element and the neighbor
316  const auto p_elem_value = p_reader(elem_dof);
317  const auto p_neighbor_value = p_reader(neighbor_dof);
318 
319  // Compute the elem matrix contributions for the face
320  const auto elem_matrix_contribution = _p_diffusion_kernel->computeElemMatrixContribution();
321  const auto neighbor_matrix_contribution =
323  const auto elem_rhs_contribution =
325 
326  // Compute the face flux from the matrix and right hand side contributions
327  p_grad_flux = (p_neighbor_value * neighbor_matrix_contribution +
328  p_elem_value * elem_matrix_contribution) -
329  elem_rhs_contribution;
330  }
331  else if (auto * bc_pointer = _p->getBoundaryCondition(*fi->boundaryIDs().begin()))
332  {
333  mooseAssert(fi->boundaryIDs().size() == 1, "We should only have one boundary on every face.");
334 
335  bc_pointer->setupFaceData(
336  fi, fi->faceType(std::make_pair(_p->number(), _global_pressure_system_number)));
337 
338  const ElemInfo & elem_info =
339  hasBlocks(fi->elemPtr()->subdomain_id()) ? *fi->elemInfo() : *fi->neighborInfo();
340  const auto p_elem_value = _p->getElemValue(elem_info, time_arg);
341  const auto matrix_contribution =
343  const auto rhs_contribution =
345 
346  // On the boundary, only the element side has a contribution
347  p_grad_flux = (p_elem_value * matrix_contribution - rhs_contribution);
348  }
349  // Compute the new face flux
350  _face_mass_flux[fi->id()] = -_HbyA_flux[fi->id()] + p_grad_flux;
351  }
352 }
353 
354 void
356 {
357  auto & pressure_gradient = _pressure_system->gradientContainer();
358 
359  // We set the dof value in the solution vector the same logic applies:
360  // u_C = -(H/A)_C - (1/A)_C*grad(p)_C where C is the cell index
361  for (const auto system_i : index_range(_momentum_implicit_systems))
362  {
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);
367  (*_momentum_implicit_systems[system_i]->solution) = *working_vector;
368  _momentum_implicit_systems[system_i]->update();
369  _momentum_systems[system_i]->setSolution(
370  *_momentum_implicit_systems[system_i]->current_local_solution);
371  }
372 }
373 
374 void
376 {
377  // We loop through the faces and populate the coupling fields (face H/A and 1/H)
378  // with 0s for now. Pressure corrector solves will always come after the
379  // momentum source so we expect these fields to change before the actual solve.
380  for (auto & fi : _fe_problem.mesh().faceInfo())
381  {
382  _Ainv[fi->id()];
383  _HbyA_flux[fi->id()];
384  }
385 }
386 
387 void
389  const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_hbya,
390  const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_Ainv)
391 {
392  // We have the raw H/A and 1/A vectors in a petsc format. This function
393  // will create face functors from them
394  using namespace Moose::FV;
395  const auto time_arg = Moose::currentState();
396 
397  // Create the petsc vector readers for faster repeated access
398  std::vector<PetscVectorReader> hbya_reader;
399  for (const auto dim_i : index_range(raw_hbya))
400  hbya_reader.emplace_back(*raw_hbya[dim_i]);
401 
402  std::vector<PetscVectorReader> ainv_reader;
403  for (const auto dim_i : index_range(raw_Ainv))
404  ainv_reader.emplace_back(*raw_Ainv[dim_i]);
405 
406  // We loop through the faces and populate the coupling fields (face H/A and 1/H)
407  for (auto & fi : _flow_face_info)
408  {
409  Real face_rho = 0;
410  RealVectorValue face_hbya;
411 
412  // We do the lookup in advance
413  auto & Ainv = _Ainv[fi->id()];
414 
415  // If it is internal, we just interpolate (using geometric weights) to the face
416  if (_vel[0]->isInternalFace(*fi))
417  {
418  // Get the dof indices for the element and the neighbor
419  const auto & elem_info = *fi->elemInfo();
420  const auto & neighbor_info = *fi->neighborInfo();
421  const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
422  const auto neighbor_dof = neighbor_info.dofIndices()[_global_momentum_system_numbers[0]][0];
423 
424  // Get the density values for the element and neighbor. We need this multiplication to make
425  // the coupling fields mass fluxes.
426  const Real elem_rho = _rho(makeElemArg(fi->elemPtr()), time_arg);
427  const Real neighbor_rho = _rho(makeElemArg(fi->neighborPtr()), time_arg);
428 
429  // Now we do the interpolation to the face
430  interpolate(Moose::FV::InterpMethod::Average, face_rho, elem_rho, neighbor_rho, *fi, true);
431  for (const auto dim_i : index_range(raw_hbya))
432  {
434  face_hbya(dim_i),
435  hbya_reader[dim_i](elem_dof),
436  hbya_reader[dim_i](neighbor_dof),
437  *fi,
438  true);
439  interpolate(InterpMethod::Average,
440  Ainv(dim_i),
441  elem_rho * ainv_reader[dim_i](elem_dof),
442  neighbor_rho * ainv_reader[dim_i](neighbor_dof),
443  *fi,
444  true);
445  }
446  }
447  else
448  {
449  const bool elem_is_fluid = hasBlocks(fi->elemPtr()->subdomain_id());
450 
451  // We need this multiplier in case the face is an internal face and
452  const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
453 
454  const ElemInfo & elem_info = elem_is_fluid ? *fi->elemInfo() : *fi->neighborInfo();
455  const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
456 
457  // If it is a Dirichlet BC, we use the dirichlet value the make sure the face flux
458  // is consistent
459  if (_vel[0]->isDirichletBoundaryFace(*fi))
460  {
461  const Moose::FaceArg boundary_face{
462  fi, Moose::FV::LimiterType::CentralDifference, true, false, elem_info.elem(), nullptr};
463  face_rho = _rho(boundary_face, Moose::currentState());
464 
465  for (const auto dim_i : make_range(_dim))
466  face_hbya(dim_i) =
467  -boundary_normal_multiplier *
468  MetaPhysicL::raw_value((*_vel[dim_i])(boundary_face, Moose::currentState()));
469  }
470  // Otherwise we just do a one-term expansion (so we just use the element value)
471  else
472  {
473  const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
474 
475  face_rho = _rho(makeElemArg(elem_info.elem()), time_arg);
476  for (const auto dim_i : make_range(_dim))
477  face_hbya(dim_i) = boundary_normal_multiplier * hbya_reader[dim_i](elem_dof);
478  }
479 
480  // We just do a one-term expansion for 1/A no matter what
481  const Real elem_rho = _rho(makeElemArg(elem_info.elem()), time_arg);
482  for (const auto dim_i : index_range(raw_Ainv))
483  Ainv(dim_i) = elem_rho * ainv_reader[dim_i](elem_dof);
484  }
485  // Lastly, we populate the face flux resulted by H/A
486  _HbyA_flux[fi->id()] = face_hbya * fi->normal() * face_rho;
487  }
488 }
489 
490 void
491 RhieChowMassFlux::computeHbyA(const bool with_updated_pressure, bool verbose)
492 {
493  if (verbose)
494  {
495  _console << "************************************" << std::endl;
496  _console << "Computing HbyA" << std::endl;
497  _console << "************************************" << std::endl;
498  }
500  "The momentum system shall be linked before calling this function!");
501 
502  auto & pressure_gradient = selectPressureGradient(with_updated_pressure);
503 
504  _HbyA_raw.clear();
505  _Ainv_raw.clear();
506  for (auto system_i : index_range(_momentum_systems))
507  {
508  LinearImplicitSystem * momentum_system = _momentum_implicit_systems[system_i];
509 
510  NumericVector<Number> & rhs = *(momentum_system->rhs);
511  NumericVector<Number> & current_local_solution = *(momentum_system->current_local_solution);
512  NumericVector<Number> & solution = *(momentum_system->solution);
513  PetscMatrix<Number> * mmat = dynamic_cast<PetscMatrix<Number> *>(momentum_system->matrix);
514  mooseAssert(mmat,
515  "The matrices used in the segregated INSFVRhieChow objects need to be convertable "
516  "to PetscMatrix!");
517 
518  if (verbose)
519  {
520  _console << "Matrix in rc object" << std::endl;
521  mmat->print();
522  }
523 
524  // First, we extract the diagonal and we will hold on to it for a little while
525  _Ainv_raw.push_back(current_local_solution.zero_clone());
526  NumericVector<Number> & Ainv = *(_Ainv_raw.back());
527 
528  mmat->get_diagonal(Ainv);
529 
530  if (verbose)
531  {
532  _console << "Velocity solution in H(u)" << std::endl;
533  solution.print();
534  }
535 
536  // Time to create H(u) = M_{offdiag} * u - b_{nonpressure}
537  _HbyA_raw.push_back(current_local_solution.zero_clone());
538  NumericVector<Number> & HbyA = *(_HbyA_raw.back());
539 
540  // We start with the matrix product part, we will do
541  // M*u - A*u for 2 reasons:
542  // 1, We assume A*u petsc operation is faster than setting the matrix diagonal to 0
543  // 2, In PISO loops we need to reuse the matrix so we can't just set the diagonals to 0
544 
545  // We create a working vector to ease some of the operations, we initialize its values
546  // with the current solution values to have something for the A*u term
547  auto working_vector = momentum_system->current_local_solution->zero_clone();
548  PetscVector<Number> * working_vector_petsc =
549  dynamic_cast<PetscVector<Number> *>(working_vector.get());
550  mooseAssert(working_vector_petsc,
551  "The vectors used in the RhieChowMassFlux objects need to be convertable "
552  "to PetscVectors!");
553 
554  mmat->vector_mult(HbyA, solution);
555  working_vector_petsc->pointwise_mult(Ainv, solution);
556  HbyA.add(-1.0, *working_vector_petsc);
557 
558  if (verbose)
559  {
560  _console << " H(u)" << std::endl;
561  HbyA.print();
562  }
563 
564  // We continue by adding the momentum right hand side contributions
565  HbyA.add(-1.0, rhs);
566 
567  // Unfortunately, the pressure forces are included in the momentum RHS
568  // so we have to correct them back
569  working_vector_petsc->pointwise_mult(*pressure_gradient[system_i], *_cell_volumes);
570  HbyA.add(-1.0, *working_vector_petsc);
571 
572  if (verbose)
573  {
574  _console << "total RHS" << std::endl;
575  rhs.print();
576  _console << "pressure RHS" << std::endl;
577  pressure_gradient[system_i]->print();
578  _console << " H(u)-rhs-relaxation_source" << std::endl;
579  HbyA.print();
580  }
581 
582  // It is time to create element-wise 1/A-s based on the the diagonal of the momentum matrix
583  *working_vector_petsc = 1.0;
584  Ainv.pointwise_divide(*working_vector_petsc, Ainv);
585 
586  // Create 1/A*(H(u)-RHS)
587  HbyA.pointwise_mult(HbyA, Ainv);
588 
589  if (verbose)
590  {
591  _console << " (H(u)-rhs)/A" << std::endl;
592  HbyA.print();
593  }
594 
595  if (_pressure_projection_method == "consistent")
596  {
597 
598  // Consistent Corrections to SIMPLE
599  // 1. Ainv_old = 1/a_p <- Ainv = 1/(a_p + \sum_n a_n)
600  // 2. H(u) <- H(u*) + H(u') = H(u*) - (Ainv - Ainv_old) * grad(p) * Vc
601 
602  if (verbose)
603  _console << "Performing SIMPLEC projection." << std::endl;
604 
605  // Lambda function to calculate the sum of diagonal and neighbor coefficients
606  auto get_row_sum = [mmat](NumericVector<Number> & sum_vector)
607  {
608  // Ensure the sum_vector is zeroed out
609  sum_vector.zero();
610 
611  // Local row size
612  const auto local_size = mmat->local_m();
613 
614  for (const auto row_i : make_range(local_size))
615  {
616  // Get all non-zero components of the row of the matrix
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);
621 
622  // Sum row elements (no absolute values)
623  const Real row_sum = std::accumulate(values.cbegin(), values.cend(), 0.0);
624 
625  // Add the sum of diagonal and elements to the sum_vector
626  sum_vector.add(global_index, row_sum);
627  }
628  sum_vector.close();
629  };
630 
631  // Create a temporary vector to store the sum of diagonal and neighbor coefficients
632  auto row_sum = current_local_solution.zero_clone();
633  get_row_sum(*row_sum);
634 
635  // Create vector with new inverse projection matrix
636  auto Ainv_full = current_local_solution.zero_clone();
637  *working_vector_petsc = 1.0;
638  Ainv_full->pointwise_divide(*working_vector_petsc, *row_sum);
639  const auto Ainv_full_old = Ainv_full->clone();
640 
641  // Correct HbyA
642  Ainv_full->add(-1.0, Ainv);
643  working_vector_petsc->pointwise_mult(*Ainv_full, *pressure_gradient[system_i]);
644  working_vector_petsc->pointwise_mult(*working_vector_petsc, *_cell_volumes);
645  HbyA.add(-1.0, *working_vector_petsc);
646 
647  // Correct Ainv
648  Ainv = *Ainv_full_old;
649  }
650 
651  Ainv.pointwise_mult(Ainv, *_cell_volumes);
652  }
653 
654  // We fill the 1/A and H/A functors
656 
657  if (verbose)
658  {
659  _console << "************************************" << std::endl;
660  _console << "DONE Computing HbyA " << std::endl;
661  _console << "************************************" << std::endl;
662  }
663 }
664 
665 std::vector<std::unique_ptr<NumericVector<Number>>> &
666 RhieChowMassFlux::selectPressureGradient(const bool updated_pressure)
667 {
668  if (updated_pressure)
669  {
670  _grad_p_current.clear();
671  for (const auto & component : _pressure_system->gradientContainer())
672  _grad_p_current.push_back(component->clone());
673  }
674 
675  return _grad_p_current;
676 }
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...
unsigned int n_threads()
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
Definition: SCM.h:35
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
Definition: NS.h:153
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)
MeshBase & mesh
static const std::string density
Definition: NS.h:33
auto raw_value(const Eigen::Map< T > &in)
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()
SubProblem & _subproblem
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 &params)
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 &param, 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
Query query()
Real getElemValue(const ElemInfo &elem_info, const StateArg &state) const
static const std::string pressure
Definition: NS.h:56
const THREAD_ID _tid
IntRange< T > make_range(T beg, T end)
virtual MooseMesh & mesh() override
dof_id_type id() const
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)
StateArg currentState()
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)
unsigned int state
virtual System & system() override
unsigned int THREAD_ID
uint8_t dof_id_type
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.