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 1495 : RhieChowMassFlux::validParams()
33 : {
34 1495 : auto params = RhieChowFaceFluxProvider::validParams();
35 1495 : params += NonADFunctorInterface::validParams();
36 :
37 1495 : params.addClassDescription("Computes H/A and 1/A together with face mass fluxes for segregated "
38 : "momentum-pressure equations using linear systems.");
39 :
40 1495 : params.addRequiredParam<VariableName>(NS::pressure, "The pressure variable.");
41 2990 : params.addRequiredParam<VariableName>("u", "The x-component of velocity");
42 2990 : params.addParam<VariableName>("v", "The y-component of velocity");
43 2990 : params.addParam<VariableName>("w", "The z-component of velocity");
44 2990 : params.addRequiredParam<std::string>("p_diffusion_kernel",
45 : "The diffusion kernel acting on the pressure.");
46 2990 : 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 1495 : 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 1495 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
57 1495 : exec_enum.addAvailableFlags(EXEC_NONE);
58 4485 : exec_enum = {EXEC_NONE};
59 1495 : params.suppressParameter<ExecFlagEnum>("execute_on");
60 :
61 : // Pressure projection
62 2990 : params.addParam<MooseEnum>("pressure_projection_method",
63 4485 : MooseEnum("standard consistent", "standard"),
64 : "The method to use in the pressure projection for Ainv - "
65 : "standard (SIMPLE) or consistent (SIMPLEC)");
66 1495 : return params;
67 1495 : }
68 :
69 672 : RhieChowMassFlux::RhieChowMassFlux(const InputParameters & params)
70 : : RhieChowFaceFluxProvider(params),
71 : NonADFunctorInterface(this),
72 1344 : _moose_mesh(UserObject::_subproblem.mesh()),
73 672 : _mesh(_moose_mesh.getMesh()),
74 672 : _dim(blocksMaxDimension()),
75 672 : _p(dynamic_cast<MooseLinearVariableFVReal *>(
76 672 : &UserObject::_subproblem.getVariable(0, getParam<VariableName>(NS::pressure)))),
77 672 : _vel(_dim, nullptr),
78 1344 : _HbyA_flux(_moose_mesh, blockIDs(), "HbyA_flux"),
79 1344 : _Ainv(_moose_mesh, blockIDs(), "Ainv"),
80 672 : _face_mass_flux(
81 1344 : declareRestartableData<FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>>(
82 : "face_flux", _moose_mesh, blockIDs(), "face_values")),
83 2016 : _body_force_kernel_names(
84 : getParam<std::vector<std::vector<std::string>>>("body_force_kernel_names")),
85 672 : _rho(getFunctor<Real>(NS::density)),
86 2688 : _pressure_projection_method(getParam<MooseEnum>("pressure_projection_method"))
87 : {
88 672 : if (!_p)
89 0 : paramError(NS::pressure, "the pressure must be a MooseLinearVariableFVReal.");
90 672 : checkBlocks(*_p);
91 :
92 672 : std::vector<std::string> vel_names = {"u", "v", "w"};
93 1970 : for (const auto i : index_range(_vel))
94 : {
95 1298 : _vel[i] = dynamic_cast<MooseLinearVariableFVReal *>(
96 1298 : &UserObject::_subproblem.getVariable(0, getParam<VariableName>(vel_names[i])));
97 :
98 1298 : if (!_vel[i])
99 0 : paramError(vel_names[i], "the velocity must be a MOOSELinearVariableFVReal.");
100 1298 : checkBlocks(*_vel[i]);
101 : }
102 :
103 : // Register the elemental/face functors which will be queried in the pressure equation
104 1344 : for (const auto tid : make_range(libMesh::n_threads()))
105 : {
106 672 : UserObject::_subproblem.addFunctor("Ainv", _Ainv, tid);
107 1344 : UserObject::_subproblem.addFunctor("HbyA", _HbyA_flux, tid);
108 : }
109 :
110 672 : if (!dynamic_cast<SIMPLE *>(getMooseApp().getExecutioner()) &&
111 70 : !dynamic_cast<PIMPLE *>(getMooseApp().getExecutioner()))
112 0 : mooseError(this->name(),
113 : " should only be used with a linear segregated thermal-hydraulics solver!");
114 672 : }
115 :
116 : void
117 671 : 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 671 : _momentum_systems = momentum_systems;
123 671 : _momentum_system_numbers = momentum_system_numbers;
124 671 : _pressure_system = &pressure_system;
125 671 : _global_pressure_system_number = _pressure_system->number();
126 :
127 671 : _momentum_implicit_systems.clear();
128 1967 : for (auto & system : _momentum_systems)
129 : {
130 1296 : _global_momentum_system_numbers.push_back(system->number());
131 1296 : _momentum_implicit_systems.push_back(dynamic_cast<LinearImplicitSystem *>(&system->system()));
132 : }
133 :
134 671 : setupMeshInformation();
135 671 : }
136 :
137 : void
138 20 : RhieChowMassFlux::meshChanged()
139 : {
140 : _HbyA_flux.clear();
141 : _Ainv.clear();
142 20 : _face_mass_flux.clear();
143 20 : setupMeshInformation();
144 20 : }
145 :
146 : void
147 671 : 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 671 : auto base_query = _fe_problem.theWarehouse()
153 671 : .query()
154 671 : .template condition<AttribThread>(_tid)
155 671 : .template condition<AttribSysNum>(_p->sys().number())
156 671 : .template condition<AttribSystem>("LinearFVFluxKernel")
157 2013 : .template condition<AttribName>(getParam<std::string>("p_diffusion_kernel"))
158 671 : .queryInto(flux_kernel);
159 671 : 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 671 : _p_diffusion_kernel = dynamic_cast<LinearFVAnisotropicDiffusion *>(flux_kernel[0]);
164 671 : 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 671 : if (!_body_force_kernel_names.empty())
174 : {
175 72 : 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 72 : _body_force_kernels.resize(_dim);
180 :
181 216 : for (const auto dim_i : make_range(_dim))
182 288 : for (const auto & force_name : _body_force_kernel_names[dim_i])
183 : {
184 : std::vector<LinearFVElementalKernel *> temp_storage;
185 144 : auto base_query_force = _fe_problem.theWarehouse()
186 144 : .query()
187 144 : .template condition<AttribThread>(_tid)
188 288 : .template condition<AttribSysNum>(_vel[dim_i]->sys().number())
189 144 : .template condition<AttribSystem>("LinearFVElementalKernel")
190 144 : .template condition<AttribName>(force_name)
191 144 : .queryInto(temp_storage);
192 144 : 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 144 : _body_force_kernels[dim_i].push_back(temp_storage[0]);
197 144 : }
198 : }
199 671 : }
200 :
201 : void
202 691 : 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 1382 : _cell_volumes = _pressure_system->currentSolution()->zero_clone();
207 81606 : 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 80915 : if (hasBlocks(elem_info->subdomain_id()))
211 : {
212 77383 : const auto elem_dof = elem_info->dofIndices()[_global_pressure_system_number][0];
213 77383 : _cell_volumes->set(elem_dof, elem_info->volume() * elem_info->coordFactor());
214 : }
215 :
216 691 : _cell_volumes->close();
217 :
218 691 : _flow_face_info.clear();
219 162058 : for (auto & fi : _fe_problem.mesh().faceInfo())
220 167638 : if (hasBlocks(fi->elemPtr()->subdomain_id()) ||
221 13364 : (fi->neighborPtr() && hasBlocks(fi->neighborPtr()->subdomain_id())))
222 154387 : _flow_face_info.push_back(fi);
223 691 : }
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 652 : RhieChowMassFlux::initFaceMassFlux()
237 : {
238 : using namespace Moose::FV;
239 :
240 652 : 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 143876 : 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 143224 : if (_vel[0]->isInternalFace(*fi))
250 : {
251 126946 : const auto & elem_info = *fi->elemInfo();
252 : const auto & neighbor_info = *fi->neighborInfo();
253 :
254 126946 : Real elem_rho = _rho(makeElemArg(fi->elemPtr()), time_arg);
255 253892 : Real neighbor_rho = _rho(makeElemArg(fi->neighborPtr()), time_arg);
256 :
257 375781 : for (const auto dim_i : index_range(_vel))
258 248835 : interpolate(InterpMethod::Average,
259 : density_times_velocity(dim_i),
260 248835 : _vel[dim_i]->getElemValue(elem_info, time_arg) * elem_rho,
261 497670 : _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 16278 : const bool elem_is_fluid = hasBlocks(fi->elemPtr()->subdomain_id());
269 16278 : 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 16278 : const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
273 16278 : const Moose::FaceArg boundary_face{
274 16278 : fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr};
275 :
276 16278 : const Real face_rho = _rho(boundary_face, time_arg);
277 49732 : for (const auto dim_i : index_range(_vel))
278 33454 : density_times_velocity(dim_i) = boundary_normal_multiplier * face_rho *
279 33454 : raw_value((*_vel[dim_i])(boundary_face, time_arg));
280 : }
281 :
282 143224 : _face_mass_flux[fi->id()] = density_times_velocity * fi->normal();
283 : }
284 652 : }
285 :
286 : Real
287 57177008 : RhieChowMassFlux::getMassFlux(const FaceInfo & fi) const
288 : {
289 57177008 : return _face_mass_flux.evaluate(&fi);
290 : }
291 :
292 : Real
293 2204230 : RhieChowMassFlux::getVolumetricFaceFlux(const FaceInfo & fi) const
294 : {
295 2204230 : 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 2204230 : /*state_limiter*/ nullptr};
301 2204230 : const Real face_rho = _rho(face_arg, Moose::currentState());
302 2204230 : return libmesh_map_find(_face_mass_flux, fi.id()) / face_rho;
303 : }
304 :
305 : Real
306 2512 : 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 2512 : if (m != Moose::FV::InterpMethod::RhieChow)
315 0 : mooseError("Interpolation methods other than Rhie-Chow are not supported!");
316 2512 : if (time.state != Moose::currentState().state)
317 0 : mooseError("Older interpolation times are not supported!");
318 :
319 2512 : return getVolumetricFaceFlux(fi);
320 : }
321 :
322 : void
323 119163 : RhieChowMassFlux::computeFaceMassFlux()
324 : {
325 : using namespace Moose::FV;
326 :
327 119163 : const auto time_arg = Moose::currentState();
328 :
329 : // Petsc vector reader to make the repeated reading from the vector faster
330 119163 : 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 21406619 : for (auto & fi : _flow_face_info)
335 : {
336 : // Making sure the kernel knows which face we are on
337 21287456 : _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 21287456 : _p_diffusion_kernel->setCurrentFaceArea(1.0);
342 :
343 : Real p_grad_flux = 0.0;
344 21287456 : if (_p->isInternalFace(*fi))
345 : {
346 18780552 : const auto & elem_info = *fi->elemInfo();
347 : const auto & neighbor_info = *fi->neighborInfo();
348 :
349 : // Fetching the dof indices for the pressure variable
350 18780552 : const auto elem_dof = elem_info.dofIndices()[_global_pressure_system_number][0];
351 18780552 : 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 18780552 : const auto p_elem_value = p_reader(elem_dof);
355 18780552 : const auto p_neighbor_value = p_reader(neighbor_dof);
356 :
357 : // Compute the elem matrix contributions for the face
358 18780552 : const auto elem_matrix_contribution = _p_diffusion_kernel->computeElemMatrixContribution();
359 : const auto neighbor_matrix_contribution =
360 18780552 : _p_diffusion_kernel->computeNeighborMatrixContribution();
361 : const auto elem_rhs_contribution =
362 18780552 : _p_diffusion_kernel->computeElemRightHandSideContribution();
363 :
364 : // Compute the face flux from the matrix and right hand side contributions
365 18780552 : p_grad_flux = (p_neighbor_value * neighbor_matrix_contribution +
366 18780552 : p_elem_value * elem_matrix_contribution) -
367 : elem_rhs_contribution;
368 : }
369 2506904 : 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 1692941 : bc_pointer->setupFaceData(
374 1692941 : fi, fi->faceType(std::make_pair(_p->number(), _global_pressure_system_number)));
375 :
376 : const ElemInfo & elem_info =
377 1692941 : hasBlocks(fi->elemPtr()->subdomain_id()) ? *fi->elemInfo() : *fi->neighborInfo();
378 1692941 : const auto p_elem_value = _p->getElemValue(elem_info, time_arg);
379 : const auto matrix_contribution =
380 1692941 : _p_diffusion_kernel->computeBoundaryMatrixContribution(*bc_pointer);
381 : const auto rhs_contribution =
382 1692941 : _p_diffusion_kernel->computeBoundaryRHSContribution(*bc_pointer);
383 :
384 : // On the boundary, only the element side has a contribution
385 1692941 : p_grad_flux = (p_elem_value * matrix_contribution - rhs_contribution);
386 : }
387 : // Compute the new face flux
388 21287456 : _face_mass_flux[fi->id()] = -_HbyA_flux[fi->id()] + p_grad_flux;
389 : }
390 119163 : }
391 :
392 : void
393 132163 : RhieChowMassFlux::computeCellVelocity()
394 : {
395 132163 : auto & pressure_gradient = _pressure_system->linearFVGradientContainer();
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 360006 : for (const auto system_i : index_range(_momentum_implicit_systems))
400 : {
401 227843 : auto working_vector = _Ainv_raw[system_i]->clone();
402 227843 : working_vector->pointwise_mult(*working_vector, *pressure_gradient[system_i]);
403 227843 : working_vector->add(*_HbyA_raw[system_i]);
404 227843 : working_vector->scale(-1.0);
405 227843 : (*_momentum_implicit_systems[system_i]->solution) = *working_vector;
406 227843 : _momentum_implicit_systems[system_i]->update();
407 227843 : _momentum_systems[system_i]->setSolution(
408 227843 : *_momentum_implicit_systems[system_i]->current_local_solution);
409 227843 : }
410 132163 : }
411 :
412 : void
413 671 : 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 161296 : for (auto & fi : _fe_problem.mesh().faceInfo())
419 : {
420 160625 : _Ainv[fi->id()];
421 160625 : _HbyA_flux[fi->id()];
422 : }
423 671 : }
424 :
425 : void
426 132163 : 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 132163 : const auto time_arg = Moose::currentState();
434 :
435 : // Create the petsc vector readers for faster repeated access
436 : std::vector<PetscVectorReader> hbya_reader;
437 360006 : for (const auto dim_i : index_range(raw_hbya))
438 227843 : hbya_reader.emplace_back(*raw_hbya[dim_i]);
439 :
440 : std::vector<PetscVectorReader> ainv_reader;
441 360006 : for (const auto dim_i : index_range(raw_Ainv))
442 227843 : 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 22900219 : 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 22768056 : auto & Ainv = _Ainv[fi->id()];
452 :
453 : // If it is internal, we just interpolate (using geometric weights) to the face
454 22768056 : if (_vel[0]->isInternalFace(*fi))
455 : {
456 : // Get the dof indices for the element and the neighbor
457 19951552 : const auto & elem_info = *fi->elemInfo();
458 : const auto & neighbor_info = *fi->neighborInfo();
459 19951552 : const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
460 19951552 : 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 19951552 : const Real elem_rho = _rho(makeElemArg(fi->elemPtr()), time_arg);
465 39903104 : const Real neighbor_rho = _rho(makeElemArg(fi->neighborPtr()), time_arg);
466 :
467 : // Now we do the interpolation to the face
468 19951552 : interpolate(Moose::FV::InterpMethod::Average, face_rho, elem_rho, neighbor_rho, *fi, true);
469 59041957 : for (const auto dim_i : index_range(raw_hbya))
470 : {
471 39090405 : interpolate(Moose::FV::InterpMethod::Average,
472 : face_hbya(dim_i),
473 39090405 : hbya_reader[dim_i](elem_dof),
474 39090405 : hbya_reader[dim_i](neighbor_dof),
475 : *fi,
476 : true);
477 39090405 : interpolate(InterpMethod::Average,
478 : Ainv(dim_i),
479 39090405 : elem_rho * ainv_reader[dim_i](elem_dof),
480 78180810 : neighbor_rho * ainv_reader[dim_i](neighbor_dof),
481 : *fi,
482 : true);
483 : }
484 : }
485 : else
486 : {
487 2816504 : 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 2816504 : const Real boundary_normal_multiplier = elem_is_fluid ? 1.0 : -1.0;
491 :
492 2816504 : const ElemInfo & elem_info = elem_is_fluid ? *fi->elemInfo() : *fi->neighborInfo();
493 2816504 : 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 2816504 : if (_vel[0]->isDirichletBoundaryFace(*fi))
498 : {
499 2093398 : const Moose::FaceArg boundary_face{
500 2093398 : fi, Moose::FV::LimiterType::CentralDifference, true, false, elem_info.elem(), nullptr};
501 2093398 : face_rho = _rho(boundary_face, Moose::currentState());
502 :
503 6281963 : for (const auto dim_i : make_range(_dim))
504 : {
505 :
506 4188565 : face_hbya(dim_i) =
507 4188565 : -MetaPhysicL::raw_value((*_vel[dim_i])(boundary_face, Moose::currentState()));
508 :
509 4188565 : if (!_body_force_kernel_names.empty())
510 639744 : for (const auto & force_kernel : _body_force_kernels[dim_i])
511 : {
512 319872 : force_kernel->setCurrentElemInfo(&elem_info);
513 319872 : face_hbya(dim_i) -=
514 319872 : force_kernel->computeRightHandSideContribution() * ainv_reader[dim_i](elem_dof) /
515 319872 : (elem_info.volume() * elem_info.coordFactor()); // zero-term expansion
516 : }
517 4188565 : 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 723106 : const auto elem_dof = elem_info.dofIndices()[_global_momentum_system_numbers[0]][0];
524 :
525 723106 : face_rho = _rho(makeElemArg(elem_info.elem()), time_arg);
526 2148383 : for (const auto dim_i : make_range(_dim))
527 1425277 : 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 2816504 : const Real elem_rho = _rho(makeElemArg(elem_info.elem()), time_arg);
532 8430346 : for (const auto dim_i : index_range(raw_Ainv))
533 5613842 : Ainv(dim_i) = elem_rho * ainv_reader[dim_i](elem_dof);
534 : }
535 : // Lastly, we populate the face flux resulted by H/A
536 22768056 : _HbyA_flux[fi->id()] = face_hbya * fi->normal() * face_rho;
537 : }
538 132163 : }
539 :
540 : void
541 132163 : RhieChowMassFlux::computeHbyA(const bool with_updated_pressure, bool verbose)
542 : {
543 132163 : 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 132163 : auto & pressure_gradient = selectPressureGradient(with_updated_pressure);
553 :
554 132163 : _HbyA_raw.clear();
555 132163 : _Ainv_raw.clear();
556 :
557 360006 : for (auto system_i : index_range(_momentum_systems))
558 : {
559 227843 : LinearImplicitSystem * momentum_system = _momentum_implicit_systems[system_i];
560 :
561 227843 : NumericVector<Number> & rhs = *(momentum_system->rhs);
562 : NumericVector<Number> & current_local_solution = *(momentum_system->current_local_solution);
563 : NumericVector<Number> & solution = *(momentum_system->solution);
564 227843 : PetscMatrix<Number> * mmat = dynamic_cast<PetscMatrix<Number> *>(momentum_system->matrix);
565 : mooseAssert(mmat,
566 : "The matrices used in the segregated INSFVRhieChow objects need to be convertable "
567 : "to PetscMatrix!");
568 :
569 227843 : if (verbose)
570 : {
571 0 : _console << "Matrix in rc object" << std::endl;
572 0 : mmat->print();
573 : }
574 :
575 : // First, we extract the diagonal and we will hold on to it for a little while
576 455686 : _Ainv_raw.push_back(current_local_solution.zero_clone());
577 : NumericVector<Number> & Ainv = *(_Ainv_raw.back());
578 :
579 227843 : mmat->get_diagonal(Ainv);
580 :
581 227843 : if (verbose)
582 : {
583 0 : _console << "Velocity solution in H(u)" << std::endl;
584 0 : solution.print();
585 : }
586 :
587 : // Time to create H(u) = M_{offdiag} * u - b_{nonpressure}
588 455686 : _HbyA_raw.push_back(current_local_solution.zero_clone());
589 : NumericVector<Number> & HbyA = *(_HbyA_raw.back());
590 :
591 : // We start with the matrix product part, we will do
592 : // M*u - A*u for 2 reasons:
593 : // 1, We assume A*u petsc operation is faster than setting the matrix diagonal to 0
594 : // 2, In PISO loops we need to reuse the matrix so we can't just set the diagonals to 0
595 :
596 : // We create a working vector to ease some of the operations, we initialize its values
597 : // with the current solution values to have something for the A*u term
598 227843 : auto working_vector = momentum_system->current_local_solution->zero_clone();
599 : PetscVector<Number> * working_vector_petsc =
600 227843 : dynamic_cast<PetscVector<Number> *>(working_vector.get());
601 : mooseAssert(working_vector_petsc,
602 : "The vectors used in the RhieChowMassFlux objects need to be convertable "
603 : "to PetscVectors!");
604 :
605 227843 : mmat->vector_mult(HbyA, solution);
606 227843 : working_vector_petsc->pointwise_mult(Ainv, solution);
607 227843 : HbyA.add(-1.0, *working_vector_petsc);
608 :
609 227843 : if (verbose)
610 : {
611 0 : _console << " H(u)" << std::endl;
612 0 : HbyA.print();
613 : }
614 :
615 : // We continue by adding the momentum right hand side contributions
616 227843 : HbyA.add(-1.0, rhs);
617 :
618 : // Unfortunately, the pressure forces are included in the momentum RHS
619 : // so we have to correct them back
620 227843 : working_vector_petsc->pointwise_mult(*pressure_gradient[system_i], *_cell_volumes);
621 227843 : HbyA.add(-1.0, *working_vector_petsc);
622 :
623 227843 : if (verbose)
624 : {
625 0 : _console << "total RHS" << std::endl;
626 0 : rhs.print();
627 0 : _console << "pressure RHS" << std::endl;
628 0 : pressure_gradient[system_i]->print();
629 0 : _console << " H(u)-rhs-relaxation_source" << std::endl;
630 0 : HbyA.print();
631 : }
632 :
633 : // It is time to create element-wise 1/A-s based on the the diagonal of the momentum matrix
634 227843 : *working_vector_petsc = 1.0;
635 227843 : Ainv.pointwise_divide(*working_vector_petsc, Ainv);
636 :
637 : // Create 1/A*(H(u)-RHS)
638 227843 : HbyA.pointwise_mult(HbyA, Ainv);
639 :
640 227843 : if (verbose)
641 : {
642 0 : _console << " (H(u)-rhs)/A" << std::endl;
643 0 : HbyA.print();
644 : }
645 :
646 227843 : if (_pressure_projection_method == "consistent")
647 : {
648 :
649 : // Consistent Corrections to SIMPLE
650 : // 1. Ainv_old = 1/a_p <- Ainv = 1/(a_p + \sum_n a_n)
651 : // 2. H(u) <- H(u*) + H(u') = H(u*) - (Ainv - Ainv_old) * grad(p) * Vc
652 :
653 4962 : if (verbose)
654 0 : _console << "Performing SIMPLEC projection." << std::endl;
655 :
656 : // Lambda function to calculate the sum of diagonal and neighbor coefficients
657 4962 : auto get_row_sum = [mmat](NumericVector<Number> & sum_vector)
658 : {
659 : // Ensure the sum_vector is zeroed out
660 4962 : sum_vector.zero();
661 :
662 : // Local row size
663 4962 : const auto local_size = mmat->local_m();
664 :
665 296370 : for (const auto row_i : make_range(local_size))
666 : {
667 : // Get all non-zero components of the row of the matrix
668 291408 : const auto global_index = mmat->row_start() + row_i;
669 : std::vector<numeric_index_type> indices;
670 : std::vector<Real> values;
671 291408 : mmat->get_row(global_index, indices, values);
672 :
673 : // Sum row elements (no absolute values)
674 : const Real row_sum = std::accumulate(values.cbegin(), values.cend(), 0.0);
675 :
676 : // Add the sum of diagonal and elements to the sum_vector
677 291408 : sum_vector.add(global_index, row_sum);
678 291408 : }
679 4962 : sum_vector.close();
680 4962 : };
681 :
682 : // Create a temporary vector to store the sum of diagonal and neighbor coefficients
683 4962 : auto row_sum = current_local_solution.zero_clone();
684 4962 : get_row_sum(*row_sum);
685 :
686 : // Create vector with new inverse projection matrix
687 4962 : auto Ainv_full = current_local_solution.zero_clone();
688 4962 : *working_vector_petsc = 1.0;
689 4962 : Ainv_full->pointwise_divide(*working_vector_petsc, *row_sum);
690 4962 : const auto Ainv_full_old = Ainv_full->clone();
691 :
692 : // Correct HbyA
693 4962 : Ainv_full->add(-1.0, Ainv);
694 4962 : working_vector_petsc->pointwise_mult(*Ainv_full, *pressure_gradient[system_i]);
695 4962 : working_vector_petsc->pointwise_mult(*working_vector_petsc, *_cell_volumes);
696 4962 : HbyA.add(-1.0, *working_vector_petsc);
697 :
698 : // Correct Ainv
699 4962 : Ainv = *Ainv_full_old;
700 4962 : }
701 :
702 227843 : Ainv.pointwise_mult(Ainv, *_cell_volumes);
703 :
704 227843 : if (verbose)
705 : {
706 0 : _console << " 1/A" << std::endl;
707 0 : Ainv.print();
708 : }
709 227843 : }
710 :
711 : // We fill the 1/A and H/A functors
712 132163 : populateCouplingFunctors(_HbyA_raw, _Ainv_raw);
713 :
714 132163 : if (verbose)
715 : {
716 0 : _console << "************************************" << std::endl;
717 0 : _console << "DONE Computing HbyA " << std::endl;
718 0 : _console << "************************************" << std::endl;
719 : }
720 132163 : }
721 :
722 : std::vector<std::unique_ptr<NumericVector<Number>>> &
723 132163 : RhieChowMassFlux::selectPressureGradient(const bool updated_pressure)
724 : {
725 132163 : if (updated_pressure)
726 : {
727 119163 : _grad_p_current.clear();
728 321006 : for (const auto & component : _pressure_system->linearFVGradientContainer())
729 403686 : _grad_p_current.push_back(component->clone());
730 : }
731 :
732 132163 : return _grad_p_current;
733 : }
|
