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