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