/home/runner/work/kynema-sgf/kynema-sgf/src/equation_systems/vof/volume_fractions.H Source File

Kynema-SGF API: /home/runner/work/kynema-sgf/kynema-sgf/src/equation_systems/vof/volume_fractions.H Source File
Kynema-SGF 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 "src/utilities/constants.H"
8#include "AMReX_REAL.H"
9
10using namespace amrex::literals;
11
12namespace kynema_sgf::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 constexpr amrex::Real 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, tiny);
304 amrex::Real V1 = m1 * m1 * m1 / pr;
305 amrex::Real V2 = V1 + (0.5_rt * (m2 - m1) / (m3 + tiny));
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 + tiny);
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(amrex::max(0.0_rt, p));
330 teta =
331 std::acos(
332 amrex::max(-1.0_rt, amrex::min(1.0_rt, q / (p * p12 + tiny)))) /
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(amrex::max(0.0_rt, p));
342 teta =
343 std::acos(
344 amrex::max(-1.0_rt, amrex::min(1.0_rt, q / (p * p12 + tiny)))) /
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 constexpr amrex::Real 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 = amrex::max<amrex::Real>(
387 tiny, (std::abs(m1) * dr0) + std::abs(m2) + std::abs(m3));
388 amrex::Real n1 = std::abs(m1) / tmp; // need positive coefficients
389 amrex::Real n2 = std::abs(m2) / tmp;
390 amrex::Real n3 = std::abs(m3) / tmp;
391 al = amrex::max<amrex::Real>(
392 0.0_rt, amrex::min<amrex::Real>(
393 1.0_rt, al / tmp)); // get new al within safe limits
394 amrex::Real al0 =
395 amrex::min<amrex::Real>(al, 1.0_rt - al); // limit to: 0 < alh < 1/2
396
397 // Order coefficients
398 amrex::Real b1 =
399 amrex::min<amrex::Real>(n1 * dr0, n2); // order coefficients
400 amrex::Real b3 = amrex::max<amrex::Real>(n1 * dr0, n2);
401 amrex::Real b2 = n3;
402
403 if (b2 < b1) {
404 tmp = b1;
405 b1 = b2;
406 b2 = tmp;
407 } else if (b2 > b3) {
408 tmp = b3;
409 b3 = b2;
410 b2 = tmp;
411 }
412
413 amrex::Real b12 = b1 + b2;
414
415 // Compute volume fraction using PLIC from Scardovelli & Zaleski (JCP 2000),
416 // adapted from code in paper by Akio Kawano (Computer & Fluids 2016)
417 amrex::Real vm1 = b1;
418 amrex::Real vm3 = b3;
419 amrex::Real vm2 = b2;
420 amrex::Real vm12 = b12;
421 amrex::Real a = al0;
422 amrex::Real v = 0.0_rt;
423
424 if (a > 0.0_rt) {
425 if (a < vm1) {
426 v = a * a * a / (6.0_rt * vm1 * vm2 * vm3);
427 } else if (a < vm2) {
428 v = (a * (a - vm1) / (2.0_rt * vm2 * vm3)) +
429 (vm1 * vm1 / (6.0_rt * vm2 * vm3 + tiny));
430 } else if (a < amrex::min(vm12, vm3)) {
431 v = (a * a * (3.0_rt * vm12 - a) + vm1 * vm1 * (vm1 - 3.0_rt * a) +
432 vm2 * vm2 * (vm2 - 3.0_rt * a)) /
433 (6.0_rt * vm1 * vm2 * vm3);
434 } else if (vm3 < vm12) {
435 v = (a * a * (3.0_rt - 2.0_rt * a) +
436 vm1 * vm1 * (vm1 - 3.0_rt * a) +
437 vm2 * vm2 * (vm2 - 3.0_rt * a) +
438 vm3 * vm3 * (vm3 - 3.0_rt * a)) /
439 (6.0_rt * vm1 * vm2 * vm3);
440 } else {
441 v = (a - 0.5_rt * vm12) / vm3;
442 }
443 }
444
445 tmp = v;
446 amrex::Real FL3D;
447 if (al <= 0.5_rt) {
448 FL3D = tmp;
449 } else {
450 FL3D = (1.0_rt - tmp);
451 }
452
453 return FL3D;
454}
455
456AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(
457 const int i,
458 const int j,
459 const int k,
460 amrex::Array4<amrex::Real const> const& volfrac,
461 amrex::Real& mx,
462 amrex::Real& my,
463 amrex::Real& mz,
464 amrex::Real& alpha)
465{
466
467 mixed_youngs_central_normal(i, j, k, volfrac, mx, my, mz);
468
469 amrex::Real invx = 1.0_rt;
470 amrex::Real invy = 1.0_rt;
471 amrex::Real invz = 1.0_rt;
472
473 if (mx < 0.0_rt) {
474 mx = -mx;
475 invx = -1.0_rt;
476 }
477 if (my < 0.0_rt) {
478 my = -my;
479 invy = -1.0_rt;
480 }
481 if (mz < 0.0_rt) {
482 mz = -mz;
483 invz = -1.0_rt;
484 }
485
486 amrex::Real mm2 = amrex::max<amrex::Real>(
487 std::numeric_limits<amrex::Real>::epsilon(), 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)
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<amrex::Real>(VOF_max, VOF);
522 VOF_min = amrex::min<amrex::Real>(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)
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 const amrex::Real norm_raw = (mx + my + mz);
545 if (norm_raw <= std::numeric_limits<amrex::Real>::epsilon()) {
546 const amrex::Real phi_ijk = phi(i, j, k);
547 if (phi_ijk < -eps) {
548 return 0.0_rt;
549 }
550 if (phi_ijk > eps) {
551 return 1.0_rt;
552 }
553 return 0.5_rt;
554 }
555 amrex::Real normL1 = norm_raw;
556 mx = mx / normL1;
557 my = my / normL1;
558 mz = mz / normL1;
559 // Make sure that alpha is negative far away from the
560 // interface
561 amrex::Real alpha;
562 if (phi(i, j, k) < -eps) {
563 alpha = -1.0_rt;
564 } else {
565 alpha = phi(i, j, k) / normL1;
566 alpha = alpha + 0.5_rt;
567 }
568 amrex::Real result;
569 if (alpha >= 1.0_rt) {
570 result = 1.0_rt;
571 } else if (alpha <= 0.0_rt) {
572 result = 0.0_rt;
573 } else {
574 result = cut_volume(mx, my, mz, alpha, 0.0_rt, 1.0_rt);
575 }
576
577 return result;
578}
579
580} // namespace kynema_sgf::multiphase
581
582#endif // VOLUME_FRACTIONS.H
kynema_sgf::constants::TIGHT_TOL
static constexpr amrex::Real TIGHT_TOL
A tight tolerance.
Definition constants.H:17
kynema_sgf::multiphase
Definition height_functions.H:8
kynema_sgf::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:456
kynema_sgf::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
kynema_sgf::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:532
kynema_sgf::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
kynema_sgf::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
kynema_sgf::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:503
kynema_sgf::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
kynema_sgf::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
kynema_sgf::InterfaceCapturingMethod::VOF
@ VOF
Volume of fluid.
Definition MultiPhase.H:26
Generated by doxygen 1.14.0