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 "INSFVRhieChowInterpolator.h"
11 : #include "INSFVAttributes.h"
12 : #include "GatherRCDataElementThread.h"
13 : #include "GatherRCDataFaceThread.h"
14 : #include "MooseMesh.h"
15 : #include "SystemBase.h"
16 : #include "NS.h"
17 : #include "Assembly.h"
18 : #include "INSFVVelocityVariable.h"
19 : #include "PiecewiseByBlockLambdaFunctor.h"
20 : #include "VectorCompositeFunctor.h"
21 : #include "FVElementalKernel.h"
22 : #include "NSFVUtils.h"
23 : #include "DisplacedProblem.h"
24 :
25 : #include "libmesh/mesh_base.h"
26 : #include "libmesh/elem_range.h"
27 : #include "libmesh/parallel_algebra.h"
28 : #include "libmesh/remote_elem.h"
29 : #include "metaphysicl/metaphysicl_version.h"
30 : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
31 : #include "metaphysicl/parallel_dualnumber.h"
32 : #if METAPHYSICL_MAJOR_VERSION < 2
33 : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
34 : #else
35 : #include "metaphysicl/parallel_dynamic_array_wrapper.h"
36 : #endif
37 : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
38 : #include "timpi/parallel_sync.h"
39 :
40 : using namespace libMesh;
41 :
42 : registerMooseObject("NavierStokesApp", INSFVRhieChowInterpolator);
43 :
44 : InputParameters
45 63545 : INSFVRhieChowInterpolator::uniqueParams()
46 : {
47 63545 : auto params = emptyInputParameters();
48 127090 : params.addParam<bool>(
49 : "pull_all_nonlocal_a",
50 127090 : false,
51 : "Whether to pull all nonlocal 'a' coefficient data to our process. Note that 'nonlocal' "
52 : "means elements that we have access to (this may not be all the elements in the mesh if the "
53 : "mesh is distributed) but that we do not own.");
54 127090 : params.addParamNamesToGroup("pull_all_nonlocal_a", "Parallel Execution Tuning");
55 :
56 127090 : params.addParam<bool>(
57 127090 : "correct_volumetric_force", false, "Flag to activate volume force corrections.");
58 : MooseEnum volume_force_correction_method("force-consistent pressure-consistent",
59 127090 : "force-consistent");
60 127090 : params.addParam<MooseEnum>(
61 : "volume_force_correction_method",
62 : volume_force_correction_method,
63 : "The method used for correcting the Rhie-Chow coefficients for a volume force.");
64 127090 : params.addParam<std::vector<MooseFunctorName>>(
65 : "volumetric_force_functors", "The names of the functors with the volumetric force sources.");
66 63545 : return params;
67 63545 : }
68 :
69 : std::vector<std::string>
70 1726 : INSFVRhieChowInterpolator::listOfCommonParams()
71 : {
72 : return {"pull_all_nonlocal_a",
73 : "correct_volumetric_force",
74 : "volume_force_correction_method",
75 1726 : "volumetric_force_functors"};
76 : }
77 :
78 : InputParameters
79 14604 : INSFVRhieChowInterpolator::validParams()
80 : {
81 14604 : auto params = RhieChowInterpolatorBase::validParams();
82 14604 : params += INSFVRhieChowInterpolator::uniqueParams();
83 :
84 14604 : params.addClassDescription(
85 : "Computes the Rhie-Chow velocity based on gathered 'a' coefficient data.");
86 :
87 14604 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
88 14604 : exec_enum.addAvailableFlags(EXEC_PRE_KERNELS);
89 43812 : exec_enum = {EXEC_PRE_KERNELS};
90 14604 : params.suppressParameter<ExecFlagEnum>("execute_on");
91 :
92 29208 : params.addParam<MooseFunctorName>(
93 : "a_u",
94 : "For simulations in which the advecting velocities are aux variables, this parameter must be "
95 : "supplied. It represents the on-diagonal coefficients for the 'x' component velocity, solved "
96 : "via the Navier-Stokes equations.");
97 29208 : params.addParam<MooseFunctorName>(
98 : "a_v",
99 : "For simulations in which the advecting velocities are aux variables, this parameter must be "
100 : "supplied when the mesh dimension is greater than 1. It represents the on-diagonal "
101 : "coefficients for the 'y' component velocity, solved via the Navier-Stokes equations.");
102 29208 : params.addParam<MooseFunctorName>(
103 : "a_w",
104 : "For simulations in which the advecting velocities are aux variables, this parameter must be "
105 : "supplied when the mesh dimension is greater than 2. It represents the on-diagonal "
106 : "coefficients for the 'z' component velocity, solved via the Navier-Stokes equations.");
107 29208 : params.addParam<VariableName>("disp_x", "The x-component of displacement");
108 29208 : params.addParam<VariableName>("disp_y", "The y-component of displacement");
109 29208 : params.addParam<VariableName>("disp_z", "The z-component of displacement");
110 14604 : return params;
111 14604 : }
112 :
113 7306 : INSFVRhieChowInterpolator::INSFVRhieChowInterpolator(const InputParameters & params)
114 : : RhieChowInterpolatorBase(params),
115 7302 : _vel(libMesh::n_threads()),
116 7302 : _a(_moose_mesh, blockIDs(), "a", /*extrapolated_boundary*/ true),
117 7302 : _ax(_a, 0),
118 7302 : _ay(_a, 1),
119 7302 : _az(_a, 2),
120 14604 : _momentum_sys_number(_fe_problem.systemNumForVariable(getParam<VariableName>("u"))),
121 7302 : _example(0),
122 7302 : _a_data_provided(false),
123 14604 : _pull_all_nonlocal(getParam<bool>("pull_all_nonlocal_a")),
124 14604 : _bool_correct_vf(getParam<bool>("correct_volumetric_force")),
125 14604 : _volume_force_correction_method(getParam<MooseEnum>("volume_force_correction_method")),
126 7302 : _volumetric_force_functors(
127 14604 : isParamValid("volumetric_force_functors")
128 7352 : ? &getParam<std::vector<MooseFunctorName>>("volumetric_force_functors")
129 14608 : : nullptr)
130 : {
131 254 : auto process_displacement = [this](const auto & disp_name, auto & disp_container)
132 : {
133 254 : if (!_displaced)
134 0 : paramError(disp_name,
135 : "Displacement provided but we are not running on the displaced mesh. If you "
136 : "really want this object to run on the displaced mesh, then set "
137 : "'use_displaced_mesh = true', otherwise remove this displacement parameter");
138 254 : disp_container.resize(libMesh::n_threads());
139 254 : fillContainer(disp_name, disp_container);
140 254 : checkBlocks(*disp_container[0]);
141 254 : };
142 :
143 14604 : if (isParamValid("disp_x"))
144 197 : process_displacement("disp_x", _disp_xs);
145 :
146 7302 : if (_dim >= 2)
147 : {
148 13064 : if (isParamValid("disp_y"))
149 57 : process_displacement("disp_y", _disp_ys);
150 12950 : else if (isParamValid("disp_x"))
151 0 : paramError("disp_y", "If 'disp_x' is provided, then 'disp_y' must be as well");
152 : }
153 :
154 7302 : if (_dim >= 3)
155 : {
156 114 : if (isParamValid("disp_z"))
157 0 : process_displacement("disp_z", _disp_zs);
158 114 : else if (isParamValid("disp_x"))
159 0 : paramError("disp_z", "If 'disp_x' is provided, then 'disp_z' must be as well");
160 : }
161 :
162 15359 : for (const auto tid : make_range(libMesh::n_threads()))
163 : {
164 16114 : _vel[tid] = std::make_unique<PiecewiseByBlockLambdaFunctor<ADRealVectorValue>>(
165 8057 : name() + std::to_string(tid),
166 3902970 : [this, tid](const auto & r, const auto & t) -> ADRealVectorValue
167 : {
168 7773712 : ADRealVectorValue velocity((*_us[tid])(r, t));
169 3886856 : if (_dim >= 2)
170 3877094 : velocity(1) = (*_vs[tid])(r, t);
171 3886856 : if (_dim >= 3)
172 49770 : velocity(2) = (*_ws[tid])(r, t);
173 3886856 : return velocity;
174 : },
175 24171 : std::set<ExecFlagType>({EXEC_ALWAYS}),
176 : _moose_mesh,
177 : blockIDs());
178 :
179 8057 : if (_disp_xs.size())
180 618 : _disps.push_back(std::make_unique<Moose::VectorCompositeFunctor<ADReal>>(
181 412 : name() + "_disp_" + std::to_string(tid),
182 : *_disp_xs[tid],
183 66 : _dim >= 2 ? static_cast<const Moose::FunctorBase<ADReal> &>(*_disp_ys[tid])
184 : : static_cast<const Moose::FunctorBase<ADReal> &>(_zero_functor),
185 206 : _dim >= 3 ? static_cast<const Moose::FunctorBase<ADReal> &>(*_disp_zs[tid])
186 : : static_cast<const Moose::FunctorBase<ADReal> &>(_zero_functor)));
187 : }
188 :
189 7344 : if (_velocity_interp_method == Moose::FV::InterpMethod::Average && isParamValid("a_u"))
190 2 : paramError("a_u",
191 : "Rhie Chow coefficients may not be specified for average velocity interpolation");
192 :
193 7300 : if (_velocity_interp_method != Moose::FV::InterpMethod::Average)
194 7281 : fillARead();
195 :
196 7296 : if (_bool_correct_vf && !_volumetric_force_functors)
197 0 : paramError("volumetric_force_functors",
198 : "At least one volumetric force functor must be specified if "
199 : "'correct_volumetric_force' is true.");
200 :
201 : // Volume correction related
202 7296 : if (_bool_correct_vf)
203 : {
204 25 : const unsigned int num_volume_forces = (*_volumetric_force_functors).size();
205 25 : _volumetric_force.resize(num_volume_forces);
206 55 : for (const auto i : make_range(num_volume_forces))
207 30 : _volumetric_force[i] = &getFunctor<Real>((*_volumetric_force_functors)[i]);
208 : }
209 15353 : }
210 :
211 : void
212 7281 : INSFVRhieChowInterpolator::fillARead()
213 : {
214 7281 : _a_read.resize(libMesh::n_threads());
215 :
216 14562 : if (isParamValid("a_u"))
217 : {
218 152 : if (_dim > 1 && !isParamValid("a_v"))
219 0 : mooseError("If a_u is provided, then a_v must be provided");
220 :
221 38 : if (_dim > 2 && !isParamValid("a_w"))
222 0 : mooseError("If a_u is provided, then a_w must be provided");
223 :
224 38 : _a_data_provided = true;
225 38 : _a_aux.resize(libMesh::n_threads());
226 : }
227 14486 : else if (isParamValid("a_v"))
228 2 : paramError("a_v", "If the a_v coefficients are provided, then a_u must be provided");
229 14482 : else if (isParamValid("a_w"))
230 2 : paramError("a_w", "If the a_w coefficients are provided, then a_u must be provided");
231 :
232 7277 : if (_a_data_provided)
233 : {
234 82 : for (const auto tid : make_range(libMesh::n_threads()))
235 : {
236 : const Moose::FunctorBase<ADReal> *v_comp, *w_comp;
237 44 : if (_dim > 1)
238 88 : v_comp = &UserObject::_subproblem.getFunctor<ADReal>(
239 88 : deduceFunctorName("a_v"), tid, name(), true);
240 : else
241 0 : v_comp = &_zero_functor;
242 44 : if (_dim > 2)
243 0 : w_comp = &UserObject::_subproblem.getFunctor<ADReal>(
244 0 : deduceFunctorName("a_w"), tid, name(), true);
245 : else
246 44 : w_comp = &_zero_functor;
247 :
248 44 : _a_aux[tid] = std::make_unique<Moose::VectorCompositeFunctor<ADReal>>(
249 : "RC_a_coeffs",
250 88 : UserObject::_subproblem.getFunctor<ADReal>(deduceFunctorName("a_u"), tid, name(), true),
251 : *v_comp,
252 44 : *w_comp);
253 44 : _a_read[tid] = _a_aux[tid].get();
254 : }
255 : }
256 : else
257 15221 : for (const auto tid : make_range(libMesh::n_threads()))
258 : {
259 7982 : _a_read[tid] = &_a;
260 :
261 : // We are the fluid flow application, so we should make sure users have the ability to
262 : // write 'a' out to aux variables for possible transfer to other applications
263 7982 : UserObject::_subproblem.addFunctor("ax", _ax, tid);
264 7982 : UserObject::_subproblem.addFunctor("ay", _ay, tid);
265 15964 : UserObject::_subproblem.addFunctor("az", _az, tid);
266 : }
267 7277 : }
268 :
269 : void
270 7216 : INSFVRhieChowInterpolator::initialSetup()
271 : {
272 7216 : insfvSetup();
273 :
274 7216 : if (_velocity_interp_method == Moose::FV::InterpMethod::Average)
275 : return;
276 20870 : for (const auto var_num : _var_numbers)
277 : {
278 : std::vector<MooseObject *> var_objects;
279 13677 : _fe_problem.theWarehouse()
280 13677 : .query()
281 13677 : .template condition<AttribVar>(static_cast<int>(var_num))
282 13677 : .template condition<AttribResidualObject>(true)
283 27354 : .template condition<AttribSysNum>(_u->sys().number())
284 : .queryInto(var_objects);
285 79533 : for (auto * const var_object : var_objects)
286 : {
287 : // Allow FVElementalKernel that are not INSFVMomentumResidualObject for now, refs #20695
288 65858 : if (!dynamic_cast<INSFVMomentumResidualObject *>(var_object) &&
289 2 : !dynamic_cast<FVElementalKernel *>(var_object))
290 0 : mooseError("Object ",
291 : var_object->name(),
292 : " is not a INSFVMomentumResidualObject. Make sure that all the objects applied "
293 : "to the momentum equation are INSFV or derived objects.");
294 65858 : else if (!dynamic_cast<INSFVMomentumResidualObject *>(var_object) &&
295 2 : dynamic_cast<FVElementalKernel *>(var_object))
296 2 : mooseWarning(
297 : "Elemental kernel ",
298 : var_object->name(),
299 : " is not a INSFVMomentumResidualObject. Make sure that all the objects applied "
300 : "to the momentum equation are INSFV or derived objects.");
301 : }
302 :
303 13675 : if (var_objects.size() == 0 && !_a_data_provided)
304 2 : mooseError("No INSFVKernels detected for the velocity variables. If you are trying to use "
305 : "auxiliary variables for advection, please specify the a_u/v/w coefficients. If "
306 : "not, please specify INSFVKernels for the momentum equations.");
307 13673 : }
308 :
309 : // Get baseline force if force-correction method is used for volumetric correction
310 7193 : if (_bool_correct_vf && _volume_force_correction_method == "force-consistent")
311 : {
312 15 : _baseline_volume_force = 1e10;
313 15 : for (const auto & loc_elem : *_elem_range)
314 : {
315 : Real elem_value = 0.0;
316 35 : for (const auto i : make_range(_volumetric_force.size()))
317 20 : elem_value += (*_volumetric_force[i])(makeElemArg(loc_elem), determineState());
318 :
319 15 : if (std::abs(elem_value) < _baseline_volume_force)
320 15 : _baseline_volume_force = std::abs(elem_value);
321 15 : if (_baseline_volume_force == 0)
322 : break;
323 : }
324 15 : _communicator.min(_baseline_volume_force);
325 : }
326 : }
327 :
328 : void
329 7216 : INSFVRhieChowInterpolator::insfvSetup()
330 : {
331 : _elem_range =
332 28864 : std::make_unique<ConstElemRange>(_mesh.active_local_subdomain_set_elements_begin(blockIDs()),
333 14432 : _mesh.active_local_subdomain_set_elements_end(blockIDs()));
334 7216 : }
335 :
336 : void
337 0 : INSFVRhieChowInterpolator::meshChanged()
338 : {
339 0 : insfvSetup();
340 :
341 : // If the mesh has been modified:
342 : // - the boundary elements may have changed
343 : // - some elements may have been refined
344 : _elements_to_push_pull.clear();
345 : _a.clear();
346 0 : }
347 :
348 : void
349 95776 : INSFVRhieChowInterpolator::initialize()
350 : {
351 : if (!needAComputation())
352 : return;
353 :
354 : // Reset map of coefficients to zero.
355 : // The keys should not have changed unless the mesh has changed
356 : // Dont reset if not in current system
357 : // IDEA: clear them derivatives
358 95594 : if (_momentum_sys_number == _fe_problem.currentNlSysNum())
359 22074856 : for (auto & pair : _a)
360 : pair.second = 0;
361 : else
362 0 : for (auto & pair : _a)
363 : {
364 : auto & a_val = pair.second;
365 0 : a_val = MetaPhysicL::raw_value(a_val);
366 : }
367 : }
368 :
369 : void
370 95776 : INSFVRhieChowInterpolator::execute()
371 : {
372 : // Either we provided the RC coefficients using aux-variable, or we are solving for another
373 : // system than the momentum equations are in, in a multi-system setup for example
374 95776 : if (_fe_problem.currentNlSysNum() != _momentum_sys_number)
375 : return;
376 :
377 : mooseAssert(!_a_data_provided,
378 : "a-coefficient data should not be provided if the velocity variables are in the "
379 : "nonlinear system and we are running kernels that compute said a-coefficients");
380 : // One might think that we should do a similar assertion for
381 : // (_velocity_interp_method == Moose::FV::InterpMethod::RhieChow). However, even if we are not
382 : // using the generated a-coefficient data in that case, some kernels have been optimized to
383 : // add their residuals into the global system during the generation of the a-coefficient data.
384 : // Hence if we were to skip the kernel execution we would drop those residuals
385 :
386 191384 : TIME_SECTION("execute", 1, "Computing Rhie-Chow coefficients");
387 :
388 : // A lot of RC data gathering leverages the automatic differentiation system, e.g. for linear
389 : // operators we pull out the 'a' coefficients by querying the ADReal residual derivatives
390 : // member at the element or neighbor dof locations. Consequently we need to enable derivative
391 : // computation. We do this here outside the threaded regions
392 95692 : const auto saved_do_derivatives = ADReal::do_derivatives;
393 95692 : ADReal::do_derivatives = true;
394 :
395 : PARALLEL_TRY
396 : {
397 95692 : GatherRCDataElementThread et(_fe_problem, _momentum_sys_number, _var_numbers);
398 95692 : Threads::parallel_reduce(*_elem_range, et);
399 95692 : }
400 95692 : PARALLEL_CATCH;
401 :
402 : PARALLEL_TRY
403 : {
404 : using FVRange = StoredRange<MooseMesh::const_face_info_iterator, const FaceInfo *>;
405 : GatherRCDataFaceThread<FVRange> fvr(
406 95692 : _fe_problem, _momentum_sys_number, _var_numbers, _displaced);
407 382768 : FVRange faces(_moose_mesh.ownedFaceInfoBegin(), _moose_mesh.ownedFaceInfoEnd());
408 95692 : Threads::parallel_reduce(faces, fvr);
409 95692 : }
410 95692 : PARALLEL_CATCH;
411 :
412 95692 : ADReal::do_derivatives = saved_do_derivatives;
413 : }
414 :
415 : void
416 95776 : INSFVRhieChowInterpolator::finalize()
417 : {
418 95594 : if (!needAComputation() || this->n_processors() == 1)
419 22716 : return;
420 :
421 : // If advecting with auxiliary variables, no need to synchronize data
422 : // Same if not solving for the velocity variables at the moment
423 73060 : if (_fe_problem.currentNlSysNum() != _momentum_sys_number)
424 : return;
425 :
426 : using Datum = std::pair<dof_id_type, VectorValue<ADReal>>;
427 : std::unordered_map<processor_id_type, std::vector<Datum>> push_data;
428 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> pull_requests;
429 78138 : static const VectorValue<ADReal> example;
430 :
431 : // Create push data
432 1420568 : for (const auto * const elem : _elements_to_push_pull)
433 : {
434 673754 : const auto id = elem->id();
435 673754 : const auto pid = elem->processor_id();
436 : auto it = _a.find(id);
437 : mooseAssert(it != _a.end(), "We definitely should have found something");
438 673754 : push_data[pid].push_back(std::make_pair(id, it->second));
439 : }
440 :
441 : // Create pull data
442 73060 : if (_pull_all_nonlocal)
443 : {
444 : for (const auto * const elem :
445 7210 : as_range(_mesh.active_not_local_elements_begin(), _mesh.active_not_local_elements_end()))
446 3500 : if (blockIDs().count(elem->subdomain_id()))
447 3570 : pull_requests[elem->processor_id()].push_back(elem->id());
448 : }
449 : else
450 : {
451 1419798 : for (const auto * const elem : _elements_to_push_pull)
452 673404 : pull_requests[elem->processor_id()].push_back(elem->id());
453 : }
454 :
455 : // First push
456 : {
457 : auto action_functor =
458 73591 : [this](const processor_id_type libmesh_dbg_var(pid), const std::vector<Datum> & sent_data)
459 : {
460 : mooseAssert(pid != this->processor_id(), "We do not send messages to ourself here");
461 747345 : for (const auto & pr : sent_data)
462 673754 : _a[pr.first] += pr.second;
463 146651 : };
464 73060 : TIMPI::push_parallel_vector_data(_communicator, push_data, action_functor);
465 : }
466 :
467 : // Then pull
468 : {
469 73626 : auto gather_functor = [this](const processor_id_type libmesh_dbg_var(pid),
470 : const std::vector<dof_id_type> & elem_ids,
471 : std::vector<VectorValue<ADReal>> & data_to_fill)
472 : {
473 : mooseAssert(pid != this->processor_id(), "We shouldn't be gathering from ourselves.");
474 73626 : data_to_fill.resize(elem_ids.size());
475 750530 : for (const auto i : index_range(elem_ids))
476 : {
477 676904 : const auto id = elem_ids[i];
478 676904 : auto it = _a.find(id);
479 : mooseAssert(it != _a.end(), "We should hold the value for this locally");
480 : data_to_fill[i] = it->second;
481 : }
482 146686 : };
483 :
484 73626 : auto action_functor = [this](const processor_id_type libmesh_dbg_var(pid),
485 : const std::vector<dof_id_type> & elem_ids,
486 : const std::vector<VectorValue<ADReal>> & filled_data)
487 : {
488 : mooseAssert(pid != this->processor_id(), "The request filler shouldn't have been ourselves");
489 : mooseAssert(elem_ids.size() == filled_data.size(), "I think these should be the same size");
490 750530 : for (const auto i : index_range(elem_ids))
491 676904 : _a[elem_ids[i]] = filled_data[i];
492 146686 : };
493 73060 : TIMPI::pull_parallel_vector_data(
494 73060 : _communicator, pull_requests, gather_functor, action_functor, &example);
495 : }
496 : }
497 :
498 : void
499 4049 : INSFVRhieChowInterpolator::ghostADataOnBoundary(const BoundaryID boundary_id)
500 : {
501 4049 : if (!needAComputation() || this->n_processors() == 1)
502 : return;
503 :
504 : // Ghost a for the elements on the boundary
505 17809 : for (auto elem_id : _moose_mesh.getBoundaryActiveSemiLocalElemIds(boundary_id))
506 : {
507 14895 : const auto & elem = _moose_mesh.elemPtr(elem_id);
508 : // no need to ghost if locally owned or far from local process
509 14895 : if (elem->processor_id() != this->processor_id() && elem->is_semilocal(this->processor_id()))
510 : // Adding to the a coefficient will make sure the final result gets communicated
511 439 : addToA(elem, 0, 0);
512 : }
513 :
514 : // Ghost a for the neighbors of the elements on the boundary
515 8263 : for (auto neighbor_id : _moose_mesh.getBoundaryActiveNeighborElemIds(boundary_id))
516 : {
517 5349 : const auto & neighbor = _moose_mesh.queryElemPtr(neighbor_id);
518 : // no need to ghost if locally owned or far from local process
519 7081 : if (neighbor->processor_id() != this->processor_id() &&
520 1732 : neighbor->is_semilocal(this->processor_id()))
521 : // Adding to the a coefficient will make sure the final result gets communicated
522 326 : addToA(neighbor, 0, 0);
523 : }
524 : }
525 :
526 : VectorValue<ADReal>
527 248933863 : INSFVRhieChowInterpolator::getVelocity(const Moose::FV::InterpMethod m,
528 : const FaceInfo & fi,
529 : const Moose::StateArg & time,
530 : const THREAD_ID tid,
531 : const bool subtract_mesh_velocity) const
532 : {
533 248933863 : const Elem * const elem = &fi.elem();
534 248933863 : const Elem * const neighbor = fi.neighborPtr();
535 248933863 : auto & vel = *_vel[tid];
536 248933863 : auto & p = *_ps[tid];
537 248933863 : auto * const u = _us[tid];
538 248933863 : MooseVariableFVReal * const v = _v ? _vs[tid] : nullptr;
539 248933863 : MooseVariableFVReal * const w = _w ? _ws[tid] : nullptr;
540 : // Check if skewness-correction is necessary
541 248933863 : const bool correct_skewness = velocitySkewCorrection(tid);
542 : auto incorporate_mesh_velocity =
543 248933863 : [this, tid, subtract_mesh_velocity, &time](const auto & space, auto & velocity)
544 : {
545 248933863 : if (_disps.size() && subtract_mesh_velocity)
546 4326420 : velocity -= _disps[tid]->dot(space, time);
547 497867726 : };
548 :
549 248933863 : if (Moose::FV::onBoundary(*this, fi))
550 : {
551 3886856 : const Elem * const boundary_elem = hasBlocks(elem->subdomain_id()) ? elem : neighbor;
552 3886856 : const Moose::FaceArg boundary_face{&fi,
553 : Moose::FV::LimiterType::CentralDifference,
554 : true,
555 : correct_skewness,
556 : boundary_elem,
557 3886856 : nullptr};
558 3886856 : auto velocity = vel(boundary_face, time);
559 3886856 : incorporate_mesh_velocity(boundary_face, velocity);
560 :
561 : // If not solving for velocity, clear derivatives
562 3886856 : if (_fe_problem.currentNlSysNum() != _momentum_sys_number)
563 2160 : return MetaPhysicL::raw_value(velocity);
564 : else
565 : return velocity;
566 : }
567 :
568 : VectorValue<ADReal> velocity;
569 :
570 245047007 : Moose::FaceArg face{
571 245047007 : &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr, nullptr};
572 : // Create the average face velocity (not corrected using RhieChow yet)
573 245047007 : velocity(0) = (*u)(face, time);
574 245047007 : if (v)
575 244447413 : velocity(1) = (*v)(face, time);
576 245047007 : if (w)
577 1223460 : velocity(2) = (*w)(face, time);
578 :
579 245047007 : incorporate_mesh_velocity(face, velocity);
580 :
581 : // If not solving for velocity, clear derivatives
582 245047007 : if (_fe_problem.currentNlSysNum() != _momentum_sys_number)
583 209520 : velocity = MetaPhysicL::raw_value(velocity);
584 :
585 : // Return if Rhie-Chow was not requested or if we have a porosity jump
586 245047007 : if (m == Moose::FV::InterpMethod::Average ||
587 245047007 : std::get<0>(NS::isPorosityJumpFace(epsilon(tid), fi, time)))
588 : return velocity;
589 :
590 : // Rhie-Chow coefficients are not available on initial
591 151246904 : if (_fe_problem.getCurrentExecuteOnFlag() == EXEC_INITIAL)
592 : {
593 0 : mooseDoOnce(mooseWarning("Cannot compute Rhie Chow coefficients on initial. Returning linearly "
594 : "interpolated velocities"););
595 : return velocity;
596 : }
597 151246904 : if (!_fe_problem.shouldSolve())
598 : {
599 0 : mooseDoOnce(mooseWarning("Cannot compute Rhie Chow coefficients if not solving. Returning "
600 : "linearly interpolated velocities"););
601 : return velocity;
602 : }
603 :
604 : mooseAssert(((m == Moose::FV::InterpMethod::RhieChow) &&
605 : (_velocity_interp_method == Moose::FV::InterpMethod::RhieChow)) ||
606 : _a_data_provided,
607 : "The 'a' coefficients have not been generated or provided for "
608 : "Rhie Chow velocity interpolation.");
609 :
610 : mooseAssert(neighbor && this->hasBlocks(neighbor->subdomain_id()),
611 : "We should be on an internal face...");
612 :
613 : // Get pressure gradient. This is the uncorrected gradient plus a correction from cell
614 : // centroid values on either side of the face
615 : const auto correct_skewness_p = pressureSkewCorrection(tid);
616 151246904 : const auto & grad_p = p.adGradSln(fi, time, correct_skewness_p);
617 :
618 : // Get uncorrected pressure gradient. This will use the element centroid gradient if we are
619 : // along a boundary face
620 151246904 : const auto & unc_grad_p = p.uncorrectedAdGradSln(fi, time, correct_skewness_p);
621 :
622 : // Volumetric Correction Method #1: pressure-based correction
623 : // Function that allows us to mark the face for which the Rhie-Chow interpolation is
624 : // inconsistent Normally, we should apply a reconstructed volume correction to the Rhie-Chow
625 : // coefficients However, since the fluxes at the face are given by the volume force we will
626 : // simply mark the face add the reverse pressure interpolation for these faces In brief, this
627 : // function is just marking the faces where the Rhie-Chow interpolation is inconsistent
628 : auto vf_indicator_pressure_based =
629 388880 : [this, &elem, &neighbor, &time, &fi, &correct_skewness](const Point & unit_basis_vector)
630 : {
631 : // Holders for the interpolated corrected and uncorrected volume force
632 : Real interp_vf;
633 : Real uncorrected_interp_vf;
634 :
635 : // Compute the corrected interpolated face value
636 388880 : Moose::FaceArg face{
637 388880 : &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr, nullptr};
638 :
639 : interp_vf = 0.0;
640 777760 : for (const auto i : make_range(_volumetric_force.size()))
641 388880 : interp_vf += (*this->_volumetric_force[i])(face, time);
642 :
643 : // Compute the uncorrected interpolated face value
644 : // For it to be consistent with the pressure gradient interpolation `uncorrectedAdGradSln`
645 : // the uncorrected volume force computation should follow the same Green-Gauss process
646 :
647 388880 : Real elem_value = 0.0;
648 388880 : Real neigh_value = 0.0;
649 :
650 : // Uncorrected interpolation - Step 1: loop over the faces of the element to compute
651 : // face-average cell value
652 : Real coord_multiplier;
653 388880 : const auto coord_type = _fe_problem.getCoordSystem(elem->subdomain_id());
654 : const unsigned int rz_radial_coord =
655 388880 : Moose::COORD_RZ ? _fe_problem.getAxisymmetricRadialCoord() : libMesh::invalid_uint;
656 :
657 1944400 : for (const auto side : make_range(elem->n_sides()))
658 : {
659 1555520 : const Elem * const loc_neighbor = elem->neighbor_ptr(side);
660 1555520 : const bool elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, loc_neighbor);
661 : const FaceInfo * const fi_loc =
662 2247728 : _moose_mesh.faceInfo(elem_has_fi ? elem : loc_neighbor,
663 692208 : elem_has_fi ? side : loc_neighbor->which_neighbor_am_i(elem));
664 :
665 1555520 : Moose::FaceArg loc_face{
666 1555520 : fi_loc, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, elem, nullptr};
667 :
668 1555520 : MooseMeshUtils::coordTransformFactor(
669 3111040 : elem->vertex_average(), coord_multiplier, coord_type, rz_radial_coord);
670 :
671 1555520 : Real face_volume_contribution = fi_loc->faceArea() *
672 1555520 : (neighbor->vertex_average() - elem->vertex_average()).norm() *
673 1555520 : coord_multiplier;
674 :
675 3111040 : for (const auto i : make_range(_volumetric_force.size()))
676 : {
677 : // Add which side (can be both, then we use a nullptr) of the face info the force is defined
678 : // on
679 1555520 : loc_face.face_side =
680 1555520 : this->_volumetric_force[i]->hasFaceSide(*fi_loc, true)
681 1555520 : ? (this->_volumetric_force[i]->hasFaceSide(*fi_loc, false) ? nullptr
682 : : fi_loc->elemPtr())
683 : : fi_loc->neighborPtr();
684 1555520 : elem_value += (*this->_volumetric_force[i])(loc_face, time) * face_volume_contribution *
685 : (fi_loc->normal() * unit_basis_vector);
686 : }
687 : }
688 388880 : elem_value = elem_value / elem->volume();
689 :
690 : // Uncorrected interpolation - Step 2: loop over the face of the neighbor to compute
691 : // face-average cell value
692 1944400 : for (const auto side : make_range(neighbor->n_sides()))
693 : {
694 1555520 : const Elem * const loc_elem = neighbor->neighbor_ptr(side);
695 1555520 : const bool elem_has_fi = Moose::FV::elemHasFaceInfo(*neighbor, loc_elem);
696 : const FaceInfo * const fi_loc =
697 2290928 : _moose_mesh.faceInfo(elem_has_fi ? neighbor : loc_elem,
698 735408 : elem_has_fi ? side : loc_elem->which_neighbor_am_i(neighbor));
699 :
700 1555520 : Moose::FaceArg loc_face{
701 1555520 : fi_loc, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, elem, nullptr};
702 :
703 1555520 : MooseMeshUtils::coordTransformFactor(
704 3111040 : neighbor->vertex_average(), coord_multiplier, coord_type, rz_radial_coord);
705 :
706 1555520 : Real face_volume_contribution = fi_loc->faceArea() *
707 1555520 : (elem->vertex_average() - neighbor->vertex_average()).norm() *
708 1555520 : coord_multiplier;
709 :
710 3111040 : for (const auto i : make_range(_volumetric_force.size()))
711 : {
712 1555520 : loc_face.face_side =
713 1555520 : this->_volumetric_force[i]->hasFaceSide(*fi_loc, true)
714 1555520 : ? (this->_volumetric_force[i]->hasFaceSide(*fi_loc, false) ? nullptr
715 : : fi_loc->elemPtr())
716 : : fi_loc->neighborPtr();
717 1555520 : neigh_value += (*this->_volumetric_force[i])(loc_face, time) * face_volume_contribution *
718 : (fi_loc->normal() * unit_basis_vector);
719 : }
720 : }
721 388880 : neigh_value = neigh_value / neighbor->volume();
722 :
723 : // Uncorrected interpolation - Step 3: interpolate element and neighbor reconstructed values
724 : // to the face
725 388880 : MooseMeshUtils::coordTransformFactor(
726 : fi.faceCentroid(), coord_multiplier, coord_type, rz_radial_coord);
727 388880 : interpolate(
728 : Moose::FV::InterpMethod::Average, uncorrected_interp_vf, elem_value, neigh_value, fi, true);
729 :
730 : // Return the flag indicator on which face the volume force correction is inconsistent
731 388880 : return MooseUtils::relativeFuzzyEqual(interp_vf, uncorrected_interp_vf, 1e-10) ? 0.0 : 1.0;
732 151246904 : };
733 :
734 : // Volumetric Correction Method #2: volume-based correction
735 : // In thery, pressure and velocity cannot be decoupled when a body force is present
736 : // Hence, we can de-activate the RC cofficient in faces that have a normal volume force
737 : // In the method we mark the faces with a non-zero volume force with recpect to the baseline
738 583320 : auto vf_indicator_force_based = [this, &time, &fi, &correct_skewness](Point & face_normal)
739 : {
740 : Real value = 0.0;
741 583320 : Moose::FaceArg loc_face{
742 583320 : &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr, nullptr};
743 :
744 1361080 : for (const auto i : make_range(_volumetric_force.size()))
745 777760 : value += (*_volumetric_force[i])(loc_face, time) * (face_normal * fi.normal());
746 583320 : if ((std::abs(value) - _baseline_volume_force) > 0)
747 : return 1.0;
748 : else
749 565548 : return 0.0;
750 151246904 : };
751 :
752 : const Point & elem_centroid = fi.elemCentroid();
753 : const Point & neighbor_centroid = fi.neighborCentroid();
754 151246904 : Real elem_volume = fi.elemVolume();
755 151246904 : Real neighbor_volume = fi.neighborVolume();
756 :
757 : // Now we need to perform the computations of D
758 151246904 : const auto elem_a = (*_a_read[tid])(makeElemArg(elem), time);
759 :
760 : mooseAssert(UserObject::_subproblem.getCoordSystem(elem->subdomain_id()) ==
761 : UserObject::_subproblem.getCoordSystem(neighbor->subdomain_id()),
762 : "Coordinate systems must be the same between the two elements");
763 :
764 : Real coord;
765 151246904 : coordTransformFactor(UserObject::_subproblem, elem->subdomain_id(), elem_centroid, coord);
766 :
767 151246904 : elem_volume *= coord;
768 :
769 151246904 : VectorValue<ADReal> elem_D = 0;
770 454060342 : for (const auto i : make_range(_dim))
771 : {
772 : mooseAssert(elem_a(i).value() != 0, "We should not be dividing by zero");
773 302813438 : elem_D(i) = elem_volume / elem_a(i);
774 : }
775 :
776 : VectorValue<ADReal> face_D;
777 :
778 151246904 : const auto neighbor_a = (*_a_read[tid])(makeElemArg(neighbor), time);
779 :
780 151246904 : coordTransformFactor(UserObject::_subproblem, neighbor->subdomain_id(), neighbor_centroid, coord);
781 151246904 : neighbor_volume *= coord;
782 :
783 151246904 : VectorValue<ADReal> neighbor_D = 0;
784 454060342 : for (const auto i : make_range(_dim))
785 : {
786 : mooseAssert(neighbor_a(i).value() != 0, "We should not be dividing by zero");
787 302813438 : neighbor_D(i) = neighbor_volume / neighbor_a(i);
788 : }
789 :
790 : // We require this to ensure that the correct interpolation weights are used.
791 : // This will change once the traditional weights are replaced by the weights
792 : // that are used by the skewness-correction.
793 151246904 : Moose::FV::InterpMethod coeff_interp_method = correct_skewness
794 151246904 : ? Moose::FV::InterpMethod::SkewCorrectedAverage
795 : : Moose::FV::InterpMethod::Average;
796 151246904 : Moose::FV::interpolate(coeff_interp_method, face_D, elem_D, neighbor_D, fi, true);
797 :
798 : // evaluate face porosity, see (18) in Hanimann 2021 or (11) in Nordlund 2016
799 151246904 : const auto face_eps = epsilon(tid)(face, time);
800 :
801 : // Perform the pressure correction. We don't use skewness-correction on the pressure since
802 : // it only influences the averaged cell gradients which cancel out in the correction
803 : // below.
804 454060342 : for (const auto i : make_range(_dim))
805 : {
806 : // "Standard" pressure-based RC interpolation
807 302813438 : velocity(i) -= face_D(i) * face_eps * (grad_p(i) - unc_grad_p(i));
808 :
809 302813438 : if (_bool_correct_vf)
810 : {
811 : // To solve the volume force incorrect interpolation, we add back the pressure gradient to the
812 : // RC-inconsistent faces regarding the marking method
813 : Point unit_basis_vector;
814 972200 : unit_basis_vector(i) = 1.0;
815 :
816 : // Get the value of the correction face indicator
817 : Real correction_indicator;
818 972200 : if (_volume_force_correction_method == "force-consistent")
819 583320 : correction_indicator = vf_indicator_force_based(unit_basis_vector);
820 : else
821 388880 : correction_indicator = vf_indicator_pressure_based(unit_basis_vector);
822 :
823 : // Correct back the velocity
824 972200 : velocity(i) += face_D(i) * face_eps * (grad_p(i) - unc_grad_p(i)) * correction_indicator;
825 : }
826 : }
827 :
828 : // If not solving for velocity, clear derivatives
829 151246904 : if (_fe_problem.currentNlSysNum() != _momentum_sys_number)
830 209520 : return MetaPhysicL::raw_value(velocity);
831 : else
832 : return velocity;
833 : }
|
