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