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 : #include "INSFVRhieChowInterpolatorSegregated.h"
11 : #include "INSFVAttributes.h"
12 : #include "SubProblem.h"
13 : #include "MooseMesh.h"
14 : #include "NS.h"
15 : #include "Assembly.h"
16 : #include "INSFVVelocityVariable.h"
17 : #include "INSFVPressureVariable.h"
18 : #include "PiecewiseByBlockLambdaFunctor.h"
19 : #include "VectorCompositeFunctor.h"
20 : #include "SIMPLENonlinearAssembly.h"
21 :
22 : #include "libmesh/mesh_base.h"
23 : #include "libmesh/elem_range.h"
24 : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
25 :
26 : #include "NonlinearSystem.h"
27 : #include "libmesh/petsc_matrix.h"
28 : #include "libmesh/petsc_vector.h"
29 :
30 : using namespace libMesh;
31 :
32 : registerMooseObject("NavierStokesApp", INSFVRhieChowInterpolatorSegregated);
33 :
34 : InputParameters
35 1484 : INSFVRhieChowInterpolatorSegregated::validParams()
36 : {
37 1484 : auto params = RhieChowInterpolatorBase::validParams();
38 :
39 1484 : params.addClassDescription("Computes H/A and 1/A together with face velocities for segregated "
40 : "momentum-pressure equations.");
41 :
42 : // We disable the execution of this, should only provide functions
43 : // for the SIMPLENonlinearAssembly executioner
44 1484 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
45 1484 : exec_enum.addAvailableFlags(EXEC_NONE);
46 4452 : exec_enum = {EXEC_NONE};
47 1484 : params.suppressParameter<ExecFlagEnum>("execute_on");
48 :
49 1484 : return params;
50 1484 : }
51 :
52 742 : INSFVRhieChowInterpolatorSegregated::INSFVRhieChowInterpolatorSegregated(
53 742 : const InputParameters & params)
54 : : RhieChowInterpolatorBase(params),
55 1484 : _HbyA(_moose_mesh, blockIDs(), "HbyA"),
56 742 : _Ainv(_moose_mesh, blockIDs(), "Ainv", false),
57 742 : _vel(std::make_unique<PiecewiseByBlockLambdaFunctor<ADRealVectorValue>>(
58 : name(),
59 1132014 : [this](const auto & r, const auto & t) -> ADRealVectorValue
60 : {
61 2262544 : ADRealVectorValue velocity((*_u)(r, t));
62 1131272 : if (_dim >= 2)
63 1131272 : velocity(1) = (*_v)(r, t);
64 1131272 : if (_dim >= 3)
65 439349 : velocity(2) = (*_w)(r, t);
66 1131272 : return velocity;
67 : },
68 2226 : std::set<ExecFlagType>({EXEC_ALWAYS}),
69 : _moose_mesh,
70 : blockIDs())),
71 2226 : _face_velocity(_moose_mesh, blockIDs(), "face_values")
72 : {
73 742 : if (_displaced)
74 0 : paramError("use_displaced_mesh",
75 : "The segregated Rhie-Chow user object does not currently support operation on a "
76 : "displaced mesh");
77 :
78 : // Register the elemental/face functors which will be queried in the pressure equation
79 1599 : for (const auto tid : make_range(libMesh::n_threads()))
80 : {
81 857 : UserObject::_subproblem.addFunctor("Ainv", _Ainv, tid);
82 1714 : UserObject::_subproblem.addFunctor("HbyA", _HbyA, tid);
83 : }
84 :
85 742 : if (_velocity_interp_method == Moose::FV::InterpMethod::Average)
86 0 : paramError("velocity_interp_method",
87 : "Segregated momentum-pressure solvers do not allow average interpolation methods!");
88 :
89 742 : if (!dynamic_cast<SIMPLENonlinearAssembly *>(getMooseApp().getExecutioner()))
90 0 : mooseError(this->name(), " should only be used with a segregated thermal-hydraulics solver!");
91 1484 : }
92 :
93 : void
94 742 : INSFVRhieChowInterpolatorSegregated::linkMomentumSystem(
95 : std::vector<NonlinearSystemBase *> momentum_systems,
96 : const std::vector<unsigned int> & momentum_system_numbers,
97 : const TagID pressure_gradient_tag)
98 : {
99 742 : _momentum_systems = momentum_systems;
100 742 : _momentum_system_numbers = momentum_system_numbers;
101 742 : _pressure_gradient_tag = pressure_gradient_tag;
102 :
103 742 : _momentum_implicit_systems.clear();
104 2283 : for (auto & system : _momentum_systems)
105 : _momentum_implicit_systems.push_back(
106 1541 : dynamic_cast<NonlinearImplicitSystem *>(&system->system()));
107 742 : }
108 :
109 : void
110 0 : INSFVRhieChowInterpolatorSegregated::meshChanged()
111 : {
112 : _HbyA.clear();
113 : _Ainv.clear();
114 : _face_velocity.clear();
115 0 : }
116 :
117 : void
118 0 : INSFVRhieChowInterpolatorSegregated::initialize()
119 : {
120 0 : for (const auto & pair : _HbyA)
121 0 : _HbyA[pair.first] = 0;
122 :
123 0 : for (const auto & pair : _Ainv)
124 0 : _Ainv[pair.first] = 0;
125 0 : }
126 :
127 : void
128 742 : INSFVRhieChowInterpolatorSegregated::initFaceVelocities()
129 : {
130 130514 : for (auto & fi : _fe_problem.mesh().faceInfo())
131 : {
132 130567 : if (hasBlocks(fi->elemPtr()->subdomain_id()) ||
133 1802 : (fi->neighborPtr() && hasBlocks(fi->neighborPtr()->subdomain_id())))
134 : {
135 : // On internal face we do a regular interpoaltion with geometric weights
136 128773 : if (_u->isInternalFace(*fi))
137 : {
138 106914 : const Moose::FaceArg face{
139 106914 : fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr};
140 :
141 213828 : _face_velocity[fi->id()] = MetaPhysicL::raw_value((*_vel)(face, Moose::currentState()));
142 : }
143 : // On the boundary, we just take the boundary values
144 : else
145 : {
146 : const Elem * const boundary_elem =
147 21859 : hasBlocks(fi->elemPtr()->subdomain_id()) ? fi->elemPtr() : fi->neighborPtr();
148 :
149 21859 : const Moose::FaceArg boundary_face{
150 21859 : fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr};
151 :
152 21859 : _face_velocity[fi->id()] =
153 43718 : MetaPhysicL::raw_value((*_vel)(boundary_face, Moose::currentState()));
154 : }
155 : }
156 : }
157 742 : }
158 :
159 : VectorValue<ADReal>
160 13460960 : INSFVRhieChowInterpolatorSegregated::getVelocity(const Moose::FV::InterpMethod m,
161 : const FaceInfo & fi,
162 : const Moose::StateArg & /*time*/,
163 : const THREAD_ID /*tid*/,
164 : const bool /*subtract_mesh_velocity*/) const
165 : {
166 13460960 : if (m != Moose::FV::InterpMethod::RhieChow)
167 0 : mooseError("Segregated solution algorithms only support Rhie-Chow interpolation!");
168 13460960 : return _face_velocity.evaluate(&fi);
169 : }
170 :
171 : void
172 43563 : INSFVRhieChowInterpolatorSegregated::computeFaceVelocity()
173 : {
174 43563 : const auto time_arg = Moose::currentState();
175 :
176 6447382 : for (auto & fi : _fe_problem.mesh().faceInfo())
177 : {
178 6465505 : if (hasBlocks(fi->elemPtr()->subdomain_id()) ||
179 142477 : (fi->neighborPtr() && hasBlocks(fi->neighborPtr()->subdomain_id())))
180 : {
181 : // On internal face we just use the interpolated H/A and the pressure face gradient
182 : // So u_f = -(H/A)_f - (1/A)_f*grad(p)_f
183 : // Notice the (-) sign on H/A which is because we use the Jacobian/Residual
184 : // computations and we get -H instead of H.
185 6324022 : if (_u->isInternalFace(*fi))
186 : {
187 5008774 : const Moose::FaceArg face{
188 5008774 : fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr};
189 :
190 : RealVectorValue Ainv;
191 5008774 : RealVectorValue HbyA = MetaPhysicL::raw_value(_HbyA(face, time_arg));
192 :
193 5008774 : interpolate(Moose::FV::InterpMethod::Average,
194 : Ainv,
195 5008774 : _Ainv(makeElemArg(fi->elemPtr()), time_arg),
196 10017548 : _Ainv(makeElemArg(fi->neighborPtr()), time_arg),
197 : *fi,
198 : true);
199 :
200 5008774 : RealVectorValue grad_p = MetaPhysicL::raw_value(_p->gradient(face, time_arg));
201 15725403 : for (const auto comp_index : make_range(_dim))
202 10716629 : _face_velocity[fi->id()](comp_index) =
203 10716629 : -HbyA(comp_index) - Ainv(comp_index) * grad_p(comp_index);
204 : }
205 : else
206 : {
207 : const Elem * const boundary_elem =
208 1315248 : hasBlocks(fi->elemPtr()->subdomain_id()) ? fi->elemPtr() : fi->neighborPtr();
209 1315248 : const Moose::FaceArg boundary_face{
210 1315248 : fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr};
211 :
212 : // If we have a dirichlet boundary conditions, this sill give us the exact value of the
213 : // velocity on the face as expected (see populateHbyA())
214 1315248 : if (_u->isDirichletBoundaryFace(*fi, boundary_elem, time_arg))
215 1791170 : _face_velocity[fi->id()] = -MetaPhysicL::raw_value(_HbyA(boundary_face, time_arg));
216 : else
217 : {
218 419663 : const RealVectorValue & Ainv = MetaPhysicL::raw_value(_Ainv(boundary_face, time_arg));
219 419663 : const RealVectorValue & HbyA = MetaPhysicL::raw_value(_HbyA(boundary_face, time_arg));
220 : const RealVectorValue & grad_p =
221 419663 : MetaPhysicL::raw_value(_p->gradient(boundary_face, time_arg));
222 1288392 : for (const auto comp_index : make_range(_dim))
223 868729 : _face_velocity[fi->id()](comp_index) =
224 868729 : -HbyA(comp_index) - Ainv(comp_index) * grad_p(comp_index);
225 : }
226 : }
227 : }
228 : }
229 43563 : }
230 :
231 : void
232 43563 : INSFVRhieChowInterpolatorSegregated::computeCellVelocity()
233 : {
234 43563 : std::vector<unsigned int> var_nums = {_momentum_implicit_systems[0]->variable_number(_u->name())};
235 43563 : if (_v)
236 43563 : var_nums.push_back(_momentum_implicit_systems[1]->variable_number(_v->name()));
237 43563 : if (_w)
238 5082 : var_nums.push_back(_momentum_implicit_systems[2]->variable_number(_w->name()));
239 :
240 43563 : const auto time_arg = Moose::currentState();
241 :
242 : for (auto & elem :
243 5412363 : as_range(_mesh.active_local_elements_begin(), _mesh.active_local_elements_end()))
244 : {
245 2640837 : if (hasBlocks(elem->subdomain_id()))
246 : {
247 2603296 : const auto elem_arg = makeElemArg(elem);
248 2603296 : const RealVectorValue Ainv = _Ainv(elem_arg, time_arg);
249 2603296 : const RealVectorValue & grad_p = raw_value(_p->gradient(elem_arg, time_arg));
250 :
251 8101974 : for (auto comp_index : make_range(_dim))
252 : {
253 : // If we are doing segregated momentum components we need to access different vector
254 : // components otherwise everything is in the same vector (with different variable names)
255 5498678 : const unsigned int system_number = _momentum_implicit_systems[comp_index]->number();
256 5498678 : const auto index = elem->dof_number(system_number, var_nums[comp_index], 0);
257 :
258 : // We set the dof value in the solution vector the same logic applies:
259 : // u_C = -(H/A)_C - (1/A)_C*grad(p)_C where C is the cell index
260 10997356 : _momentum_implicit_systems[comp_index]->solution->set(
261 5498678 : index, -(*_HbyA_raw[comp_index])(index)-Ainv(comp_index) * grad_p(comp_index));
262 : }
263 : }
264 43563 : }
265 :
266 135771 : for (auto system_i : index_range(_momentum_implicit_systems))
267 : {
268 92208 : _momentum_implicit_systems[system_i]->solution->close();
269 92208 : _momentum_implicit_systems[system_i]->update();
270 92208 : _momentum_systems[system_i]->setSolution(
271 92208 : *_momentum_implicit_systems[system_i]->current_local_solution);
272 : }
273 43563 : }
274 :
275 : void
276 43563 : INSFVRhieChowInterpolatorSegregated::populateHbyA(
277 : const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_hbya,
278 : const std::vector<unsigned int> & var_nums)
279 : {
280 6447382 : for (auto & fi : _fe_problem.mesh().faceInfo())
281 : {
282 6465505 : if (hasBlocks(fi->elemPtr()->subdomain_id()) ||
283 142477 : (fi->neighborPtr() && hasBlocks(fi->neighborPtr()->subdomain_id())))
284 : {
285 : // If we are on an internal face, we just interpolate the values to the faces.
286 : // Otherwise, depending on the boundary type, we take the velocity value or
287 : // extrapolated HbyA values.
288 6324022 : if (_u->isInternalFace(*fi))
289 : {
290 5008774 : const Elem * elem = fi->elemPtr();
291 : const Elem * neighbor = fi->neighborPtr();
292 15725403 : for (auto comp_index : make_range(_dim))
293 : {
294 10716629 : unsigned int system_number = _momentum_implicit_systems[comp_index]->number();
295 10716629 : const auto dof_index_elem = elem->dof_number(system_number, var_nums[comp_index], 0);
296 : const auto dof_index_neighbor =
297 10716629 : neighbor->dof_number(system_number, var_nums[comp_index], 0);
298 :
299 10716629 : interpolate(Moose::FV::InterpMethod::Average,
300 21433258 : _HbyA[fi->id()](comp_index),
301 10716629 : (*raw_hbya[comp_index])(dof_index_elem),
302 21433258 : (*raw_hbya[comp_index])(dof_index_neighbor),
303 : *fi,
304 : true);
305 : }
306 : }
307 : else
308 : {
309 : const Elem * const boundary_elem =
310 1315248 : hasBlocks(fi->elemPtr()->subdomain_id()) ? fi->elemPtr() : fi->neighborPtr();
311 :
312 1315248 : const Moose::FaceArg boundary_face{
313 1315248 : fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr};
314 :
315 1315248 : if (_u->isDirichletBoundaryFace(*fi, boundary_elem, Moose::currentState()))
316 1791170 : _HbyA[fi->id()] = -MetaPhysicL::raw_value((*_vel)(boundary_face, Moose::currentState()));
317 : else
318 1288392 : for (const auto comp_index : make_range(_dim))
319 : {
320 868729 : unsigned int system_number = _momentum_implicit_systems[comp_index]->number();
321 : const auto dof_index_elem =
322 868729 : boundary_elem->dof_number(system_number, var_nums[comp_index], 0);
323 868729 : _HbyA[fi->id()](comp_index) = (*raw_hbya[comp_index])(dof_index_elem);
324 : }
325 : }
326 : }
327 : }
328 43563 : }
329 :
330 : void
331 43563 : INSFVRhieChowInterpolatorSegregated::computeHbyA(bool verbose)
332 : {
333 43563 : if (verbose)
334 : {
335 0 : _console << "************************************" << std::endl;
336 0 : _console << "Computing HbyA" << std::endl;
337 0 : _console << "************************************" << std::endl;
338 : }
339 : mooseAssert(_momentum_implicit_systems.size() && _momentum_implicit_systems[0],
340 : "The momentum system shall be linked before calling this function!");
341 :
342 : NonlinearImplicitSystem * momentum_system = _momentum_implicit_systems[0];
343 :
344 43563 : std::vector<unsigned int> var_nums = {_momentum_implicit_systems[0]->variable_number(_u->name())};
345 43563 : if (_v)
346 43563 : var_nums.push_back(_momentum_implicit_systems[1]->variable_number(_v->name()));
347 43563 : if (_w)
348 5082 : var_nums.push_back(_momentum_implicit_systems[2]->variable_number(_w->name()));
349 :
350 43563 : _HbyA_raw.clear();
351 135771 : for (auto system_i : index_range(_momentum_systems))
352 : {
353 92208 : _fe_problem.setCurrentNonlinearSystem(_momentum_system_numbers[system_i]);
354 :
355 92208 : momentum_system = _momentum_implicit_systems[system_i];
356 :
357 92208 : NumericVector<Number> & rhs = *(momentum_system->rhs);
358 : NumericVector<Number> & current_local_solution = *(momentum_system->current_local_solution);
359 : NumericVector<Number> & solution = *(momentum_system->solution);
360 92208 : PetscMatrix<Number> * mmat = dynamic_cast<PetscMatrix<Number> *>(momentum_system->matrix);
361 : mooseAssert(mmat,
362 : "The matrices used in the segregated INSFVRhieChow objects need to be convertable "
363 : "to PetscMAtrix!");
364 :
365 92208 : if (verbose)
366 : {
367 0 : _console << "Matrix in rc object" << std::endl;
368 0 : mmat->print();
369 : }
370 :
371 92208 : auto Ainv = current_local_solution.zero_clone();
372 92208 : PetscVector<Number> * Ainv_petsc = dynamic_cast<PetscVector<Number> *>(Ainv.get());
373 :
374 92208 : mmat->get_diagonal(*Ainv_petsc);
375 :
376 92208 : auto working_vector = momentum_system->current_local_solution->zero_clone();
377 : PetscVector<Number> * working_vector_petsc =
378 92208 : dynamic_cast<PetscVector<Number> *>(working_vector.get());
379 : mooseAssert(working_vector_petsc,
380 : "The vectors used in the segregated INSFVRhieChow objects need to be convertable "
381 : "to PetscVectors!");
382 :
383 92208 : *working_vector_petsc = 1.0;
384 92208 : Ainv_petsc->pointwise_divide(*working_vector_petsc, *Ainv_petsc);
385 :
386 184416 : _HbyA_raw.push_back(current_local_solution.zero_clone());
387 : NumericVector<Number> & HbyA = *(_HbyA_raw.back());
388 :
389 92208 : if (verbose)
390 : {
391 0 : _console << "Velocity solution in H(u)" << std::endl;
392 0 : solution.print();
393 : }
394 :
395 : // We fill the 1/A functor
396 : auto active_local_begin =
397 92208 : _mesh.evaluable_elements_begin(momentum_system->get_dof_map(), var_nums[system_i]);
398 : auto active_local_end =
399 92208 : _mesh.evaluable_elements_end(momentum_system->get_dof_map(), var_nums[system_i]);
400 :
401 92208 : const auto & state = Moose::currentState();
402 13868740 : for (auto it = active_local_begin; it != active_local_end; ++it)
403 : {
404 6888266 : const Elem * elem = *it;
405 6888266 : if (this->hasBlocks(elem->subdomain_id()))
406 : {
407 6779088 : const Moose::ElemArg & elem_arg = makeElemArg(elem);
408 : Real coord_multiplier;
409 6779088 : const auto coord_type = _fe_problem.getCoordSystem(elem->subdomain_id());
410 : const unsigned int rz_radial_coord =
411 6779088 : Moose::COORD_RZ ? _fe_problem.getAxisymmetricRadialCoord() : libMesh::invalid_uint;
412 :
413 6779088 : MooseMeshUtils::coordTransformFactor(
414 6779088 : elem->vertex_average(), coord_multiplier, coord_type, rz_radial_coord);
415 :
416 6779088 : const Real volume = _moose_mesh.elemInfo(elem->id()).volume();
417 6779088 : const auto dof_index = elem->dof_number(momentum_system->number(), var_nums[system_i], 0);
418 13558176 : _Ainv[elem->id()](system_i) = MetaPhysicL::raw_value(epsilon(_tid)(elem_arg, state)) *
419 6779088 : (*Ainv_petsc)(dof_index)*volume * coord_multiplier;
420 : }
421 : }
422 :
423 : // Now we set the diagonal of our system matrix to 0 so we can create H*u
424 : // TODO: Add a function for this in libmesh
425 92208 : *working_vector_petsc = 0.0;
426 92208 : LibmeshPetscCall(MatDiagonalSet(mmat->mat(), working_vector_petsc->vec(), INSERT_VALUES));
427 :
428 92208 : if (verbose)
429 : {
430 0 : _console << "H" << std::endl;
431 0 : mmat->print();
432 : }
433 :
434 : // We need to subtract the contribution of the pressure gradient from the residual (right hand
435 : // side). We plug working_vector=0 in here to get the right hand side contribution
436 92208 : _fe_problem.computeResidualTag(*working_vector, HbyA, _pressure_gradient_tag);
437 :
438 : // Now we reorganize the association between tags and vectors to make sure
439 : // that the next solve is perfect
440 92208 : _momentum_systems[system_i]->associateVectorToTag(
441 184416 : momentum_system->get_vector(_fe_problem.vectorTagName(0)), 0);
442 92208 : _momentum_systems[system_i]->associateVectorToTag(
443 184416 : momentum_system->get_vector(_fe_problem.vectorTagName(_pressure_gradient_tag)),
444 : _pressure_gradient_tag);
445 92208 : _momentum_systems[system_i]->setSolution(current_local_solution);
446 :
447 92208 : if (verbose)
448 : {
449 0 : _console << "total RHS" << std::endl;
450 0 : rhs.print();
451 0 : _console << "pressure RHS" << std::endl;
452 0 : HbyA.print();
453 : }
454 :
455 : // We correct the right hand side to exclude the pressure contribution
456 92208 : HbyA.scale(-1.0);
457 92208 : HbyA.add(-1.0, rhs);
458 :
459 92208 : if (verbose)
460 : {
461 0 : _console << "H RHS" << std::endl;
462 0 : HbyA.print();
463 : }
464 :
465 : // Create H(u)
466 92208 : mmat->vector_mult(*working_vector_petsc, solution);
467 :
468 92208 : if (verbose)
469 : {
470 0 : _console << " H(u)" << std::endl;
471 0 : working_vector_petsc->print();
472 : }
473 :
474 : // Create H(u) - RHS
475 92208 : HbyA.add(*working_vector_petsc);
476 :
477 92208 : if (verbose)
478 : {
479 0 : _console << " H(u)-rhs-relaxation_source" << std::endl;
480 0 : HbyA.print();
481 : }
482 :
483 : // Create 1/A*(H(u)-RHS)
484 92208 : HbyA.pointwise_mult(HbyA, *Ainv);
485 :
486 92208 : if (verbose)
487 : {
488 0 : _console << " (H(u)-rhs)/A" << std::endl;
489 0 : HbyA.print();
490 : }
491 92208 : }
492 :
493 43563 : populateHbyA(_HbyA_raw, var_nums);
494 :
495 43563 : if (verbose)
496 : {
497 0 : _console << "************************************" << std::endl;
498 0 : _console << "DONE Computing HbyA " << std::endl;
499 0 : _console << "************************************" << std::endl;
500 : }
501 43563 : }
|
