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 "NavierStokesMethods.h"
11 : #include "MooseError.h"
12 : #include "libmesh/vector_value.h"
13 : #include "NS.h"
14 :
15 : namespace NS
16 : {
17 : int
18 0 : delta(unsigned int i, unsigned int j)
19 : {
20 0 : if (i == j)
21 : return 1;
22 : else
23 0 : return 0;
24 : }
25 :
26 : int
27 0 : computeSign(const Real & a)
28 : {
29 0 : return a > 0 ? 1 : (a < 0 ? -1 : 0);
30 : }
31 :
32 : unsigned int
33 0 : getIndex(const Real & p, const std::vector<Real> & bounds)
34 : {
35 0 : if (p < bounds.front() || p > bounds.back())
36 0 : mooseError("Point exceeds bounds of domain!");
37 :
38 0 : for (unsigned int i = 1; i < bounds.size(); ++i)
39 0 : if (p <= bounds[i])
40 0 : return i - 1;
41 :
42 0 : return bounds.size() - 2;
43 : }
44 :
45 : Real
46 239094 : reynoldsPropertyDerivative(
47 : const Real & Re, const Real & rho, const Real & mu, const Real & drho, const Real & dmu)
48 : {
49 239094 : return Re * (drho / std::max(rho, 1e-6) - dmu / std::max(mu, 1e-8));
50 : }
51 :
52 : Real
53 239094 : prandtlPropertyDerivative(const Real & mu,
54 : const Real & cp,
55 : const Real & k,
56 : const Real & dmu,
57 : const Real & dcp,
58 : const Real & dk)
59 : {
60 239094 : return (k * (mu * dcp + cp * dmu) - mu * cp * dk) / std::max(k * k, 1e-8);
61 : }
62 :
63 : template <typename T>
64 : T
65 1394376 : findUStar(const T & mu, const T & rho, const T & u, const Real dist)
66 : {
67 : using std::log, std::sqrt, std::pow, std::max, std::abs;
68 :
69 : // usually takes about 3-4 iterations
70 : constexpr int MAX_ITERS{50};
71 : constexpr Real REL_TOLERANCE{1e-6};
72 :
73 : // Check inputs
74 : mooseAssert(mu > 0, "Need a strictly positive viscosity");
75 : mooseAssert(rho > 0, "Need a strictly positive density");
76 : mooseAssert(u > 0, "Need a strictly positive velocity");
77 : mooseAssert(dist > 0, "Need a strictly positive wall distance");
78 :
79 959008 : const T nu = mu / rho;
80 :
81 : // Wall-function linearized guess
82 : const Real a_c = 1 / NS::von_karman_constant;
83 2700480 : const T b_c = 1.0 / NS::von_karman_constant * (log(NS::E_turb_constant * dist / mu) + 1.0);
84 : const T & c_c = u;
85 :
86 : /// This is important to reduce the number of nonlinear iterations
87 2700480 : T u_star = max(1e-20, (-b_c + sqrt(pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c));
88 :
89 : // Newton-Raphson method to solve for u_star (friction velocity).
90 5531276 : for (int i = 0; i < MAX_ITERS; ++i)
91 : {
92 7231916 : T residual =
93 3830636 : u_star / NS::von_karman_constant * log(NS::E_turb_constant * u_star * dist / nu) - u;
94 7231916 : T deriv = (1.0 + log(NS::E_turb_constant * u_star * dist / nu)) / NS::von_karman_constant;
95 5531276 : T new_u_star = max(1e-20, u_star - residual / deriv);
96 :
97 : Real rel_err =
98 5531276 : abs(MetaPhysicL::raw_value(new_u_star - u_star) / MetaPhysicL::raw_value(new_u_star));
99 :
100 1700640 : u_star = new_u_star;
101 5531276 : if (rel_err < REL_TOLERANCE)
102 1394376 : return u_star;
103 : }
104 :
105 0 : mooseException("Could not find the wall friction velocity (mu: ",
106 : mu,
107 : " rho: ",
108 : rho,
109 : " velocity: ",
110 : u,
111 : " wall distance: ",
112 : dist,
113 : ")");
114 : }
115 : template Real findUStar<Real>(const Real & mu, const Real & rho, const Real & u, const Real dist);
116 : template ADReal
117 : findUStar<ADReal>(const ADReal & mu, const ADReal & rho, const ADReal & u, const Real dist);
118 :
119 : template <typename T>
120 : T
121 1332116 : findyPlus(const T & mu, const T & rho, const T & u, const Real dist)
122 : {
123 : using std::max, std::log, std::abs;
124 :
125 : // Fixed point iteration method to find y_plus
126 : // It should take 3 or 4 iterations
127 : constexpr int MAX_ITERS{10};
128 : constexpr Real REL_TOLERANCE{1e-2};
129 :
130 : // Check inputs
131 : mooseAssert(mu > 0, "Need a strictly positive viscosity");
132 : mooseAssert(u > 0, "Need a strictly positive velocity");
133 : mooseAssert(rho > 0, "Need a strictly positive density");
134 : mooseAssert(dist > 0, "Need a strictly positive wall distance");
135 :
136 : Real yPlusLast = 0.0;
137 1332116 : T yPlus = dist * u * rho / mu; // Assign initial value to laminar
138 1332116 : const Real rev_yPlusLam = 1.0 / MetaPhysicL::raw_value(yPlus);
139 1332116 : const T kappa_time_Re = NS::von_karman_constant * u * dist / (mu / rho);
140 : unsigned int iters = 0;
141 :
142 : do
143 : {
144 : yPlusLast = MetaPhysicL::raw_value(yPlus);
145 : // Negative y plus does not make sense
146 3998198 : yPlus = max(NS::min_y_plus, yPlus);
147 4564684 : yPlus = (kappa_time_Re + yPlus) / (1.0 + log(NS::E_turb_constant * yPlus));
148 3998198 : } while (abs(rev_yPlusLam * (MetaPhysicL::raw_value(yPlus) - yPlusLast)) > REL_TOLERANCE &&
149 : ++iters < MAX_ITERS);
150 :
151 1332116 : return max(NS::min_y_plus, yPlus);
152 : }
153 :
154 : template Real findyPlus<Real>(const Real & mu, const Real & rho, const Real & u, Real dist);
155 : template ADReal
156 : findyPlus<ADReal>(const ADReal & mu, const ADReal & rho, const ADReal & u, Real dist);
157 :
158 : template <typename T>
159 : T
160 57912582 : computeSpeed(const libMesh::VectorValue<T> & velocity)
161 : {
162 : // if the velocity is zero, then the norm function call fails because AD tries to calculate the
163 : // derivatives which causes a divide by zero - because d/dx(sqrt(f(x))) = 1/2/sqrt(f(x))*df/dx.
164 : // So add a bit of noise (based on hitchhiker's guide to the galaxy's meaning of life number) to
165 : // avoid this failure mode.
166 57912582 : return isZero(velocity) ? 1e-42 : velocity.norm();
167 : }
168 : template Real computeSpeed<Real>(const libMesh::VectorValue<Real> & velocity);
169 : template ADReal computeSpeed<ADReal>(const libMesh::VectorValue<ADReal> & velocity);
170 :
171 : template <typename T>
172 : T
173 1474324 : computeShearStrainRateNormSquared(const Moose::Functor<T> & u,
174 : const Moose::Functor<T> * v,
175 : const Moose::Functor<T> * w,
176 : const Moose::ElemArg & elem_arg,
177 : const Moose::StateArg & state,
178 : const Moose::CoordinateSystemType coord_sys,
179 : const unsigned int rz_radial_coord)
180 : {
181 1474324 : const auto & grad_u = u.gradient(elem_arg, state);
182 1474324 : const T Sij_xx = 2.0 * grad_u(0);
183 536308 : T Sij_xy = 0.0;
184 536308 : T Sij_xz = 0.0;
185 536308 : T Sij_yy = 0.0;
186 536308 : T Sij_yz = 0.0;
187 536308 : T Sij_zz = 0.0;
188 :
189 536308 : const T grad_xx = grad_u(0);
190 536308 : T grad_xy = 0.0;
191 536308 : T grad_xz = 0.0;
192 536308 : T grad_yx = 0.0;
193 536308 : T grad_yy = 0.0;
194 536308 : T grad_yz = 0.0;
195 536308 : T grad_zx = 0.0;
196 536308 : T grad_zy = 0.0;
197 536308 : T grad_zz = 0.0;
198 :
199 1474324 : T trace = Sij_xx / 3.0;
200 :
201 1474324 : if (v) // dim >= 2
202 : {
203 1474324 : const auto & grad_v = (*v).gradient(elem_arg, state);
204 1474324 : Sij_xy = grad_u(1) + grad_v(0);
205 2010632 : Sij_yy = 2.0 * grad_v(1);
206 :
207 536308 : grad_xy = grad_u(1);
208 536308 : grad_yx = grad_v(0);
209 536308 : grad_yy = grad_v(1);
210 :
211 2010632 : trace += Sij_yy / 3.0;
212 :
213 1474324 : if (w) // dim >= 3
214 : {
215 0 : const auto & grad_w = (*w).gradient(elem_arg, state);
216 :
217 0 : Sij_xz = grad_u(2) + grad_w(0);
218 0 : Sij_yz = grad_v(2) + grad_w(1);
219 0 : Sij_zz = 2.0 * grad_w(2);
220 :
221 0 : grad_xz = grad_u(2);
222 0 : grad_yz = grad_v(2);
223 0 : grad_zx = grad_w(0);
224 0 : grad_zy = grad_w(1);
225 0 : grad_zz = grad_w(2);
226 :
227 0 : trace += Sij_zz / 3.0;
228 : }
229 : }
230 :
231 1474324 : if (coord_sys == Moose::COORD_RZ)
232 : {
233 : mooseAssert(elem_arg.elem, "ElemArg must reference an element for cylindrical calculations.");
234 0 : const auto radius = elem_arg.elem->vertex_average()(rz_radial_coord);
235 : mooseAssert(radius > 0.0, "Axisymmetric radial coordinate should be positive.");
236 :
237 : const Moose::Functor<T> * radial_functor = nullptr;
238 0 : switch (rz_radial_coord)
239 : {
240 : case 0:
241 : radial_functor = &u;
242 : break;
243 0 : case 1:
244 : radial_functor = v;
245 0 : break;
246 0 : default:
247 0 : mooseError("Unsupported axisymmetric radial coordinate index: ", rz_radial_coord);
248 : }
249 :
250 : mooseAssert(radial_functor,
251 : "The functor corresponding to the axisymmetric radial velocity must be provided.");
252 0 : const T radial_velocity = (*radial_functor)(elem_arg, state);
253 0 : Sij_zz = 2.0 * radial_velocity / radius;
254 0 : grad_zz = radial_velocity / radius;
255 0 : trace += Sij_zz / 3.0;
256 : }
257 :
258 : // Note that in RZ we can only do X or Y axis. However, due to the hoop stress the
259 : // Z direction (polar coordinate) tensors are not 0 in this case.
260 938016 : return (Sij_xx - trace) * grad_xx + Sij_xy * grad_xy + Sij_xz * grad_xz + Sij_xy * grad_yx +
261 938016 : (Sij_yy - trace) * grad_yy + Sij_yz * grad_yz + Sij_xz * grad_zx + Sij_yz * grad_zy +
262 1474324 : (Sij_zz - trace) * grad_zz;
263 : }
264 : template Real computeShearStrainRateNormSquared<Real>(const Moose::Functor<Real> & u,
265 : const Moose::Functor<Real> * v,
266 : const Moose::Functor<Real> * w,
267 : const Moose::ElemArg & elem_arg,
268 : const Moose::StateArg & state,
269 : const Moose::CoordinateSystemType coord_sys,
270 : const unsigned int rz_radial_coord);
271 : template ADReal
272 : computeShearStrainRateNormSquared<ADReal>(const Moose::Functor<ADReal> & u,
273 : const Moose::Functor<ADReal> * v,
274 : const Moose::Functor<ADReal> * w,
275 : const Moose::ElemArg & elem_arg,
276 : const Moose::StateArg & state,
277 : const Moose::CoordinateSystemType coord_sys,
278 : const unsigned int rz_radial_coord);
279 :
280 : /// Bounded element maps for wall treatment
281 : void
282 10417 : getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_names,
283 : const FEProblemBase & fe_problem,
284 : const SubProblem & subproblem,
285 : const std::set<SubdomainID> & block_ids,
286 : std::unordered_set<const Elem *> & wall_bounded)
287 : {
288 :
289 : wall_bounded.clear();
290 10417 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_names);
291 :
292 : // We define these lambdas so that we can fetch the bounded elements from other
293 : // processors.
294 5504 : auto gather_functor = [&subproblem, &wall_bounded](const processor_id_type libmesh_dbg_var(pid),
295 : const std::vector<dof_id_type> & elem_ids,
296 : std::vector<unsigned char> & data_to_fill)
297 : {
298 : mooseAssert(pid != subproblem.processor_id(), "We shouldn't be gathering from ourselves.");
299 5504 : data_to_fill.resize(elem_ids.size());
300 :
301 5504 : const auto & mesh = subproblem.mesh().getMesh();
302 :
303 58204 : for (const auto i : index_range(elem_ids))
304 : {
305 52700 : const auto elem = mesh.elem_ptr(elem_ids[i]);
306 105400 : data_to_fill[i] = wall_bounded.count(elem) != 0;
307 : }
308 15921 : };
309 :
310 5504 : auto action_functor = [&subproblem, &wall_bounded](const processor_id_type libmesh_dbg_var(pid),
311 : const std::vector<dof_id_type> & elem_ids,
312 : const std::vector<unsigned char> & filled_data)
313 : {
314 : mooseAssert(pid != subproblem.processor_id(),
315 : "The request filler shouldn't have been ourselves");
316 : mooseAssert(elem_ids.size() == filled_data.size(), "I think these should be the same size");
317 :
318 5504 : const auto & mesh = subproblem.mesh().getMesh();
319 :
320 58204 : for (const auto i : index_range(elem_ids))
321 : {
322 52700 : const auto elem = mesh.elem_ptr(elem_ids[i]);
323 52700 : if (filled_data[i])
324 4804 : wall_bounded.insert(elem);
325 : }
326 15921 : };
327 :
328 : // We need these elements from other processors
329 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> elem_ids_requested;
330 :
331 2798538 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
332 1388852 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
333 6924100 : for (const auto i_side : elem->side_index_range())
334 : {
335 : // This is needed because in some cases the internal boundary is registered
336 : // to the neighbor element
337 : std::set<BoundaryID> combined_side_bds;
338 5539280 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
339 5539280 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
340 5539280 : if (const auto neighbor = elem->neighbor_ptr(i_side))
341 : {
342 5166296 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
343 5166296 : const auto & neighbor_bnds = subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
344 5166296 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
345 :
346 : // If the neighbor lives on the first layer of the ghost region then we would
347 : // like to grab its value as well (if it exists)
348 5166296 : if (neighbor->processor_id() != subproblem.processor_id() &&
349 : block_ids.find(neighbor->subdomain_id()) != block_ids.end())
350 52700 : elem_ids_requested[neighbor->processor_id()].push_back(neighbor->id());
351 5166296 : }
352 :
353 15151804 : for (const auto wall_id : wall_boundary_ids)
354 : if (combined_side_bds.count(wall_id))
355 : {
356 : wall_bounded.insert(elem);
357 : break;
358 : }
359 5549697 : }
360 :
361 : unsigned char * bool_ex = nullptr;
362 10417 : TIMPI::pull_parallel_vector_data(
363 : subproblem.comm(), elem_ids_requested, gather_functor, action_functor, bool_ex);
364 10417 : }
365 :
366 : /// Bounded element face distances for wall treatment
367 : void
368 2941 : getWallDistance(const std::vector<BoundaryName> & wall_boundary_name,
369 : const FEProblemBase & fe_problem,
370 : const SubProblem & subproblem,
371 : const std::set<SubdomainID> & block_ids,
372 : std::map<const Elem *, std::vector<Real>> & dist_map)
373 : {
374 : dist_map.clear();
375 2941 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_name);
376 :
377 786290 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
378 390204 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
379 1945260 : for (const auto i_side : elem->side_index_range())
380 : {
381 : // This is needed because in some cases the internal boundary is registered
382 : // to the neighbor element
383 : std::set<BoundaryID> combined_side_bds;
384 1556208 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
385 1556208 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
386 1556208 : if (const auto neighbor = elem->neighbor_ptr(i_side))
387 : {
388 1452192 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
389 : const std::vector<BoundaryID> & neighbor_bnds =
390 1452192 : subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
391 1452192 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
392 1452192 : }
393 :
394 5927344 : for (const auto wall_id : wall_boundary_ids)
395 : if (combined_side_bds.count(wall_id))
396 : {
397 : // The list below stores the face infos with respect to their owning elements,
398 : // depending on the block restriction we might encounter situations where the
399 : // element outside of the block owns the face info.
400 82628 : const auto & neighbor = elem->neighbor_ptr(i_side);
401 82628 : const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
402 82628 : const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
403 192 : const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
404 :
405 82628 : const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
406 82628 : const auto & elem_centroid = elem_has_fi ? fi->elemCentroid() : fi->neighborCentroid();
407 82628 : const Real dist = std::abs((elem_centroid - fi->faceCentroid()) * fi->normal());
408 82628 : dist_map[elem].push_back(dist);
409 : }
410 1559149 : }
411 2941 : }
412 :
413 : /// Face arguments to wall-bounded faces for wall treatment
414 : void
415 2941 : getElementFaceArgs(const std::vector<BoundaryName> & wall_boundary_name,
416 : const FEProblemBase & fe_problem,
417 : const SubProblem & subproblem,
418 : const std::set<SubdomainID> & block_ids,
419 : std::map<const Elem *, std::vector<const FaceInfo *>> & face_info_map)
420 : {
421 : face_info_map.clear();
422 2941 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_name);
423 :
424 786290 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
425 390204 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
426 1945260 : for (const auto i_side : elem->side_index_range())
427 : {
428 : // This is needed because in some cases the internal boundary is registered
429 : // to the neighbor element
430 : std::set<BoundaryID> combined_side_bds;
431 1556208 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
432 1556208 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
433 1556208 : if (elem->neighbor_ptr(i_side) && !elem->neighbor_ptr(i_side)->is_remote())
434 : {
435 1452192 : const auto neighbor = elem->neighbor_ptr(i_side);
436 1452192 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
437 : const std::vector<BoundaryID> & neighbor_bnds =
438 1452192 : subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
439 1452192 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
440 1452192 : }
441 :
442 5927344 : for (const auto wall_id : wall_boundary_ids)
443 : if (combined_side_bds.count(wall_id))
444 : {
445 : // The list below stores the face infos with respect to their owning elements,
446 : // depending on the block restriction we might encounter situations where the
447 : // element outside of the block owns the face info.
448 82628 : const auto & neighbor = elem->neighbor_ptr(i_side);
449 82628 : const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
450 82628 : const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
451 192 : const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
452 :
453 82628 : const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
454 82628 : face_info_map[elem].push_back(fi);
455 : }
456 1559149 : }
457 2941 : }
458 : }
|
