/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)) +
311 (m1 * m1 * (m1 - (3.0_rt * m3))) +
312 (m2 * m2 * (m2 - (3.0_rt * m3)))) /
313 pr;
314 } else {
315 mm = m12;
316 V3 = 0.5_rt * mm / (m3 + tiny);
317 }
318
319 amrex::Real ch = amrex::min<amrex::Real>(volF, 1.0_rt - 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_rt *
326 (m1 + std::sqrt((m1 * m1) + (8.0_rt * m2 * m3 * (ch - V1))));
327 } else if (ch < V3) {
328 p = 2.0_rt * m1 * m2;
329 q = 1.5_rt * m1 * m2 * (m12 - (2.0_rt * m3 * ch));
330 p12 = std::sqrt(amrex::max(0.0_rt, p));
331 teta = std::acos(
332 amrex::max(
333 -1.0_rt, amrex::min(1.0_rt, q / ((p * p12) + tiny)))) /
334 3.0_rt;
335 cs = std::cos(teta);
336 alpha = (p12 * (std::sqrt(3.0_rt * (1.0_rt - (cs * cs))) - cs)) + m12;
337 } else if (m12 < m3) {
338 alpha = (m3 * ch) + (0.5_rt * mm);
339 } else {
340 p = (m1 * (m2 + m3)) + (m2 * m3) - 0.25_rt;
341 q = 1.5_rt * m1 * m2 * m3 * (0.5_rt - ch);
342 p12 = std::sqrt(amrex::max(0.0_rt, p));
343 teta = std::acos(
344 amrex::max(
345 -1.0_rt, amrex::min(1.0_rt, q / ((p * p12) + tiny)))) /
346 3.0_rt;
347 cs = std::cos(teta);
348 alpha =
349 (p12 * (std::sqrt(3.0_rt * (1.0_rt - (cs * cs))) - cs)) + 0.5_rt;
350 }
351
352 if (volF > 0.5_rt) {
353 alpha = 1.0_rt - alpha;
354 }
355
356 return alpha;
357}
358
367AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real cut_volume(
368 const amrex::Real m1,
369 const amrex::Real m2,
370 const amrex::Real m3,
371 const amrex::Real alpha,
372 const amrex::Real r0,
373 const amrex::Real dr0)
374{
375 constexpr amrex::Real tiny = std::numeric_limits<amrex::Real>::epsilon();
376 amrex::Real al;
377
378 // move origin to x0
379 al = alpha - (m1 * r0);
380
381 // reflect the figure when negative coefficients
382 al = al + amrex::max<amrex::Real>(0.0_rt, -m1 * dr0) +
383 amrex::max<amrex::Real>(0.0_rt, -m2) +
384 amrex::max<amrex::Real>(0.0_rt, -m3);
385
386 // normalized equation: m1*y1 + m2*y2 + m3*y3 = alh, with 0 <= m1 <= m2 <=
387 // m3 the problem is then solved again in the unit cube
388 amrex::Real tmp = amrex::max<amrex::Real>(
389 tiny, (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 =
401 amrex::min<amrex::Real>(n1 * dr0, n2); // order coefficients
402 amrex::Real b3 = amrex::max<amrex::Real>(n1 * dr0, n2);
403 amrex::Real b2 = n3;
404
405 if (b2 < b1) {
406 tmp = b1;
407 b1 = b2;
408 b2 = tmp;
409 } else if (b2 > b3) {
410 tmp = b3;
411 b3 = b2;
412 b2 = tmp;
413 }
414
415 amrex::Real b12 = b1 + b2;
416
417 // Compute volume fraction using PLIC from Scardovelli & Zaleski (JCP 2000),
418 // adapted from code in paper by Akio Kawano (Computer & Fluids 2016)
419 amrex::Real vm1 = b1;
420 amrex::Real vm3 = b3;
421 amrex::Real vm2 = b2;
422 amrex::Real vm12 = b12;
423 amrex::Real a = al0;
424 amrex::Real v = 0.0_rt;
425
426 if (a > 0.0_rt) {
427 if (a < vm1) {
428 v = a * a * a / (6.0_rt * vm1 * vm2 * vm3);
429 } else if (a < vm2) {
430 v = (a * (a - vm1) / (2.0_rt * vm2 * vm3)) +
431 (vm1 * vm1 / ((6.0_rt * vm2 * vm3) + tiny));
432 } else if (a < amrex::min(vm12, vm3)) {
433 v = ((a * a * ((3.0_rt * vm12) - a)) +
434 (vm1 * vm1 * (vm1 - (3.0_rt * a))) +
435 (vm2 * vm2 * (vm2 - (3.0_rt * a)))) /
436 (6.0_rt * vm1 * vm2 * vm3);
437 } else if (vm3 < vm12) {
438 v = ((a * a * (3.0_rt - (2.0_rt * a))) +
439 (vm1 * vm1 * (vm1 - (3.0_rt * a))) +
440 (vm2 * vm2 * (vm2 - (3.0_rt * a))) +
441 (vm3 * vm3 * (vm3 - (3.0_rt * a)))) /
442 (6.0_rt * vm1 * vm2 * vm3);
443 } else {
444 v = (a - (0.5_rt * vm12)) / vm3;
445 }
446 }
447
448 tmp = v;
449 amrex::Real FL3D;
450 if (al <= 0.5_rt) {
451 FL3D = tmp;
452 } else {
453 FL3D = (1.0_rt - tmp);
454 }
455
456 return FL3D;
457}
458
459AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(
460 const int i,
461 const int j,
462 const int k,
463 amrex::Array4<amrex::Real const> const& volfrac,
464 amrex::Real& mx,
465 amrex::Real& my,
466 amrex::Real& mz,
467 amrex::Real& alpha)
468{
469
470 mixed_youngs_central_normal(i, j, k, volfrac, mx, my, mz);
471
472 amrex::Real invx = 1.0_rt;
473 amrex::Real invy = 1.0_rt;
474 amrex::Real invz = 1.0_rt;
475
476 if (mx < 0.0_rt) {
477 mx = -mx;
478 invx = -1.0_rt;
479 }
480 if (my < 0.0_rt) {
481 my = -my;
482 invy = -1.0_rt;
483 }
484 if (mz < 0.0_rt) {
485 mz = -mz;
486 invz = -1.0_rt;
487 }
488
489 amrex::Real mm2 = amrex::max<amrex::Real>(
490 std::numeric_limits<amrex::Real>::epsilon(), mx + my + mz);
491 mx = mx / mm2;
492 my = my / mm2;
493 mz = mz / mm2;
494
495 alpha = volume_intercept(mx, my, mz, volfrac(i, j, k));
496
497 // Back to the original plane
498 mx = invx * mx;
499 my = invy * my;
500 mz = invz * mz;
501 alpha = alpha + amrex::min<amrex::Real>(0.0_rt, mx) +
502 amrex::min<amrex::Real>(0.0_rt, my) +
503 amrex::min<amrex::Real>(0.0_rt, mz);
504}
505
506AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(
507 const int i,
508 const int j,
509 const int k,
510 const int ibdy,
511 const int jbdy,
512 const int kbdy,
513 amrex::Array4<amrex::Real const> const& volfrac,
514 amrex::Real& mx,
515 amrex::Real& my,
516 amrex::Real& mz,
517 amrex::Real& alpha)
518{
519 youngs_finite_difference_normal_neumann(
520 i, j, k, ibdy, jbdy, kbdy, volfrac, mx, my, mz);
521
522 amrex::Real norm = std::abs(mx) + std::abs(my) + std::abs(mz);
523 if (norm <= std::numeric_limits<amrex::Real>::epsilon()) {
524 mx = -1.0_rt;
525 my = 0.0_rt;
526 mz = 0.0_rt;
527 norm = 1.0_rt;
528 }
529 mx /= norm;
530 my /= norm;
531 mz /= norm;
532 alpha = volume_intercept(
533 std::abs(mx), std::abs(my), std::abs(mz), volfrac(i, j, k));
534 alpha += amrex::min(0.0_rt, mx) + amrex::min(0.0_rt, my) +
535 amrex::min(0.0_rt, mz);
536}
537
538AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
539 const int i,
540 const int j,
541 const int k,
542 amrex::Array4<amrex::Real const> const& volfrac,
543 const int n_band = 1,
544 const amrex::Real vof_tol = constants::VOF_TOL)
545{
546 // n_band must be <= number of vof ghost cells (3)
547
548 // Check if near interface
549 amrex::Real VOF_max = 0.0_rt;
550 amrex::Real VOF_min = 1.0_rt;
551 bool VOF_mid = false;
552 for (int ii = -n_band; ii <= n_band; ++ii) {
553 for (int jj = -n_band; jj <= n_band; ++jj) {
554 for (int kk = -n_band; kk <= n_band; ++kk) {
555 amrex::Real VOF = volfrac(i + ii, j + jj, k + kk);
556 VOF_max = amrex::max<amrex::Real>(VOF_max, VOF);
557 VOF_min = amrex::min<amrex::Real>(VOF_min, VOF);
558 if (VOF < 1.0_rt - vof_tol && VOF > vof_tol) {
559 VOF_mid = true;
560 }
561 }
562 }
563 }
564 return (std::abs(VOF_max - VOF_min) > vof_tol || VOF_mid);
565}
566
567AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real levelset_to_vof(
568 const int i,
569 const int j,
570 const int k,
571 const amrex::Real eps,
572 amrex::Array4<amrex::Real const> const& phi)
573{
574 amrex::Real mx, my, mz;
575 youngs_finite_difference_normal(i, j, k, phi, mx, my, mz);
576 mx = std::abs(mx / 32.0_rt);
577 my = std::abs(my / 32.0_rt);
578 mz = std::abs(mz / 32.0_rt);
579 const amrex::Real norm_raw = (mx + my + mz);
580 if (norm_raw <= std::numeric_limits<amrex::Real>::epsilon()) {
581 const amrex::Real phi_ijk = phi(i, j, k);
582 if (phi_ijk < -eps) {
583 return 0.0_rt;
584 }
585 if (phi_ijk > eps) {
586 return 1.0_rt;
587 }
588 return 0.5_rt;
589 }
590 amrex::Real normL1 = norm_raw;
591 mx = mx / normL1;
592 my = my / normL1;
593 mz = mz / normL1;
594 // Make sure that alpha is negative far away from the
595 // interface
596 amrex::Real alpha;
597 if (phi(i, j, k) < -eps) {
598 alpha = -1.0_rt;
599 } else {
600 alpha = phi(i, j, k) / normL1;
601 alpha = alpha + 0.5_rt;
602 }
603 amrex::Real result;
604 if (alpha >= 1.0_rt) {
605 result = 1.0_rt;
606 } else if (alpha <= 0.0_rt) {
607 result = 0.0_rt;
608 } else {
609 result = cut_volume(mx, my, mz, alpha, 0.0_rt, 1.0_rt);
610 }
611
612 return result;
613}
614
615} // namespace kynema_sgf::multiphase
616
617#endif // VOLUME_FRACTIONS.H
kynema_sgf::constants::VOF_TOL
static constexpr amrex::Real VOF_TOL
Tolerance for vof determining presence of liquid.
Definition constants.H:34
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:459
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:567
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:367
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 vof_tol=constants::VOF_TOL)
Definition volume_fractions.H:538
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