/home/runner/work/amr-wind/amr-wind/amr-wind/equation_systems/vof/volume_fractions.H Source File

AMR-Wind API: /home/runner/work/amr-wind/amr-wind/amr-wind/equation_systems/vof/volume_fractions.H Source File
AMR-Wind API v0.1.0
CFD solver for wind plant simulations
Loading...
Searching...
No Matches
volume_fractions.H
Go to the documentation of this file.
1#ifndef VOLUME_FRACTIONS_H_
2#define VOLUME_FRACTIONS_H_
3
4#include <AMReX_FArrayBox.H>
5#include <cmath>
6#include "amr-wind/utilities/constants.H"
7#include "AMReX_REAL.H"
8
9using namespace amrex::literals;
10
11namespace amr_wind::multiphase {
12
22
23AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_finite_difference_normal(
24 const int i,
25 const int j,
26 const int k,
27 amrex::Array4<amrex::Real const> const& volfrac,
28 amrex::Real& mx,
29 amrex::Real& my,
30 amrex::Real& mz) noexcept
31{
32 amrex::Real mm1 =
33 volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
34 volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
35 2.0_rt * (volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
36 volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1)) +
37 4.0_rt * volfrac(i - 1, j, k);
38 amrex::Real mm2 =
39 volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
40 volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
41 2.0_rt * (volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
42 volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1)) +
43 4.0_rt * volfrac(i + 1, j, k);
44 mx = mm1 - mm2;
45
46 mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
47 volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
48 2.0_rt * (volfrac(i - 1, j - 1, k) + volfrac(i + 1, j - 1, k) +
49 volfrac(i, j - 1, k - 1) + volfrac(i, j - 1, k + 1)) +
50 4.0_rt * volfrac(i, j - 1, k);
51 mm2 = volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
52 volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
53 2.0_rt * (volfrac(i - 1, j + 1, k) + volfrac(i + 1, j + 1, k) +
54 volfrac(i, j + 1, k - 1) + volfrac(i, j + 1, k + 1)) +
55 4.0_rt * volfrac(i, j + 1, k);
56 my = mm1 - mm2;
57
58 mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j + 1, k - 1) +
59 volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j + 1, k - 1) +
60 2.0_rt * (volfrac(i - 1, j, k - 1) + volfrac(i + 1, j, k - 1) +
61 volfrac(i, j - 1, k - 1) + volfrac(i, j + 1, k - 1)) +
62 4.0_rt * volfrac(i, j, k - 1);
63 mm2 = volfrac(i - 1, j - 1, k + 1) + volfrac(i - 1, j + 1, k + 1) +
64 volfrac(i + 1, j - 1, k + 1) + volfrac(i + 1, j + 1, k + 1) +
65 2.0_rt * (volfrac(i - 1, j, k + 1) + volfrac(i + 1, j, k + 1) +
66 volfrac(i, j - 1, k + 1) + volfrac(i, j + 1, k + 1)) +
67 4.0_rt * volfrac(i, j, k + 1);
68 mz = mm1 - mm2;
69}
70
71AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
72youngs_finite_difference_normal_neumann(
73 const int i,
74 const int j,
75 const int k,
76 const int ibdy,
77 const int jbdy,
78 const int kbdy,
79 amrex::Array4<amrex::Real const> const& volfrac,
80 amrex::Real& mx,
81 amrex::Real& my,
82 amrex::Real& mz) noexcept
83{
84 // Do neumann condition via indices
85 const int im1 = ibdy == -1 ? i : i - 1;
86 const int jm1 = jbdy == -1 ? j : j - 1;
87 const int km1 = kbdy == -1 ? k : k - 1;
88 const int ip1 = ibdy == +1 ? i : i + 1;
89 const int jp1 = jbdy == +1 ? j : j + 1;
90 const int kp1 = kbdy == +1 ? k : k + 1;
91
92 amrex::Real mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
93 volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
94 2.0_rt * (volfrac(im1, jm1, k) + volfrac(im1, jp1, k) +
95 volfrac(im1, j, km1) + volfrac(im1, j, kp1)) +
96 4.0_rt * volfrac(im1, j, k);
97 amrex::Real mm2 = volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
98 volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
99 2.0_rt * (volfrac(ip1, jm1, k) + volfrac(ip1, jp1, k) +
100 volfrac(ip1, j, km1) + volfrac(ip1, j, kp1)) +
101 4.0_rt * volfrac(ip1, j, k);
102 mx = mm1 - mm2;
103
104 mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
105 volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
106 2.0_rt * (volfrac(im1, jm1, k) + volfrac(ip1, jm1, k) +
107 volfrac(i, jm1, km1) + volfrac(i, jm1, kp1)) +
108 4.0_rt * volfrac(i, jm1, k);
109 mm2 = volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
110 volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
111 2.0_rt * (volfrac(im1, jp1, k) + volfrac(ip1, jp1, k) +
112 volfrac(i, jp1, km1) + volfrac(i, jp1, kp1)) +
113 4.0_rt * volfrac(i, jp1, k);
114 my = mm1 - mm2;
115
116 mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jp1, km1) +
117 volfrac(ip1, jm1, km1) + volfrac(ip1, jp1, km1) +
118 2.0_rt * (volfrac(im1, j, km1) + volfrac(ip1, j, km1) +
119 volfrac(i, jm1, km1) + volfrac(i, jp1, km1)) +
120 4.0_rt * volfrac(i, j, km1);
121 mm2 = volfrac(im1, jm1, kp1) + volfrac(im1, jp1, kp1) +
122 volfrac(ip1, jm1, kp1) + volfrac(ip1, jp1, kp1) +
123 2.0_rt * (volfrac(im1, j, kp1) + volfrac(ip1, j, kp1) +
124 volfrac(i, jm1, kp1) + volfrac(i, jp1, kp1)) +
125 4.0_rt * volfrac(i, j, kp1);
126 mz = mm1 - mm2;
127}
128
129AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
130 const int i,
131 const int j,
132 const int k,
133 amrex::Array4<amrex::Real const> const& volfrac,
134 amrex::Real& mx,
135 amrex::Real& my,
136 amrex::Real& mz) noexcept
137{
138 amrex::Array2D<amrex::Real, 0, AMREX_SPACEDIM + 1, 0, AMREX_SPACEDIM> m;
139 // write the plane as: sgn(mx) X = my Y + mz Z + alpha
140 // m00 X = m01 Y + m02 Z + alpha
141 amrex::Real m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
142 volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
143 volfrac(i - 1, j, k);
144 amrex::Real m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
145 volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
146 volfrac(i + 1, j, k);
147
148 m(0, 0) = (m1 > m2) ? 1.0_rt : -1.0_rt;
149
150 m1 = volfrac(i - 1, j - 1, k) + volfrac(i + 1, j - 1, k) +
151 volfrac(i, j - 1, k);
152 m2 = volfrac(i - 1, j + 1, k) + volfrac(i + 1, j + 1, k) +
153 volfrac(i, j + 1, k);
154 m(0, 1) = 0.5_rt * (m1 - m2);
155
156 m1 = volfrac(i - 1, j, k - 1) + volfrac(i + 1, j, k - 1) +
157 volfrac(i, j, k - 1);
158 m2 = volfrac(i - 1, j, k + 1) + volfrac(i + 1, j, k + 1) +
159 volfrac(i, j, k + 1);
160 m(0, 2) = 0.5_rt * (m1 - m2);
161
162 // write the plane as: sgn(my) Y = mx X + mz Z + alpha,
163 // m11 Y = m10 X + m12 Z + alpha.
164
165 m1 = volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
166 volfrac(i - 1, j, k);
167 m2 = volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
168 volfrac(i + 1, j, k);
169 m(1, 0) = 0.5_rt * (m1 - m2);
170
171 m1 = volfrac(i, j - 1, k - 1) + volfrac(i, j - 1, k + 1) +
172 volfrac(i + 1, j - 1, k) + volfrac(i - 1, j - 1, k) +
173 volfrac(i, j - 1, k);
174 m2 = volfrac(i, j + 1, k - 1) + volfrac(i, j + 1, k + 1) +
175 volfrac(i + 1, j + 1, k) + volfrac(i - 1, j + 1, k) +
176 volfrac(i, j + 1, k);
177
178 m(1, 1) = (m1 > m2) ? 1.0_rt : -1.0_rt;
179
180 m1 = volfrac(i, j - 1, k - 1) + volfrac(i, j, k - 1) +
181 volfrac(i, j + 1, k - 1);
182 m2 = volfrac(i, j - 1, k + 1) + volfrac(i, j, k + 1) +
183 volfrac(i, j + 1, k + 1);
184 m(1, 2) = 0.5_rt * (m1 - m2);
185
186 // write the plane as: sgn(mz) Z = mx X + my Y + alpha
187 // m22 Z = m20 X + m21 Y + alpha
188
189 m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
190 volfrac(i - 1, j, k);
191 m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
192 volfrac(i + 1, j, k);
193 m(2, 0) = 0.5_rt * (m1 - m2);
194
195 m1 = volfrac(i, j - 1, k - 1) + volfrac(i, j - 1, k + 1) +
196 volfrac(i, j - 1, k);
197 m2 = volfrac(i, j + 1, k - 1) + volfrac(i, j + 1, k + 1) +
198 volfrac(i, j + 1, k);
199 m(2, 1) = 0.5_rt * (m1 - m2);
200
201 m1 = volfrac(i - 1, j, k - 1) + volfrac(i + 1, j, k - 1) +
202 volfrac(i, j - 1, k - 1) + volfrac(i, j + 1, k - 1) +
203 volfrac(i, j, k - 1);
204 m2 = volfrac(i - 1, j, k + 1) + volfrac(i + 1, j, k + 1) +
205 volfrac(i, j - 1, k + 1) + volfrac(i, j + 1, k + 1) +
206 volfrac(i, j, k + 1);
207
208 m(2, 2) = (m1 > m2) ? 1.0_rt : -1.0_rt;
209
210 // normalize each set (mx,my,mz): |mx|+|my|+|mz| = 1
211 amrex::Real t0, t1, t2;
212
213 t0 = std::abs(m(0, 0)) + std::abs(m(0, 1)) + std::abs(m(0, 2));
214 m(0, 0) = m(0, 0) / t0;
215 m(0, 1) = m(0, 1) / t0;
216 m(0, 2) = m(0, 2) / t0;
217
218 t0 = std::abs(m(1, 0)) + std::abs(m(1, 1)) + std::abs(m(1, 2));
219 m(1, 0) = m(1, 0) / t0;
220 m(1, 1) = m(1, 1) / t0;
221 m(1, 2) = m(1, 2) / t0;
222
223 t0 = std::abs(m(2, 0)) + std::abs(m(2, 1)) + std::abs(m(2, 2));
224 m(2, 0) = m(2, 0) / t0;
225 m(2, 1) = m(2, 1) / t0;
226 m(2, 2) = m(2, 2) / t0;
227
228 // choose among the three central schemes */
229 t0 = std::abs(m(0, 0));
230 t1 = std::abs(m(1, 1));
231 t2 = std::abs(m(2, 2));
232
233 int cn = 0;
234 if (t1 > t0) {
235 t0 = t1;
236 cn = 1;
237 }
238
239 if (t2 > t0) {
240 cn = 2;
241 }
242
243 // Youngs-CIAM scheme */
244 youngs_finite_difference_normal(
245 i, j, k, volfrac, m(3, 0), m(3, 1), m(3, 2));
246 // normalize the set (mx,my,mz): |mx|+|my|+|mz| = 1
247 constexpr amrex::Real tiny = std::numeric_limits<amrex::Real>::epsilon();
248 t0 = std::abs(m(3, 0)) + std::abs(m(3, 1)) + std::abs(m(3, 2)) + tiny;
249 m(3, 0) = m(3, 0) / t0;
250 m(3, 1) = m(3, 1) / t0;
251 m(3, 2) = m(3, 2) / t0;
252
253 // choose between the previous choice and Youngs-CIAM
254 t0 = std::abs(m(3, 0));
255 t1 = std::abs(m(3, 1));
256 t2 = std::abs(m(3, 2));
257 if (t1 > t0) {
258 t0 = t1;
259 }
260 if (t2 > t0) {
261 t0 = t2;
262 }
263
264 if (std::abs(m(cn, cn)) > t0 && t0 > 0.0_rt) {
265 // second t0 condition is to ensure nonzero normal magnitude
266 cn = 3;
267 }
268
269 // components of the normal vector */
270 mx = m(cn, 0);
271 my = m(cn, 1);
272 mz = m(cn, 2);
273}
274
275/* Computes alpha: m1*x1 + m2* x2 + m3*x3 = alpha
276 * given that m1+m2+m3=1 (m1,m2,m3>0) and the volumetric fraction volF
277 */
278AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(
279 const amrex::Real b1,
280 const amrex::Real b2,
281 const amrex::Real b3,
282 const amrex::Real volF) noexcept
283{
284 amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
285 // (1) "order coefficients: m1<m2<m3"
286 // (2) "get ranges: V1<V2<v3"
287 // (3) "limit ch (0.d0 < ch < 0.5d0)"
288 // (4) "calculate alpha"
289
290 amrex::Real m1 = amrex::min(b1, b2);
291 amrex::Real m3 = amrex::max(b1, b2);
292 amrex::Real m2 = b3;
293
294 amrex::Real tmp;
295 if (m2 < m1) {
296 tmp = m1;
297 m1 = m2;
298 m2 = tmp;
299 } else if (m2 > m3) {
300 tmp = m3;
301 m3 = m2;
302 m2 = tmp;
303 }
304
305 amrex::Real m12 = m1 + m2;
306 amrex::Real pr = amrex::max<amrex::Real>(6.0_rt * m1 * m2 * m3, const_tiny);
307 amrex::Real V1 = m1 * m1 * m1 / pr;
308 amrex::Real V2 = V1 + 0.5_rt * (m2 - m1) / m3;
309
310 amrex::Real mm, V3;
311 if (m3 < m12) {
312 mm = m3;
313 V3 = (m3 * m3 * (3.0_rt * m12 - m3) + m1 * m1 * (m1 - 3.0_rt * m3) +
314 m2 * m2 * (m2 - 3.0_rt * m3)) /
315 pr;
316 } else {
317 mm = m12;
318 V3 = 0.5_rt * mm / m3;
319 }
320
321 amrex::Real ch = amrex::min<amrex::Real>(volF, 1.0_rt - volF);
322
323 amrex::Real alpha, p, q, p12, teta, cs;
324 if (ch < V1) {
325 alpha = std::cbrt(pr * ch); // case (1)
326 } else if (ch < V2) {
327 alpha =
328 0.50_rt * (m1 + std::sqrt(m1 * m1 + 8.0_rt * m2 * m3 * (ch - V1)));
329 } else if (ch < V3) {
330 p = 2.0_rt * m1 * m2;
331 q = 1.5_rt * m1 * m2 * (m12 - 2.0_rt * m3 * ch);
332 p12 = std::sqrt(p);
333 teta =
334 std::acos(
335 q / (p * p12 + std::numeric_limits<amrex::Real>::epsilon())) /
336 3.0_rt;
337 cs = std::cos(teta);
338 alpha = p12 * (std::sqrt(3.0_rt * (1.0_rt - cs * cs)) - cs) + m12;
339 } else if (m12 < m3) {
340 alpha = m3 * ch + 0.5_rt * mm;
341 } else {
342 p = m1 * (m2 + m3) + m2 * m3 - 0.25_rt;
343 q = 1.5_rt * m1 * m2 * m3 * (0.5_rt - ch);
344 p12 = std::sqrt(p);
345 teta =
346 std::acos(
347 q / (p * p12 + std::numeric_limits<amrex::Real>::epsilon())) /
348 3.0_rt;
349 cs = std::cos(teta);
350 alpha = p12 * (std::sqrt(3.0_rt * (1.0_rt - cs * cs)) - cs) + 0.5_rt;
351 }
352
353 if (volF > 0.5_rt) {
354 alpha = 1.0_rt - alpha;
355 }
356
357 return alpha;
358}
359
368AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real cut_volume(
369 const amrex::Real m1,
370 const amrex::Real m2,
371 const amrex::Real m3,
372 const amrex::Real alpha,
373 const amrex::Real r0,
374 const amrex::Real dr0) noexcept
375{
376 const amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
377 amrex::Real al;
378
379 // move origin to x0
380 al = alpha - m1 * r0;
381
382 // reflect the figure when negative coefficients
383 al = al + amrex::max<amrex::Real>(0.0_rt, -m1 * dr0) +
384 amrex::max<amrex::Real>(0.0_rt, -m2) +
385 amrex::max<amrex::Real>(0.0_rt, -m3);
386
387 // normalized equation: m1*y1 + m2*y2 + m3*y3 = alh, with 0 <= m1 <= m2 <=
388 // m3 the problem is then solved again in the unit cube
389 amrex::Real tmp = std::abs(m1) * dr0 + std::abs(m2) + std::abs(m3);
390 amrex::Real n1 = std::abs(m1) / tmp; // need positive coefficients
391 amrex::Real n2 = std::abs(m2) / tmp;
392 amrex::Real n3 = std::abs(m3) / tmp;
393 al = amrex::max<amrex::Real>(
394 0.0_rt, amrex::min<amrex::Real>(
395 1.0_rt, al / tmp)); // get new al within safe limits
396 amrex::Real al0 =
397 amrex::min<amrex::Real>(al, 1.0_rt - al); // limit to: 0 < alh < 1/2
398
399 // Order coefficients
400 amrex::Real b1 = amrex::min(n1 * dr0, n2); // order coefficients
401 amrex::Real b3 = amrex::max(n1 * dr0, n2);
402 amrex::Real b2 = n3;
403
404 if (b2 < b1) {
405 tmp = b1;
406 b1 = b2;
407 b2 = tmp;
408 } else if (b2 > b3) {
409 tmp = b3;
410 b3 = b2;
411 b2 = tmp;
412 }
413
414 amrex::Real b12 = b1 + b2;
415
416 // Compute volume fraction using PLIC from Scardovelli & Zaleski (JCP 2000),
417 // adapted from code in paper by Akio Kawano (Computer & Fluids 2016)
418 amrex::Real vm1 = b1;
419 amrex::Real vm3 = b3;
420 amrex::Real vm2 = b2;
421 amrex::Real vm12 = b12;
422 amrex::Real a = al0;
423 amrex::Real v = 0.0_rt;
424
425 if (a > 0.0_rt) {
426 if (a < vm1) {
427 v = a * a * a / (6.0_rt * vm1 * vm2 * vm3);
428 } else if (a < vm2) {
429 v = a * (a - vm1) / (2.0_rt * vm2 * vm3) +
430 vm1 * vm1 / (6.0_rt * vm2 * vm3 + const_tiny);
431 } else if (a < amrex::min(vm12, vm3)) {
432 v = (a * a * (3.0_rt * vm12 - a) + vm1 * vm1 * (vm1 - 3.0_rt * a) +
433 vm2 * vm2 * (vm2 - 3.0_rt * a)) /
434 (6.0_rt * vm1 * vm2 * vm3);
435 } else if (vm3 < vm12) {
436 v = (a * a * (3.0_rt - 2.0_rt * a) +
437 vm1 * vm1 * (vm1 - 3.0_rt * a) +
438 vm2 * vm2 * (vm2 - 3.0_rt * a) +
439 vm3 * vm3 * (vm3 - 3.0_rt * a)) /
440 (6.0_rt * vm1 * vm2 * vm3);
441 } else {
442 v = (a - 0.5_rt * vm12) / vm3;
443 }
444 }
445
446 tmp = v;
447 amrex::Real FL3D;
448 if (al <= 0.5_rt) {
449 FL3D = tmp;
450 } else {
451 FL3D = (1.0_rt - tmp);
452 }
453
454 return FL3D;
455}
456
457AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(
458 const int i,
459 const int j,
460 const int k,
461 amrex::Array4<amrex::Real const> const& volfrac,
462 amrex::Real& mx,
463 amrex::Real& my,
464 amrex::Real& mz,
465 amrex::Real& alpha)
466{
467
468 mixed_youngs_central_normal(i, j, k, volfrac, mx, my, mz);
469
470 amrex::Real invx = 1.0_rt;
471 amrex::Real invy = 1.0_rt;
472 amrex::Real invz = 1.0_rt;
473
474 if (mx < 0.0_rt) {
475 mx = -mx;
476 invx = -1.0_rt;
477 }
478 if (my < 0.0_rt) {
479 my = -my;
480 invy = -1.0_rt;
481 }
482 if (mz < 0.0_rt) {
483 mz = -mz;
484 invz = -1.0_rt;
485 }
486
487 amrex::Real mm2 = mx + my + mz;
488 mx = mx / mm2;
489 my = my / mm2;
490 mz = mz / mm2;
491
492 alpha = volume_intercept(mx, my, mz, volfrac(i, j, k));
493
494 // Back to the original plane
495 mx = invx * mx;
496 my = invy * my;
497 mz = invz * mz;
498 alpha = alpha + amrex::min<amrex::Real>(0.0_rt, mx) +
499 amrex::min<amrex::Real>(0.0_rt, my) +
500 amrex::min<amrex::Real>(0.0_rt, mz);
501}
502
503AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
504 const int i,
505 const int j,
506 const int k,
507 amrex::Array4<amrex::Real const> const& volfrac,
508 const int n_band = 1,
509 const amrex::Real tiny = constants::TIGHT_TOL) noexcept
510{
511 // n_band must be <= number of vof ghost cells (3)
512
513 // Check if near interface
514 amrex::Real VOF_max = 0.0_rt;
515 amrex::Real VOF_min = 1.0_rt;
516 bool VOF_mid = false;
517 for (int ii = -n_band; ii <= n_band; ++ii) {
518 for (int jj = -n_band; jj <= n_band; ++jj) {
519 for (int kk = -n_band; kk <= n_band; ++kk) {
520 amrex::Real VOF = volfrac(i + ii, j + jj, k + kk);
521 VOF_max = amrex::max(VOF_max, VOF);
522 VOF_min = amrex::min(VOF_min, VOF);
523 if (VOF < 1.0_rt - tiny && VOF > tiny) {
524 VOF_mid = true;
525 }
526 }
527 }
528 }
529 return (std::abs(VOF_max - VOF_min) > tiny || VOF_mid);
530}
531
532AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real levelset_to_vof(
533 const int i,
534 const int j,
535 const int k,
536 const amrex::Real eps,
537 amrex::Array4<amrex::Real const> const& phi) noexcept
538{
539 amrex::Real mx, my, mz;
540 youngs_finite_difference_normal(i, j, k, phi, mx, my, mz);
541 mx = std::abs(mx / 32.0_rt);
542 my = std::abs(my / 32.0_rt);
543 mz = std::abs(mz / 32.0_rt);
544 amrex::Real normL1 = (mx + my + mz);
545 mx = mx / normL1;
546 my = my / normL1;
547 mz = mz / normL1;
548 // Make sure that alpha is negative far away from the
549 // interface
550 amrex::Real alpha;
551 if (phi(i, j, k) < -eps) {
552 alpha = -1.0_rt;
553 } else {
554 alpha = phi(i, j, k) / normL1;
555 alpha = alpha + 0.5_rt;
556 }
557 amrex::Real result;
558 if (alpha >= 1.0_rt) {
559 result = 1.0_rt;
560 } else if (alpha <= 0.0_rt) {
561 result = 0.0_rt;
562 } else {
563 result = cut_volume(mx, my, mz, alpha, 0.0_rt, 1.0_rt);
564 }
565
566 return result;
567}
568
569} // namespace amr_wind::multiphase
570
571#endif // VOLUME_FRACTIONS.H
amr_wind::constants::TIGHT_TOL
static constexpr amrex::Real TIGHT_TOL
A tight tolerance.
Definition constants.H:17
amr_wind::multiphase
Definition height_functions.H:8
amr_wind::multiphase::volume_intercept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(const amrex::Real b1, const amrex::Real b2, const amrex::Real b3, const amrex::Real volF) noexcept
Definition volume_fractions.H:278
amr_wind::multiphase::fit_plane
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(const int i, const int j, const int k, amrex::Array4< amrex::Real const > const &volfrac, amrex::Real &mx, amrex::Real &my, amrex::Real &mz, amrex::Real &alpha)
Definition volume_fractions.H:457
amr_wind::multiphase::interface_band
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(const int i, const int j, const int k, amrex::Array4< amrex::Real const > const &volfrac, const int n_band=1, const amrex::Real tiny=constants::TIGHT_TOL) noexcept
Definition volume_fractions.H:503
amr_wind::multiphase::levelset_to_vof
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real levelset_to_vof(const int i, const int j, const int k, const amrex::Real eps, amrex::Array4< amrex::Real const > const &phi) noexcept
Definition volume_fractions.H:532
amr_wind::multiphase::cut_volume
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real cut_volume(const amrex::Real m1, const amrex::Real m2, const amrex::Real m3, const amrex::Real alpha, const amrex::Real r0, const amrex::Real dr0) noexcept
Definition volume_fractions.H:368
amr_wind::multiphase::youngs_finite_difference_normal
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_finite_difference_normal(const int i, const int j, const int k, amrex::Array4< amrex::Real const > const &volfrac, amrex::Real &mx, amrex::Real &my, amrex::Real &mz) noexcept
Definition volume_fractions.H:23
amr_wind::multiphase::mixed_youngs_central_normal
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(const int i, const int j, const int k, amrex::Array4< amrex::Real const > const &volfrac, amrex::Real &mx, amrex::Real &my, amrex::Real &mz) noexcept
Definition volume_fractions.H:129
amr_wind::multiphase::youngs_finite_difference_normal_neumann
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_finite_difference_normal_neumann(const int i, const int j, const int k, const int ibdy, const int jbdy, const int kbdy, amrex::Array4< amrex::Real const > const &volfrac, amrex::Real &mx, amrex::Real &my, amrex::Real &mz) noexcept
Definition volume_fractions.H:72
amr_wind::InterfaceCapturingMethod::VOF
@ VOF
Volume of fluid.
Definition MultiPhase.H:25
Generated by doxygen 1.14.0