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