Line data Source code
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"
20 : #include "LinearFVBoundaryCondition.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 :
31 : InputParameters
32 1784 : RhieChowMassFlux::validParams()
33 : {
34 1784 : auto params = RhieChowFaceFluxProvider::validParams();
35 1784 : params += NonADFunctorInterface::validParams();
36 :
37 1784 : params.addClassDescription("Computes H/A and 1/A together with face mass fluxes for segregated "
38 : "momentum-pressure equations using linear systems.");
39 :
40 1784 : params.addRequiredParam<VariableName>(NS::pressure, "The pressure variable.");
41 3568 : params.addRequiredParam<VariableName>("u", "The x-component of velocity");
42 3568 : params.addParam<VariableName>("v", "The y-component of velocity");
43 3568 : params.addParam<VariableName>("w", "The z-component of velocity");
44 3568 : params.addRequiredParam<std::string>("p_diffusion_kernel",
45 : "The diffusion kernel acting on the pressure.");
46 3568 : params.addParam<std::vector<std::vector<std::string>>>(
47 : "body_force_kernel_names",
48 : {},
49 : "The body force kernel names."
50 : "this double vector would have size index_x_dim: 'f1x f2x; f1y f2y; f1z f2z'");
51 :
52 1784 : params.addRequiredParam<MooseFunctorName>(NS::density, "Density functor");
53 :
54 : // We disable the execution of this, should only provide functions
55 : // for the SIMPLE executioner
56 1784 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
57 1784 : exec_enum.addAvailableFlags(EXEC_NONE);
58 5352 : exec_enum = {EXEC_NONE};
59 1784 : params.suppressParameter<ExecFlagEnum>("execute_on");
60 :
61 : // Pressure projection
62 3568 : params.addParam<MooseEnum>("pressure_projection_method",
63 5352 : MooseEnum("standard consistent", "standard"),
64 : "The method to use in the pressure projection for Ainv - "
65 : "standard (SIMPLE) or consistent (SIMPLEC)");
66 1784 : return params;
67 1784 : }
68 :
69 892 : RhieChowMassFlux::RhieChowMassFlux(const InputParameters & params)
70 : : RhieChowFaceFluxProvider(params),
71 : NonADFunctorInterface(this),
72 1784 : _moose_mesh(UserObject::_subproblem.mesh()),
73 892 : _mesh(_moose_mesh.getMesh()),
74 892 : _dim(blocksMaxDimension()),
75 892 : _p(dynamic_cast<MooseLinearVariableFVReal *>(
76 892 : &UserObject::_subproblem.getVariable(0, getParam<VariableName>(NS::pressure)))),
77 892 : _vel(_dim, nullptr),
78 1784 : _HbyA_flux(_moose_mesh, blockIDs(), "HbyA_flux"),
79 1784 : _Ainv(_moose_mesh, blockIDs(), "Ainv"),
80 892 : _face_mass_flux(
81 1784 : declareRestartableData<FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>>(
82 : "face_flux", _moose_mesh, blockIDs(), "face_values")),
83 2676 : _body_force_kernel_names(
84 : getParam<std::vector<std::vector<std::string>>>("body_force_kernel_names")),
85 892 : _rho(getFunctor<Real>(NS::density)),
86 3568 : _pressure_projection_method(getParam<MooseEnum>("pressure_projection_method"))
87 : {
88 892 : if (!_p)
89 0 : paramError(NS::pressure, "the pressure must be a MooseLinearVariableFVReal.");
90 892 : checkBlocks(*_p);
91 :
92 892 : std::vector<std::string> vel_names = {"u", "v", "w"};
93 2637 : for (const auto i : index_range(_vel))
94 : {
95 1745 : _vel[i] = dynamic_cast<MooseLinearVariableFVReal *>(
96 1745 : &UserObject::_subproblem.getVariable(0, getParam<VariableName>(vel_names[i])));
97 :
98 1745 : if (!_vel[i])
99 0 : paramError(vel_names[i], "the velocity must be a MOOSELinearVariableFVReal.");
100 1745 : checkBlocks(*_vel[i]);
101 : }
102 :
103 : // Register the elemental/face functors which will be queried in the pressure equation
104 1784 : for (const auto tid : make_range(libMesh::n_threads()))
105 : {
106 892 : UserObject::_subproblem.addFunctor("Ainv", _Ainv, tid);
107 1784 : UserObject::_subproblem.addFunctor("HbyA", _HbyA_flux, tid);
108 : }
109 :
110 892 : if (!dynamic_cast<SIMPLE *>(getMooseApp().getExecutioner()) &&
111 143 : !dynamic_cast<PIMPLE *>(getMooseApp().getExecutioner()))
112 0 : mooseError(this->name(),
113 : " should only be used with a linear segregated thermal-hydraulics solver!");
114 892 : }
115 :
116 : void
117 890 : RhieChowMassFlux::linkMomentumPressureSystems(
118 : const std::vector<LinearSystem *> & momentum_systems,
119 : const LinearSystem & pressure_system,
120 : const std::vector<unsigned int> & momentum_system_numbers)
121 : {
122 890 : _momentum_systems = momentum_systems;
123 890 : _momentum_system_numbers = momentum_system_numbers;
124 890 : _pressure_system = &pressure_system;
125 890 : _global_pressure_system_number = _pressure_system->number();
126 :
127 890 : _momentum_implicit_systems.clear();
128 2631 : for (auto & system : _momentum_systems)
129 : {
130 1741 : _global_momentum_system_numbers.push_back(system->number());
131 1741 : _momentum_implicit_systems.push_back(dynamic_cast<LinearImplicitSystem *>(&system->system()));
132 : }
133 :
134 890 : setupMeshInformation();
135 890 : }
136 :
137 : void
138 9 : RhieChowMassFlux::meshChanged()
139 : {
140 : _HbyA_flux.clear();
141 : _Ainv.clear();
142 9 : _face_mass_flux.clear();
143 9 : setupMeshInformation();
144 9 : }
145 :
146 : void
147 890 : RhieChowMassFlux::initialSetup()
148 : {
149 : // We fetch the pressure diffusion kernel to ensure that the face flux correction
150 : // is consistent with the pressure discretization in the Poisson equation.
151 : std::vector<LinearFVFluxKernel *> flux_kernel;
152 890 : auto base_query = _fe_problem.theWarehouse()
153 890 : .query()
154 890 : .template condition<AttribThread>(_tid)
155 890 : .template condition<AttribSysNum>(_p->sys().number())
156 890 : .template condition<AttribSystem>("LinearFVFluxKernel")
157 2670 : .template condition<AttribName>(getParam<std::string>("p_diffusion_kernel"))
158 890 : .queryInto(flux_kernel);
159 890 : if (flux_kernel.size() != 1)
160 0 : paramError(
161 : "p_diffusion_kernel",
162 : "The kernel with the given name could not be found or multiple instances were identified.");
163 890 : _p_diffusion_kernel = dynamic_cast<LinearFVAnisotropicDiffusion *>(flux_kernel[0]);
164 890 : if (!_p_diffusion_kernel)
165 0 : paramError("p_diffusion_kernel",
166 : "The provided diffusion kernel should of type LinearFVAnisotropicDiffusion!");
167 :
168 : // We fetch the body forces kernel to ensure that the face flux correction
169 : // is accurate.
170 :
171 : // Check if components match the dimension.
172 :
173 890 : if (!_body_force_kernel_names.empty())
174 : {
175 57 : if (_body_force_kernel_names.size() != _dim)
176 0 : paramError("body_force_kernel_names",
177 : "The dimension of the body force vector does not match the problem dimension.");
178 :
179 57 : _body_force_kernels.resize(_dim);
180 :
181 171 : for (const auto dim_i : make_range(_dim))
182 228 : for (const auto & force_name : _body_force_kernel_names[dim_i])
183 : {
184 : std::vector<LinearFVElementalKernel *> temp_storage;
185 114 : auto base_query_force = _fe_problem.theWarehouse()
186 114 : .query()
187 114 : .template condition<AttribThread>(_tid)
188 228 : .template condition<AttribSysNum>(_vel[dim_i]->sys().number())
189 114 : .template condition<AttribSystem>("LinearFVElementalKernel")
190 114 : .template condition<AttribName>(force_name)
191 114 : .queryInto(temp_storage);
192 114 : if (temp_storage.size() != 1)
193 0 : paramError("body_force_kernel_names",
194 0 : "The kernel with the given name: " + force_name +
195 : " could not be found or multiple instances were identified.");
196 114 : _body_force_kernels[dim_i].push_back(temp_storage[0]);
197 114 : }
198 : }
199 890 : }
200 :
201 : void
202 899 : RhieChowMassFlux::setupMeshInformation()
203 : {
204 : // We cache the cell volumes into a petsc vector for corrections here so we can use
205 : // the optimized petsc operations for the normalization
206 1798 : _cell_volumes = _pressure_system->currentSolution()->zero_clone();
207 127866 : for (const auto & elem_info : _fe_problem.mesh().elemInfoVector())
208 : // We have to check this because the variable might not be defined on the given
209 : // block
210 126967 : if (hasBlocks(elem_info->subdomain_id()))
211 : {
212 121999 : const auto elem_dof = elem_info->dofIndices()[_global_pressure_system_number][0];
213 121999 : _cell_volumes->set(elem_dof, elem_info->volume() * elem_info->coordFactor());
214 : }
215 :
216 899 : _cell_volumes->close();
217 :
218 899 : _flow_face_info.clear();
219 251598 : for (auto & fi : _fe_problem.mesh().faceInfo())
220 259832 : if (hasBlocks(fi->elemPtr()->subdomain_id()) ||
221 19466 : (fi->neighborPtr() && hasBlocks(fi->neighborPtr()->subdomain_id())))
222 240559 : _flow_face_info.push_back(fi);
223 899 : }
224 :
225 : void
226 0 : RhieChowMassFlux::initialize()
227 : {
228 0 : for (const auto & pair : _HbyA_flux)
229 0 : _HbyA_flux[pair.first] = 0;
230 :
231 0 : for (const auto & pair : _Ainv)
232 0 : _Ainv[pair.first] = 0;
233 0 : }
234 :
235 : void
236 812 : RhieChowMassFlux::initFaceMassFlux()
237 : {
238 : using namespace Moose::FV;
239 :
240 812 : const auto time_arg = Moose::currentState();
241 :
242 : // We loop through the faces and compute the resulting face fluxes from the
243 : // initial conditions for velocity
244 212750 : for (auto & fi : _flow_face_info)
245 : {
246 : RealVectorValue density_times_velocity;
247 :
248 : // On internal face we do a regular interpolation with geometric weights
249 211938 : if (_vel[0]->isInternalFace(*fi))
250 : {
251 189988 : const auto & elem_info = *fi->elemInfo();
252 : const auto & neighbor_info = *fi->neighborInfo();
253 :
254 189988 : Real elem_rho = _rho(makeElemArg(fi->elemPtr()), time_arg);
255 379976 : Real neighbor_rho = _rho(makeElemArg(fi->neighborPtr()), time_arg);
256 :
257 560870 : for (const auto dim_i : index_range(_vel))
258 370882 : interpolate(InterpMethod::Average,
259 : density_times_velocity(dim_i),
260 370882 : _vel[dim_i]->getElemValue(elem_info, time_arg) * elem_rho,
261 741764 : _vel[dim_i]->getElemValue(neighbor_info, time_arg) * neighbor_rho,
262 : *fi,
263 : true);
264 : }
265 : // On the boundary, we just take the boundary values
266 : else
267 : {
268 21950 : const bool elem_is_fluid = hasBlocks(fi->elemPtr()->subdomain_id());
269 21950 : const Elem * const boundary_elem = elem_is_fluid ? fi->elemPtr() : fi->neighborPtr();
270 :
271 : // We need this multiplier in case the face is an internal face and
272 21950 : const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
273 21950 : const Moose::FaceArg boundary_face{
274 21950 : fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr};
275 :
276 21950 : const Real face_rho = _rho(boundary_face, time_arg);
277 67260 : for (const auto dim_i : index_range(_vel))
278 45310 : density_times_velocity(dim_i) = boundary_normal_multiplier * face_rho *
279 45310 : raw_value((*_vel[dim_i])(boundary_face, time_arg));
280 : }
281 :
282 211938 : _face_mass_flux[fi->id()] = density_times_velocity * fi->normal();
283 : }
284 812 : }
285 :
286 : Real
287 80663466 : RhieChowMassFlux::getMassFlux(const FaceInfo & fi) const
288 : {
289 80663466 : return _face_mass_flux.evaluate(&fi);
290 : }
291 :
292 : Real
293 3813098 : RhieChowMassFlux::getVolumetricFaceFlux(const FaceInfo & fi) const
294 : {
295 3813098 : const Moose::FaceArg face_arg{&fi,
296 : /*limiter_type=*/Moose::FV::LimiterType::CentralDifference,
297 : /*elem_is_upwind=*/true,
298 : /*correct_skewness=*/false,
299 : &fi.elem(),
300 3813098 : /*state_limiter*/ nullptr};
301 3813098 : const Real face_rho = _rho(face_arg, Moose::currentState());
302 3813098 : return libmesh_map_find(_face_mass_flux, fi.id()) / face_rho;
303 : }
304 :
305 : Real
306 3652 : RhieChowMassFlux::getVolumetricFaceFlux(const Moose::FV::InterpMethod m,
307 : const FaceInfo & fi,
308 : const Moose::StateArg & time,
309 : const THREAD_ID /*tid*/,
310 : bool libmesh_dbg_var(subtract_mesh_velocity)) const
311 : {
312 : mooseAssert(!subtract_mesh_velocity, "RhieChowMassFlux does not support moving meshes yet!");
313 :
314 3652 : if (m != Moose::FV::InterpMethod::RhieChow)
315 0 : mooseError("Interpolation methods other than Rhie-Chow are not supported!");
316 3652 : if (time.state != Moose::currentState().state)
317 0 : mooseError("Older interpolation times are not supported!");
318 :
319 3652 : return getVolumetricFaceFlux(fi);
320 : }
321 :
322 : void
323 171133 : RhieChowMassFlux::computeFaceMassFlux()
324 : {
325 : using namespace Moose::FV;
326 :
327 171133 : const auto time_arg = Moose::currentState();
328 :
329 : // Petsc vector reader to make the repeated reading from the vector faster
330 171133 : PetscVectorReader p_reader(*_pressure_system->system().current_local_solution);
331 :
332 : // We loop through the faces and compute the face fluxes using the pressure gradient
333 : // and the momentum matrix/right hand side
334 30456876 : for (auto & fi : _flow_face_info)
335 : {
336 : // Making sure the kernel knows which face we are on
337 30285743 : _p_diffusion_kernel->setupFaceData(fi);
338 :
339 : // We are setting this to 1.0 because we don't want to multiply the kernel contributions
340 : // with the surface area yet. The surface area will be factored in in the advection kernels.
341 30285743 : _p_diffusion_kernel->setCurrentFaceArea(1.0);
342 :
343 : Real p_grad_flux = 0.0;
344 30285743 : if (_p->isInternalFace(*fi))
345 : {
346 27095581 : const auto & elem_info = *fi->elemInfo();
347 : const auto & neighbor_info = *fi->neighborInfo();
348 :
349 : // Fetching the dof indices for the pressure variable
350 27095581 : const auto elem_dof = elem_info.dofIndices()[_global_pressure_system_number][0];
351 27095581 : const auto neighbor_dof = neighbor_info.dofIndices()[_global_pressure_system_number][0];
352 :
353 : // Fetching the values of the pressure for the element and the neighbor
354 27095581 : const auto p_elem_value = p_reader(elem_dof);
355 27095581 : const auto p_neighbor_value = p_reader(neighbor_dof);
356 :
357 : // Compute the elem matrix contributions for the face
358 27095581 : const auto elem_matrix_contribution = _p_diffusion_kernel->computeElemMatrixContribution();
359 : const auto neighbor_matrix_contribution =
360 27095581 : _p_diffusion_kernel->computeNeighborMatrixContribution();
361 : const auto elem_rhs_contribution =
362 27095581 : _p_diffusion_kernel->computeElemRightHandSideContribution();
363 :
364 : // Compute the face flux from the matrix and right hand side contributions
365 27095581 : p_grad_flux = (p_neighbor_value * neighbor_matrix_contribution +
366 27095581 : p_elem_value * elem_matrix_contribution) -
367 : elem_rhs_contribution;
368 : }
369 3190162 : else if (auto * bc_pointer = _p->getBoundaryCondition(*fi->boundaryIDs().begin()))
370 : {
371 : mooseAssert(fi->boundaryIDs().size() == 1, "We should only have one boundary on every face.");
372 :
373 2135848 : bc_pointer->setupFaceData(
374 2135848 : fi, fi->faceType(std::make_pair(_p->number(), _global_pressure_system_number)));
375 :
376 : const ElemInfo & elem_info =
377 2135848 : hasBlocks(fi->elemPtr()->subdomain_id()) ? *fi->elemInfo() : *fi->neighborInfo();
378 2135848 : const auto p_elem_value = _p->getElemValue(elem_info, time_arg);
379 : const auto matrix_contribution =
380 2135848 : _p_diffusion_kernel->computeBoundaryMatrixContribution(*bc_pointer);
381 : const auto rhs_contribution =
382 2135848 : _p_diffusion_kernel->computeBoundaryRHSContribution(*bc_pointer);
383 :
384 : // On the boundary, only the element side has a contribution
385 2135848 : p_grad_flux = (p_elem_value * matrix_contribution - rhs_contribution);
386 : }
387 : // Compute the new face flux
388 30285743 : _face_mass_flux[fi->id()] = -_HbyA_flux[fi->id()] + p_grad_flux;
389 : }
390 171133 : }
391 :
392 : void
393 192833 : RhieChowMassFlux::computeCellVelocity()
394 : {
395 192833 : auto & pressure_gradient = _pressure_system->gradientContainer();
396 :
397 : // We set the dof value in the solution vector the same logic applies:
398 : // u_C = -(H/A)_C - (1/A)_C*grad(p)_C where C is the cell index
399 515686 : for (const auto system_i : index_range(_momentum_implicit_systems))
400 : {
401 322853 : auto working_vector = _Ainv_raw[system_i]->clone();
402 322853 : working_vector->pointwise_mult(*working_vector, *pressure_gradient[system_i]);
403 322853 : working_vector->add(*_HbyA_raw[system_i]);
404 322853 : working_vector->scale(-1.0);
405 322853 : (*_momentum_implicit_systems[system_i]->solution) = *working_vector;
406 322853 : _momentum_implicit_systems[system_i]->update();
407 322853 : _momentum_systems[system_i]->setSolution(
408 322853 : *_momentum_implicit_systems[system_i]->current_local_solution);
409 322853 : }
410 192833 : }
411 :
412 : void
413 890 : RhieChowMassFlux::initCouplingField()
414 : {
415 : // We loop through the faces and populate the coupling fields (face H/A and 1/H)
416 : // with 0s for now. Pressure corrector solves will always come after the
417 : // momentum source so we expect these fields to change before the actual solve.
418 250797 : for (auto & fi : _fe_problem.mesh().faceInfo())
419 : {
420 249907 : _Ainv[fi->id()];
421 249907 : _HbyA_flux[fi->id()];
422 : }
423 890 : }
424 :
425 : void
426 192833 : RhieChowMassFlux::populateCouplingFunctors(
427 : const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_hbya,
428 : const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_Ainv)
429 : {
430 : // We have the raw H/A and 1/A vectors in a petsc format. This function
431 : // will create face functors from them
432 : using namespace Moose::FV;
433 192833 : const auto time_arg = Moose::currentState();
434 :
435 : // Create the petsc vector readers for faster repeated access
436 : std::vector<PetscVectorReader> hbya_reader;
437 515686 : for (const auto dim_i : index_range(raw_hbya))
438 322853 : hbya_reader.emplace_back(*raw_hbya[dim_i]);
439 :
440 : std::vector<PetscVectorReader> ainv_reader;
441 515686 : for (const auto dim_i : index_range(raw_Ainv))
442 322853 : ainv_reader.emplace_back(*raw_Ainv[dim_i]);
443 :
444 : // We loop through the faces and populate the coupling fields (face H/A and 1/H)
445 32895076 : for (auto & fi : _flow_face_info)
446 : {
447 : Real face_rho = 0;
448 : RealVectorValue face_hbya;
449 :
450 : // We do the lookup in advance
451 32702243 : auto & Ainv = _Ainv[fi->id()];
452 :
453 : // If it is internal, we just interpolate (using geometric weights) to the face
454 32702243 : if (_vel[0]->isInternalFace(*fi))
455 : {
456 : // Get the dof indices for the element and the neighbor
457 29010081 : const auto & elem_info = *fi->elemInfo();
458 : const auto & neighbor_info = *fi->neighborInfo();
459 29010081 : const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
460 29010081 : const auto neighbor_dof = neighbor_info.dofIndices()[_global_momentum_system_numbers[0]][0];
461 :
462 : // Get the density values for the element and neighbor. We need this multiplication to make
463 : // the coupling fields mass fluxes.
464 29010081 : const Real elem_rho = _rho(makeElemArg(fi->elemPtr()), time_arg);
465 58020162 : const Real neighbor_rho = _rho(makeElemArg(fi->neighborPtr()), time_arg);
466 :
467 : // Now we do the interpolation to the face
468 29010081 : interpolate(Moose::FV::InterpMethod::Average, face_rho, elem_rho, neighbor_rho, *fi, true);
469 85613768 : for (const auto dim_i : index_range(raw_hbya))
470 : {
471 56603687 : interpolate(Moose::FV::InterpMethod::Average,
472 : face_hbya(dim_i),
473 56603687 : hbya_reader[dim_i](elem_dof),
474 56603687 : hbya_reader[dim_i](neighbor_dof),
475 : *fi,
476 : true);
477 56603687 : interpolate(InterpMethod::Average,
478 : Ainv(dim_i),
479 56603687 : elem_rho * ainv_reader[dim_i](elem_dof),
480 113207374 : neighbor_rho * ainv_reader[dim_i](neighbor_dof),
481 : *fi,
482 : true);
483 : }
484 : }
485 : else
486 : {
487 3692162 : const bool elem_is_fluid = hasBlocks(fi->elemPtr()->subdomain_id());
488 :
489 : // We need this multiplier in case the face is an internal face and
490 3692162 : const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
491 :
492 3692162 : const ElemInfo & elem_info = elem_is_fluid ? *fi->elemInfo() : *fi->neighborInfo();
493 3692162 : const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
494 :
495 : // If it is a Dirichlet BC, we use the dirichlet value the make sure the face flux
496 : // is consistent
497 3692162 : if (_vel[0]->isDirichletBoundaryFace(*fi))
498 : {
499 3073787 : const Moose::FaceArg boundary_face{
500 3073787 : fi, Moose::FV::LimiterType::CentralDifference, true, false, elem_info.elem(), nullptr};
501 3073787 : face_rho = _rho(boundary_face, Moose::currentState());
502 :
503 9225200 : for (const auto dim_i : make_range(_dim))
504 : {
505 :
506 6151413 : face_hbya(dim_i) =
507 6151413 : -MetaPhysicL::raw_value((*_vel[dim_i])(boundary_face, Moose::currentState()));
508 :
509 6151413 : if (!_body_force_kernel_names.empty())
510 725760 : for (const auto & force_kernel : _body_force_kernels[dim_i])
511 : {
512 362880 : force_kernel->setCurrentElemInfo(&elem_info);
513 725760 : face_hbya(dim_i) -= force_kernel->computeRightHandSideContribution() *
514 362880 : ainv_reader[dim_i](elem_dof) /
515 : elem_info.volume(); // zero-term expansion
516 : }
517 6151413 : face_hbya(dim_i) *= boundary_normal_multiplier;
518 : }
519 : }
520 : // Otherwise we just do a one-term expansion (so we just use the element value)
521 : else
522 : {
523 618375 : const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
524 :
525 618375 : face_rho = _rho(makeElemArg(elem_info.elem()), time_arg);
526 1820372 : for (const auto dim_i : make_range(_dim))
527 1201997 : face_hbya(dim_i) = boundary_normal_multiplier * hbya_reader[dim_i](elem_dof);
528 : }
529 :
530 : // We just do a one-term expansion for 1/A no matter what
531 3692162 : const Real elem_rho = _rho(makeElemArg(elem_info.elem()), time_arg);
532 11045572 : for (const auto dim_i : index_range(raw_Ainv))
533 7353410 : Ainv(dim_i) = elem_rho * ainv_reader[dim_i](elem_dof);
534 : }
535 : // Lastly, we populate the face flux resulted by H/A
536 32702243 : _HbyA_flux[fi->id()] = face_hbya * fi->normal() * face_rho;
537 : }
538 192833 : }
539 :
540 : void
541 192833 : RhieChowMassFlux::computeHbyA(const bool with_updated_pressure, bool verbose)
542 : {
543 192833 : if (verbose)
544 : {
545 0 : _console << "************************************" << std::endl;
546 0 : _console << "Computing HbyA" << std::endl;
547 0 : _console << "************************************" << std::endl;
548 : }
549 : mooseAssert(_momentum_implicit_systems.size() && _momentum_implicit_systems[0],
550 : "The momentum system shall be linked before calling this function!");
551 :
552 192833 : auto & pressure_gradient = selectPressureGradient(with_updated_pressure);
553 :
554 192833 : _HbyA_raw.clear();
555 192833 : _Ainv_raw.clear();
556 515686 : for (auto system_i : index_range(_momentum_systems))
557 : {
558 322853 : LinearImplicitSystem * momentum_system = _momentum_implicit_systems[system_i];
559 :
560 322853 : NumericVector<Number> & rhs = *(momentum_system->rhs);
561 : NumericVector<Number> & current_local_solution = *(momentum_system->current_local_solution);
562 : NumericVector<Number> & solution = *(momentum_system->solution);
563 322853 : PetscMatrix<Number> * mmat = dynamic_cast<PetscMatrix<Number> *>(momentum_system->matrix);
564 : mooseAssert(mmat,
565 : "The matrices used in the segregated INSFVRhieChow objects need to be convertable "
566 : "to PetscMatrix!");
567 :
568 322853 : if (verbose)
569 : {
570 0 : _console << "Matrix in rc object" << std::endl;
571 0 : mmat->print();
572 : }
573 :
574 : // First, we extract the diagonal and we will hold on to it for a little while
575 645706 : _Ainv_raw.push_back(current_local_solution.zero_clone());
576 : NumericVector<Number> & Ainv = *(_Ainv_raw.back());
577 :
578 322853 : mmat->get_diagonal(Ainv);
579 :
580 322853 : if (verbose)
581 : {
582 0 : _console << "Velocity solution in H(u)" << std::endl;
583 0 : solution.print();
584 : }
585 :
586 : // Time to create H(u) = M_{offdiag} * u - b_{nonpressure}
587 645706 : _HbyA_raw.push_back(current_local_solution.zero_clone());
588 : NumericVector<Number> & HbyA = *(_HbyA_raw.back());
589 :
590 : // We start with the matrix product part, we will do
591 : // M*u - A*u for 2 reasons:
592 : // 1, We assume A*u petsc operation is faster than setting the matrix diagonal to 0
593 : // 2, In PISO loops we need to reuse the matrix so we can't just set the diagonals to 0
594 :
595 : // We create a working vector to ease some of the operations, we initialize its values
596 : // with the current solution values to have something for the A*u term
597 322853 : auto working_vector = momentum_system->current_local_solution->zero_clone();
598 : PetscVector<Number> * working_vector_petsc =
599 322853 : dynamic_cast<PetscVector<Number> *>(working_vector.get());
600 : mooseAssert(working_vector_petsc,
601 : "The vectors used in the RhieChowMassFlux objects need to be convertable "
602 : "to PetscVectors!");
603 :
604 322853 : mmat->vector_mult(HbyA, solution);
605 322853 : working_vector_petsc->pointwise_mult(Ainv, solution);
606 322853 : HbyA.add(-1.0, *working_vector_petsc);
607 :
608 322853 : if (verbose)
609 : {
610 0 : _console << " H(u)" << std::endl;
611 0 : HbyA.print();
612 : }
613 :
614 : // We continue by adding the momentum right hand side contributions
615 322853 : HbyA.add(-1.0, rhs);
616 :
617 : // Unfortunately, the pressure forces are included in the momentum RHS
618 : // so we have to correct them back
619 322853 : working_vector_petsc->pointwise_mult(*pressure_gradient[system_i], *_cell_volumes);
620 322853 : HbyA.add(-1.0, *working_vector_petsc);
621 :
622 322853 : if (verbose)
623 : {
624 0 : _console << "total RHS" << std::endl;
625 0 : rhs.print();
626 0 : _console << "pressure RHS" << std::endl;
627 0 : pressure_gradient[system_i]->print();
628 0 : _console << " H(u)-rhs-relaxation_source" << std::endl;
629 0 : HbyA.print();
630 : }
631 :
632 : // It is time to create element-wise 1/A-s based on the the diagonal of the momentum matrix
633 322853 : *working_vector_petsc = 1.0;
634 322853 : Ainv.pointwise_divide(*working_vector_petsc, Ainv);
635 :
636 : // Create 1/A*(H(u)-RHS)
637 322853 : HbyA.pointwise_mult(HbyA, Ainv);
638 :
639 322853 : if (verbose)
640 : {
641 0 : _console << " (H(u)-rhs)/A" << std::endl;
642 0 : HbyA.print();
643 : }
644 :
645 322853 : if (_pressure_projection_method == "consistent")
646 : {
647 :
648 : // Consistent Corrections to SIMPLE
649 : // 1. Ainv_old = 1/a_p <- Ainv = 1/(a_p + \sum_n a_n)
650 : // 2. H(u) <- H(u*) + H(u') = H(u*) - (Ainv - Ainv_old) * grad(p) * Vc
651 :
652 1053 : if (verbose)
653 0 : _console << "Performing SIMPLEC projection." << std::endl;
654 :
655 : // Lambda function to calculate the sum of diagonal and neighbor coefficients
656 1053 : auto get_row_sum = [mmat](NumericVector<Number> & sum_vector)
657 : {
658 : // Ensure the sum_vector is zeroed out
659 1053 : sum_vector.zero();
660 :
661 : // Local row size
662 1053 : const auto local_size = mmat->local_m();
663 :
664 55593 : for (const auto row_i : make_range(local_size))
665 : {
666 : // Get all non-zero components of the row of the matrix
667 54540 : const auto global_index = mmat->row_start() + row_i;
668 : std::vector<numeric_index_type> indices;
669 : std::vector<Real> values;
670 54540 : mmat->get_row(global_index, indices, values);
671 :
672 : // Sum row elements (no absolute values)
673 : const Real row_sum = std::accumulate(values.cbegin(), values.cend(), 0.0);
674 :
675 : // Add the sum of diagonal and elements to the sum_vector
676 54540 : sum_vector.add(global_index, row_sum);
677 54540 : }
678 1053 : sum_vector.close();
679 1053 : };
680 :
681 : // Create a temporary vector to store the sum of diagonal and neighbor coefficients
682 1053 : auto row_sum = current_local_solution.zero_clone();
683 1053 : get_row_sum(*row_sum);
684 :
685 : // Create vector with new inverse projection matrix
686 1053 : auto Ainv_full = current_local_solution.zero_clone();
687 1053 : *working_vector_petsc = 1.0;
688 1053 : Ainv_full->pointwise_divide(*working_vector_petsc, *row_sum);
689 1053 : const auto Ainv_full_old = Ainv_full->clone();
690 :
691 : // Correct HbyA
692 1053 : Ainv_full->add(-1.0, Ainv);
693 1053 : working_vector_petsc->pointwise_mult(*Ainv_full, *pressure_gradient[system_i]);
694 1053 : working_vector_petsc->pointwise_mult(*working_vector_petsc, *_cell_volumes);
695 1053 : HbyA.add(-1.0, *working_vector_petsc);
696 :
697 : // Correct Ainv
698 1053 : Ainv = *Ainv_full_old;
699 1053 : }
700 :
701 322853 : Ainv.pointwise_mult(Ainv, *_cell_volumes);
702 322853 : }
703 :
704 : // We fill the 1/A and H/A functors
705 192833 : populateCouplingFunctors(_HbyA_raw, _Ainv_raw);
706 :
707 192833 : if (verbose)
708 : {
709 0 : _console << "************************************" << std::endl;
710 0 : _console << "DONE Computing HbyA " << std::endl;
711 0 : _console << "************************************" << std::endl;
712 : }
713 192833 : }
714 :
715 : std::vector<std::unique_ptr<NumericVector<Number>>> &
716 192833 : RhieChowMassFlux::selectPressureGradient(const bool updated_pressure)
717 : {
718 192833 : if (updated_pressure)
719 : {
720 171133 : _grad_p_current.clear();
721 450586 : for (const auto & component : _pressure_system->gradientContainer())
722 558906 : _grad_p_current.push_back(component->clone());
723 : }
724 :
725 192833 : return _grad_p_current;
726 : }
|
