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