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