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 356838 : reynoldsPropertyDerivative(
47 : const Real & Re, const Real & rho, const Real & mu, const Real & drho, const Real & dmu)
48 : {
49 356838 : return Re * (drho / std::max(rho, 1e-6) - dmu / std::max(mu, 1e-8));
50 : }
51 :
52 : Real
53 356838 : 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 356838 : return (k * (mu * dcp + cp * dmu) - mu * cp * dk) / std::max(k * k, 1e-8);
61 : }
62 :
63 : template <typename T>
64 : T
65 1320472 : findUStar(const T & mu, const T & rho, const T & u, const Real dist)
66 : {
67 : // usually takes about 3-4 iterations
68 : constexpr int MAX_ITERS{50};
69 : constexpr Real REL_TOLERANCE{1e-6};
70 :
71 : // Check inputs
72 : mooseAssert(mu > 0, "Need a strictly positive viscosity");
73 : mooseAssert(rho > 0, "Need a strictly positive density");
74 : mooseAssert(u > 0, "Need a strictly positive velocity");
75 : mooseAssert(dist > 0, "Need a strictly positive wall distance");
76 :
77 790312 : const T nu = mu / rho;
78 :
79 : // Wall-function linearized guess
80 : const Real a_c = 1 / NS::von_karman_constant;
81 2910952 : const T b_c = 1.0 / NS::von_karman_constant * (std::log(NS::E_turb_constant * dist / mu) + 1.0);
82 : const T & c_c = u;
83 :
84 : /// This is important to reduce the number of nonlinear iterations
85 2910952 : T u_star = std::max(1e-20, (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c));
86 :
87 : // Newton-Raphson method to solve for u_star (friction velocity).
88 5223294 : for (int i = 0; i < MAX_ITERS; ++i)
89 : {
90 7291834 : T residual =
91 3154754 : u_star / NS::von_karman_constant * std::log(NS::E_turb_constant * u_star * dist / nu) - u;
92 7291834 : T deriv = (1.0 + std::log(NS::E_turb_constant * u_star * dist / nu)) / NS::von_karman_constant;
93 5223294 : T new_u_star = std::max(1e-20, u_star - residual / deriv);
94 :
95 : Real rel_err =
96 5223294 : std::abs(MetaPhysicL::raw_value(new_u_star - u_star) / MetaPhysicL::raw_value(new_u_star));
97 :
98 2068540 : u_star = new_u_star;
99 5223294 : if (rel_err < REL_TOLERANCE)
100 1320472 : return u_star;
101 : }
102 :
103 0 : mooseException("Could not find the wall friction velocity (mu: ",
104 : mu,
105 : " rho: ",
106 : rho,
107 : " velocity: ",
108 : u,
109 : " wall distance: ",
110 : dist,
111 : ")");
112 : }
113 : template Real findUStar<Real>(const Real & mu, const Real & rho, const Real & u, const Real dist);
114 : template ADReal
115 : findUStar<ADReal>(const ADReal & mu, const ADReal & rho, const ADReal & u, const Real dist);
116 :
117 : template <typename T>
118 : T
119 1206728 : findyPlus(const T & mu, const T & rho, const T & u, const Real dist)
120 : {
121 : // Fixed point iteration method to find y_plus
122 : // It should take 3 or 4 iterations
123 : constexpr int MAX_ITERS{10};
124 : constexpr Real REL_TOLERANCE{1e-2};
125 :
126 : // Check inputs
127 : mooseAssert(mu > 0, "Need a strictly positive viscosity");
128 : mooseAssert(u > 0, "Need a strictly positive velocity");
129 : mooseAssert(rho > 0, "Need a strictly positive density");
130 : mooseAssert(dist > 0, "Need a strictly positive wall distance");
131 :
132 : Real yPlusLast = 0.0;
133 1206728 : T yPlus = dist * u * rho / mu; // Assign initial value to laminar
134 1206728 : const Real rev_yPlusLam = 1.0 / MetaPhysicL::raw_value(yPlus);
135 1206728 : const T kappa_time_Re = NS::von_karman_constant * u * dist / (mu / rho);
136 : unsigned int iters = 0;
137 :
138 : do
139 : {
140 : yPlusLast = MetaPhysicL::raw_value(yPlus);
141 : // Negative y plus does not make sense
142 3622951 : yPlus = std::max(NS::min_y_plus, yPlus);
143 4381838 : yPlus = (kappa_time_Re + yPlus) / (1.0 + std::log(NS::E_turb_constant * yPlus));
144 3622951 : } while (std::abs(rev_yPlusLam * (MetaPhysicL::raw_value(yPlus) - yPlusLast)) > REL_TOLERANCE &&
145 : ++iters < MAX_ITERS);
146 :
147 1206728 : return std::max(NS::min_y_plus, yPlus);
148 : }
149 : template Real findyPlus<Real>(const Real & mu, const Real & rho, const Real & u, Real dist);
150 : template ADReal
151 : findyPlus<ADReal>(const ADReal & mu, const ADReal & rho, const ADReal & u, Real dist);
152 :
153 : template <typename T>
154 : T
155 81827794 : computeSpeed(const libMesh::VectorValue<T> & velocity)
156 : {
157 : // if the velocity is zero, then the norm function call fails because AD tries to calculate the
158 : // derivatives which causes a divide by zero - because d/dx(sqrt(f(x))) = 1/2/sqrt(f(x))*df/dx.
159 : // So add a bit of noise (based on hitchhiker's guide to the galaxy's meaning of life number) to
160 : // avoid this failure mode.
161 81827794 : return isZero(velocity) ? 1e-42 : velocity.norm();
162 : }
163 : template Real computeSpeed<Real>(const libMesh::VectorValue<Real> & velocity);
164 : template ADReal computeSpeed<ADReal>(const libMesh::VectorValue<ADReal> & velocity);
165 :
166 : template <typename T>
167 : T
168 904024 : computeShearStrainRateNormSquared(const Moose::Functor<T> & u,
169 : const Moose::Functor<T> * v,
170 : const Moose::Functor<T> * w,
171 : const Moose::ElemArg & elem_arg,
172 : const Moose::StateArg & state)
173 : {
174 904024 : const auto & grad_u = u.gradient(elem_arg, state);
175 904024 : const T Sij_xx = 2.0 * grad_u(0);
176 0 : T Sij_xy = 0.0;
177 0 : T Sij_xz = 0.0;
178 0 : T Sij_yy = 0.0;
179 0 : T Sij_yz = 0.0;
180 0 : T Sij_zz = 0.0;
181 :
182 0 : const T grad_xx = grad_u(0);
183 0 : T grad_xy = 0.0;
184 0 : T grad_xz = 0.0;
185 0 : T grad_yx = 0.0;
186 0 : T grad_yy = 0.0;
187 0 : T grad_yz = 0.0;
188 0 : T grad_zx = 0.0;
189 0 : T grad_zy = 0.0;
190 0 : T grad_zz = 0.0;
191 :
192 904024 : T trace = Sij_xx / 3.0;
193 :
194 904024 : if (v) // dim >= 2
195 : {
196 904024 : const auto & grad_v = (*v).gradient(elem_arg, state);
197 904024 : Sij_xy = grad_u(1) + grad_v(0);
198 904024 : Sij_yy = 2.0 * grad_v(1);
199 :
200 0 : grad_xy = grad_u(1);
201 0 : grad_yx = grad_v(0);
202 0 : grad_yy = grad_v(1);
203 :
204 904024 : trace += Sij_yy / 3.0;
205 :
206 904024 : if (w) // dim >= 3
207 : {
208 0 : const auto & grad_w = (*w).gradient(elem_arg, state);
209 :
210 0 : Sij_xz = grad_u(2) + grad_w(0);
211 0 : Sij_yz = grad_v(2) + grad_w(1);
212 0 : Sij_zz = 2.0 * grad_w(2);
213 :
214 0 : grad_xz = grad_u(2);
215 0 : grad_yz = grad_v(2);
216 0 : grad_zx = grad_w(0);
217 0 : grad_zy = grad_w(1);
218 0 : grad_zz = grad_w(2);
219 :
220 0 : trace += Sij_zz / 3.0;
221 : }
222 : }
223 :
224 904024 : return (Sij_xx - trace) * grad_xx + Sij_xy * grad_xy + Sij_xz * grad_xz + Sij_xy * grad_yx +
225 904024 : (Sij_yy - trace) * grad_yy + Sij_yz * grad_yz + Sij_xz * grad_zx + Sij_yz * grad_zy +
226 904024 : (Sij_zz - trace) * grad_zz;
227 : }
228 : template Real computeShearStrainRateNormSquared<Real>(const Moose::Functor<Real> & u,
229 : const Moose::Functor<Real> * v,
230 : const Moose::Functor<Real> * w,
231 : const Moose::ElemArg & elem_arg,
232 : const Moose::StateArg & state);
233 : template ADReal computeShearStrainRateNormSquared<ADReal>(const Moose::Functor<ADReal> & u,
234 : const Moose::Functor<ADReal> * v,
235 : const Moose::Functor<ADReal> * w,
236 : const Moose::ElemArg & elem_arg,
237 : const Moose::StateArg & state);
238 :
239 : /// Bounded element maps for wall treatment
240 : void
241 18621 : getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_names,
242 : const FEProblemBase & fe_problem,
243 : const SubProblem & subproblem,
244 : const std::set<SubdomainID> & block_ids,
245 : std::unordered_set<const Elem *> & wall_bounded)
246 : {
247 :
248 : wall_bounded.clear();
249 18621 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_names);
250 :
251 : // We define these lambdas so that we can fetch the bounded elements from other
252 : // processors.
253 14120 : auto gather_functor = [&subproblem, &wall_bounded](const processor_id_type libmesh_dbg_var(pid),
254 : const std::vector<dof_id_type> & elem_ids,
255 : std::vector<unsigned char> & data_to_fill)
256 : {
257 : mooseAssert(pid != subproblem.processor_id(), "We shouldn't be gathering from ourselves.");
258 14120 : data_to_fill.resize(elem_ids.size());
259 :
260 14120 : const auto & mesh = subproblem.mesh().getMesh();
261 :
262 172472 : for (const auto i : index_range(elem_ids))
263 : {
264 158352 : const auto elem = mesh.elem_ptr(elem_ids[i]);
265 316704 : data_to_fill[i] = wall_bounded.count(elem) != 0;
266 : }
267 32741 : };
268 :
269 14120 : auto action_functor = [&subproblem, &wall_bounded](const processor_id_type libmesh_dbg_var(pid),
270 : const std::vector<dof_id_type> & elem_ids,
271 : const std::vector<unsigned char> & filled_data)
272 : {
273 : mooseAssert(pid != subproblem.processor_id(),
274 : "The request filler shouldn't have been ourselves");
275 : mooseAssert(elem_ids.size() == filled_data.size(), "I think these should be the same size");
276 :
277 14120 : const auto & mesh = subproblem.mesh().getMesh();
278 :
279 172472 : for (const auto i : index_range(elem_ids))
280 : {
281 158352 : const auto elem = mesh.elem_ptr(elem_ids[i]);
282 158352 : if (filled_data[i])
283 14432 : wall_bounded.insert(elem);
284 : }
285 32741 : };
286 :
287 : // We need these elements from other processors
288 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> elem_ids_requested;
289 :
290 4727410 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
291 2345084 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
292 11695180 : for (const auto i_side : elem->side_index_range())
293 : {
294 : // This is needed because in some cases the internal boundary is registered
295 : // to the neighbor element
296 : std::set<BoundaryID> combined_side_bds;
297 9356144 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
298 9356144 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
299 9356144 : if (const auto neighbor = elem->neighbor_ptr(i_side))
300 : {
301 8752224 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
302 8752224 : const auto & neighbor_bnds = subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
303 8752224 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
304 :
305 : // If the neighbor lives on the first layer of the ghost region then we would
306 : // like to grab its value as well (if it exists)
307 8752224 : if (neighbor->processor_id() != subproblem.processor_id() &&
308 : block_ids.find(neighbor->subdomain_id()) != block_ids.end())
309 158352 : elem_ids_requested[neighbor->processor_id()].push_back(neighbor->id());
310 8752224 : }
311 :
312 25680252 : for (const auto wall_id : wall_boundary_ids)
313 : if (combined_side_bds.count(wall_id))
314 : {
315 : wall_bounded.insert(elem);
316 : break;
317 : }
318 9374765 : }
319 :
320 : unsigned char * bool_ex = nullptr;
321 18621 : TIMPI::pull_parallel_vector_data(
322 : subproblem.comm(), elem_ids_requested, gather_functor, action_functor, bool_ex);
323 18621 : }
324 :
325 : /// Bounded element face distances for wall treatment
326 : void
327 5049 : getWallDistance(const std::vector<BoundaryName> & wall_boundary_name,
328 : const FEProblemBase & fe_problem,
329 : const SubProblem & subproblem,
330 : const std::set<SubdomainID> & block_ids,
331 : std::map<const Elem *, std::vector<Real>> & dist_map)
332 : {
333 : dist_map.clear();
334 5049 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_name);
335 :
336 1321018 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
337 655460 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
338 3268660 : for (const auto i_side : elem->side_index_range())
339 : {
340 : // This is needed because in some cases the internal boundary is registered
341 : // to the neighbor element
342 : std::set<BoundaryID> combined_side_bds;
343 2614928 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
344 2614928 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
345 2614928 : if (const auto neighbor = elem->neighbor_ptr(i_side))
346 : {
347 2449272 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
348 : const std::vector<BoundaryID> & neighbor_bnds =
349 2449272 : subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
350 2449272 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
351 2449272 : }
352 :
353 10016208 : for (const auto wall_id : wall_boundary_ids)
354 : if (combined_side_bds.count(wall_id))
355 : {
356 : // The list below stores the face infos with respect to their owning elements,
357 : // depending on the block restriction we might encounter situations where the
358 : // element outside of the block owns the face info.
359 138016 : const auto & neighbor = elem->neighbor_ptr(i_side);
360 138016 : const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
361 138016 : const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
362 288 : const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
363 :
364 138016 : const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
365 138016 : const auto & elem_centroid = elem_has_fi ? fi->elemCentroid() : fi->neighborCentroid();
366 138016 : const Real dist = std::abs((elem_centroid - fi->faceCentroid()) * fi->normal());
367 138016 : dist_map[elem].push_back(dist);
368 : }
369 2619977 : }
370 5049 : }
371 :
372 : /// Face arguments to wall-bounded faces for wall treatment
373 : void
374 5049 : getElementFaceArgs(const std::vector<BoundaryName> & wall_boundary_name,
375 : const FEProblemBase & fe_problem,
376 : const SubProblem & subproblem,
377 : const std::set<SubdomainID> & block_ids,
378 : std::map<const Elem *, std::vector<const FaceInfo *>> & face_info_map)
379 : {
380 : face_info_map.clear();
381 5049 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_name);
382 :
383 1321018 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
384 655460 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
385 3268660 : for (const auto i_side : elem->side_index_range())
386 : {
387 : // This is needed because in some cases the internal boundary is registered
388 : // to the neighbor element
389 : std::set<BoundaryID> combined_side_bds;
390 2614928 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
391 2614928 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
392 2614928 : if (elem->neighbor_ptr(i_side) && !elem->neighbor_ptr(i_side)->is_remote())
393 : {
394 2449272 : const auto neighbor = elem->neighbor_ptr(i_side);
395 2449272 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
396 : const std::vector<BoundaryID> & neighbor_bnds =
397 2449272 : subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
398 2449272 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
399 2449272 : }
400 :
401 10016208 : for (const auto wall_id : wall_boundary_ids)
402 : if (combined_side_bds.count(wall_id))
403 : {
404 : // The list below stores the face infos with respect to their owning elements,
405 : // depending on the block restriction we might encounter situations where the
406 : // element outside of the block owns the face info.
407 138016 : const auto & neighbor = elem->neighbor_ptr(i_side);
408 138016 : const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
409 138016 : const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
410 288 : const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
411 :
412 138016 : const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
413 138016 : face_info_map[elem].push_back(fi);
414 : }
415 2619977 : }
416 5049 : }
417 : }
|
