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 : 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 790312 : const T nu = mu / rho;
80 :
81 : // Wall-function linearized guess
82 : const Real a_c = 1 / NS::von_karman_constant;
83 2910952 : 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 2910952 : 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 5223294 : for (int i = 0; i < MAX_ITERS; ++i)
91 : {
92 7291834 : T residual =
93 3154754 : u_star / NS::von_karman_constant * log(NS::E_turb_constant * u_star * dist / nu) - u;
94 7291834 : T deriv = (1.0 + log(NS::E_turb_constant * u_star * dist / nu)) / NS::von_karman_constant;
95 5223294 : T new_u_star = max(1e-20, u_star - residual / deriv);
96 :
97 : Real rel_err =
98 5223294 : abs(MetaPhysicL::raw_value(new_u_star - u_star) / MetaPhysicL::raw_value(new_u_star));
99 :
100 2068540 : u_star = new_u_star;
101 5223294 : if (rel_err < REL_TOLERANCE)
102 1320472 : 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 1206728 : 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 1206728 : T yPlus = dist * u * rho / mu; // Assign initial value to laminar
138 1206728 : const Real rev_yPlusLam = 1.0 / MetaPhysicL::raw_value(yPlus);
139 1206728 : 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 3622951 : yPlus = max(NS::min_y_plus, yPlus);
147 4381838 : yPlus = (kappa_time_Re + yPlus) / (1.0 + log(NS::E_turb_constant * yPlus));
148 3622951 : } while (abs(rev_yPlusLam * (MetaPhysicL::raw_value(yPlus) - yPlusLast)) > REL_TOLERANCE &&
149 : ++iters < MAX_ITERS);
150 :
151 1206728 : 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 82916346 : 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 82916346 : 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 904024 : 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 : {
179 904024 : const auto & grad_u = u.gradient(elem_arg, state);
180 904024 : const T Sij_xx = 2.0 * grad_u(0);
181 0 : T Sij_xy = 0.0;
182 0 : T Sij_xz = 0.0;
183 0 : T Sij_yy = 0.0;
184 0 : T Sij_yz = 0.0;
185 0 : T Sij_zz = 0.0;
186 :
187 0 : const T grad_xx = grad_u(0);
188 0 : T grad_xy = 0.0;
189 0 : T grad_xz = 0.0;
190 0 : T grad_yx = 0.0;
191 0 : T grad_yy = 0.0;
192 0 : T grad_yz = 0.0;
193 0 : T grad_zx = 0.0;
194 0 : T grad_zy = 0.0;
195 0 : T grad_zz = 0.0;
196 :
197 904024 : T trace = Sij_xx / 3.0;
198 :
199 904024 : if (v) // dim >= 2
200 : {
201 904024 : const auto & grad_v = (*v).gradient(elem_arg, state);
202 904024 : Sij_xy = grad_u(1) + grad_v(0);
203 904024 : Sij_yy = 2.0 * grad_v(1);
204 :
205 0 : grad_xy = grad_u(1);
206 0 : grad_yx = grad_v(0);
207 0 : grad_yy = grad_v(1);
208 :
209 904024 : trace += Sij_yy / 3.0;
210 :
211 904024 : if (w) // dim >= 3
212 : {
213 0 : const auto & grad_w = (*w).gradient(elem_arg, state);
214 :
215 0 : Sij_xz = grad_u(2) + grad_w(0);
216 0 : Sij_yz = grad_v(2) + grad_w(1);
217 0 : Sij_zz = 2.0 * grad_w(2);
218 :
219 0 : grad_xz = grad_u(2);
220 0 : grad_yz = grad_v(2);
221 0 : grad_zx = grad_w(0);
222 0 : grad_zy = grad_w(1);
223 0 : grad_zz = grad_w(2);
224 :
225 0 : trace += Sij_zz / 3.0;
226 : }
227 : }
228 :
229 904024 : return (Sij_xx - trace) * grad_xx + Sij_xy * grad_xy + Sij_xz * grad_xz + Sij_xy * grad_yx +
230 904024 : (Sij_yy - trace) * grad_yy + Sij_yz * grad_yz + Sij_xz * grad_zx + Sij_yz * grad_zy +
231 904024 : (Sij_zz - trace) * grad_zz;
232 : }
233 : template Real computeShearStrainRateNormSquared<Real>(const Moose::Functor<Real> & u,
234 : const Moose::Functor<Real> * v,
235 : const Moose::Functor<Real> * w,
236 : const Moose::ElemArg & elem_arg,
237 : const Moose::StateArg & state);
238 : template ADReal computeShearStrainRateNormSquared<ADReal>(const Moose::Functor<ADReal> & u,
239 : const Moose::Functor<ADReal> * v,
240 : const Moose::Functor<ADReal> * w,
241 : const Moose::ElemArg & elem_arg,
242 : const Moose::StateArg & state);
243 :
244 : /// Bounded element maps for wall treatment
245 : void
246 18621 : getWallBoundedElements(const std::vector<BoundaryName> & wall_boundary_names,
247 : const FEProblemBase & fe_problem,
248 : const SubProblem & subproblem,
249 : const std::set<SubdomainID> & block_ids,
250 : std::unordered_set<const Elem *> & wall_bounded)
251 : {
252 :
253 : wall_bounded.clear();
254 18621 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_names);
255 :
256 : // We define these lambdas so that we can fetch the bounded elements from other
257 : // processors.
258 14120 : auto gather_functor = [&subproblem, &wall_bounded](const processor_id_type libmesh_dbg_var(pid),
259 : const std::vector<dof_id_type> & elem_ids,
260 : std::vector<unsigned char> & data_to_fill)
261 : {
262 : mooseAssert(pid != subproblem.processor_id(), "We shouldn't be gathering from ourselves.");
263 14120 : data_to_fill.resize(elem_ids.size());
264 :
265 14120 : const auto & mesh = subproblem.mesh().getMesh();
266 :
267 172472 : for (const auto i : index_range(elem_ids))
268 : {
269 158352 : const auto elem = mesh.elem_ptr(elem_ids[i]);
270 316704 : data_to_fill[i] = wall_bounded.count(elem) != 0;
271 : }
272 32741 : };
273 :
274 14120 : auto action_functor = [&subproblem, &wall_bounded](const processor_id_type libmesh_dbg_var(pid),
275 : const std::vector<dof_id_type> & elem_ids,
276 : const std::vector<unsigned char> & filled_data)
277 : {
278 : mooseAssert(pid != subproblem.processor_id(),
279 : "The request filler shouldn't have been ourselves");
280 : mooseAssert(elem_ids.size() == filled_data.size(), "I think these should be the same size");
281 :
282 14120 : const auto & mesh = subproblem.mesh().getMesh();
283 :
284 172472 : for (const auto i : index_range(elem_ids))
285 : {
286 158352 : const auto elem = mesh.elem_ptr(elem_ids[i]);
287 158352 : if (filled_data[i])
288 14432 : wall_bounded.insert(elem);
289 : }
290 32741 : };
291 :
292 : // We need these elements from other processors
293 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> elem_ids_requested;
294 :
295 4727410 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
296 2345084 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
297 11695180 : for (const auto i_side : elem->side_index_range())
298 : {
299 : // This is needed because in some cases the internal boundary is registered
300 : // to the neighbor element
301 : std::set<BoundaryID> combined_side_bds;
302 9356144 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
303 9356144 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
304 9356144 : if (const auto neighbor = elem->neighbor_ptr(i_side))
305 : {
306 8752224 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
307 8752224 : const auto & neighbor_bnds = subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
308 8752224 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
309 :
310 : // If the neighbor lives on the first layer of the ghost region then we would
311 : // like to grab its value as well (if it exists)
312 8752224 : if (neighbor->processor_id() != subproblem.processor_id() &&
313 : block_ids.find(neighbor->subdomain_id()) != block_ids.end())
314 158352 : elem_ids_requested[neighbor->processor_id()].push_back(neighbor->id());
315 8752224 : }
316 :
317 25680252 : for (const auto wall_id : wall_boundary_ids)
318 : if (combined_side_bds.count(wall_id))
319 : {
320 : wall_bounded.insert(elem);
321 : break;
322 : }
323 9374765 : }
324 :
325 : unsigned char * bool_ex = nullptr;
326 18621 : TIMPI::pull_parallel_vector_data(
327 : subproblem.comm(), elem_ids_requested, gather_functor, action_functor, bool_ex);
328 18621 : }
329 :
330 : /// Bounded element face distances for wall treatment
331 : void
332 5049 : getWallDistance(const std::vector<BoundaryName> & wall_boundary_name,
333 : const FEProblemBase & fe_problem,
334 : const SubProblem & subproblem,
335 : const std::set<SubdomainID> & block_ids,
336 : std::map<const Elem *, std::vector<Real>> & dist_map)
337 : {
338 : dist_map.clear();
339 5049 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_name);
340 :
341 1321018 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
342 655460 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
343 3268660 : for (const auto i_side : elem->side_index_range())
344 : {
345 : // This is needed because in some cases the internal boundary is registered
346 : // to the neighbor element
347 : std::set<BoundaryID> combined_side_bds;
348 2614928 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
349 2614928 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
350 2614928 : if (const auto neighbor = elem->neighbor_ptr(i_side))
351 : {
352 2449272 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
353 : const std::vector<BoundaryID> & neighbor_bnds =
354 2449272 : subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
355 2449272 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
356 2449272 : }
357 :
358 10016208 : for (const auto wall_id : wall_boundary_ids)
359 : if (combined_side_bds.count(wall_id))
360 : {
361 : // The list below stores the face infos with respect to their owning elements,
362 : // depending on the block restriction we might encounter situations where the
363 : // element outside of the block owns the face info.
364 138016 : const auto & neighbor = elem->neighbor_ptr(i_side);
365 138016 : const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
366 138016 : const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
367 288 : const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
368 :
369 138016 : const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
370 138016 : const auto & elem_centroid = elem_has_fi ? fi->elemCentroid() : fi->neighborCentroid();
371 138016 : const Real dist = std::abs((elem_centroid - fi->faceCentroid()) * fi->normal());
372 138016 : dist_map[elem].push_back(dist);
373 : }
374 2619977 : }
375 5049 : }
376 :
377 : /// Face arguments to wall-bounded faces for wall treatment
378 : void
379 5049 : getElementFaceArgs(const std::vector<BoundaryName> & wall_boundary_name,
380 : const FEProblemBase & fe_problem,
381 : const SubProblem & subproblem,
382 : const std::set<SubdomainID> & block_ids,
383 : std::map<const Elem *, std::vector<const FaceInfo *>> & face_info_map)
384 : {
385 : face_info_map.clear();
386 5049 : const auto wall_boundary_ids = subproblem.mesh().getBoundaryIDs(wall_boundary_name);
387 :
388 1321018 : for (const auto & elem : fe_problem.mesh().getMesh().active_local_element_ptr_range())
389 655460 : if (block_ids.find(elem->subdomain_id()) != block_ids.end())
390 3268660 : for (const auto i_side : elem->side_index_range())
391 : {
392 : // This is needed because in some cases the internal boundary is registered
393 : // to the neighbor element
394 : std::set<BoundaryID> combined_side_bds;
395 2614928 : const auto & side_bnds = subproblem.mesh().getBoundaryIDs(elem, i_side);
396 2614928 : combined_side_bds.insert(side_bnds.begin(), side_bnds.end());
397 2614928 : if (elem->neighbor_ptr(i_side) && !elem->neighbor_ptr(i_side)->is_remote())
398 : {
399 2449272 : const auto neighbor = elem->neighbor_ptr(i_side);
400 2449272 : const auto neighbor_side = neighbor->which_neighbor_am_i(elem);
401 : const std::vector<BoundaryID> & neighbor_bnds =
402 2449272 : subproblem.mesh().getBoundaryIDs(neighbor, neighbor_side);
403 2449272 : combined_side_bds.insert(neighbor_bnds.begin(), neighbor_bnds.end());
404 2449272 : }
405 :
406 10016208 : for (const auto wall_id : wall_boundary_ids)
407 : if (combined_side_bds.count(wall_id))
408 : {
409 : // The list below stores the face infos with respect to their owning elements,
410 : // depending on the block restriction we might encounter situations where the
411 : // element outside of the block owns the face info.
412 138016 : const auto & neighbor = elem->neighbor_ptr(i_side);
413 138016 : const auto elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
414 138016 : const auto & elem_for_fi = elem_has_fi ? elem : neighbor;
415 288 : const auto side = elem_has_fi ? i_side : neighbor->which_neighbor_am_i(elem);
416 :
417 138016 : const FaceInfo * const fi = subproblem.mesh().faceInfo(elem_for_fi, side);
418 138016 : face_info_map[elem].push_back(fi);
419 : }
420 2619977 : }
421 5049 : }
422 : }
|
