flow 1.3.0
Kokkos cut-cell IBM incompressible Navier-Stokes solver + pnm pore extraction
Loading...
Searching...
No Matches
interface_area.hpp File Reference

flow — WO-P3c: the interfacial AREA taken from the CURVATURE cascade's geometry. More...

#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
#include "vof/curvature.hpp"
#include "vof/phase_change.hpp"
#include "vof/plic.hpp"

Go to the source code of this file.

Namespaces

namespace  peclet
 
namespace  peclet::flow
 
namespace  peclet::flow::vof
 

Enumerations

enum  peclet::flow::vof::InterfaceAreaMode : int {
  peclet::flow::vof::kAreaPlic = 0 , peclet::flow::vof::kAreaMetric = 1 , peclet::flow::vof::kAreaNormal = 2 , peclet::flow::vof::kAreaFootprint = 3 ,
  peclet::flow::vof::kAreaMcColour = 4 , peclet::flow::vof::kAreaMcColourSplit = 5 , peclet::flow::vof::kAreaMcPlic = 6 , peclet::flow::vof::kAreaMcPlicSplit = 7
}
 Which construction turns the cascade's normal into an area. More...
 

Functions

double peclet::flow::vof::hfAreaElement (const double h[9])
 The height function's area element sqrt(1 + h_x^2 + h_y^2), with the same central differences hfPatchKappa uses — so the metric and the curvature are two derivatives of one interpolant, not two independent estimates.
 
double peclet::flow::vof::hfAreaElement (const double h[9], const VofMetric &g, int d)
 The same on an ANISOTROPIC cell: the graph is the PHYSICAL f = h_d * h_idx over the physical transverse coordinates, so the two slopes carry the ratios of hfPatchKappa(h, g, d) (V2.3).
 
double peclet::flow::vof::hfFootprintArea (const double h[9])
 Variant C — the height function's own footprint times its own metric (kAreaFootprint).
 
double peclet::flow::vof::hfFootprintArea (const double h[9], const VofMetric &g, int d)
 The same on an ANISOTROPIC cell.
 
void peclet::flow::vof::hfSurfaceNormal (const double h[9], int d, double m[3])
 The unit surface normal implied by the SAME 3x3 patch of heights, in world components.
 
void peclet::flow::vof::hfSurfaceNormal (const double h[9], int d, const VofMetric &g, double m[3])
 The same on an ANISOTROPIC cell — the PHYSICAL unit normal of the graph, i.e.
 
void peclet::flow::vof::pvSurfaceNormal (const double a[6], const double t1[3], const double t2[3], const double nn[3], double m[3])
 The unit normal of the fitted paraboloid z = a0 + a1 x + a2 y + ... AT THE ORIGIN of its own frame (the target cell's PLIC centroid), in world components: (-a1, -a2, 1) rotated by (t1, t2, nn) and normalized.
 
double peclet::flow::vof::interfaceAreaFromNormal (int mode, double mx, double my, double mz, double alpha, double c, const double ns[3], const VofMetric &g)
 The projection axis of kAreaMetric is the one ns is most aligned with — the column direction on tiers 1/2 by construction, and the same choice on tier 3 — because the footprint factor |n.e_d| is the one that stays farthest from zero there.
 
double peclet::flow::vof::interfaceAreaFromNormal (int mode, double mx, double my, double mz, double alpha, double c, const double ns[3])
 Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.
 

Detailed Description

flow — WO-P3c: the interfacial AREA taken from the CURVATURE cascade's geometry.

Container-free by contract, exactly like plic.hpp and phase_change.hpp: scalars and small local arrays only. The view-level driver (the block walk, the height columns, the PV fit) is vof/interface_area_field.hpp, the same split curvature.hpp / curvature_field.hpp already uses.

Why this file exists (the measurement, WO-P3b)

Rung P0/P1 takes the interfacial area of a cell from the PLIC polygon of the MYC normal, A = |m|_2 dV/dalpha (vof::plicArea). On a PLANE that is exact to the last bit — which is every gate P0a/P0b/P1/P2 runs on. On a SPHERE it is 5.5-9.3 % LOW at R = 4...28 and does not converge (a-priori probe on exact fractions; marching cubes on the same colour field returns 4 pi R^2 to 0.4 %). plicArea is linear in the normal, and rung V0 measured MYC's normal error as non-convergent (order 0.83, ~28 % of mixed cells taking the Youngs branch), so a non-convergent normal gives a non-convergent area. That deficit is exactly the P3 Scriven growth-rate miss: the flux PER UNIT AREA is within +-0.5 % of Scriven's exact mdot while the bubble grows as int mdot dA, and beta_eff sits 2.5 % low at every mesh.

The cure is to take A_Gamma from the SAME geometry the curvature comes from — the height function on tiers 1/2 and the PLIC-volumetric paraboloid on tier 3 (core/vof/curvature.hpp). Both are built from column sums / polygon volumes and never differentiate a normal, which is precisely why the V3 cascade converges where MYC does not. It is also the consistent choice: kappa and A_Gamma then describe one surface.

The two constructions, and what each one assumes

Write the interface locally as the graph x_d = f(x_d1, x_d2) over the column direction d. The area of the piece inside one cell is int_R sqrt(1 + f_x^2 + f_y^2) over that piece's FOOTPRINT R in the transverse plane, and the footprints of the cells of one column tile the unit square exactly. So an area needs two things: a metric (a slope) and a footprint.

  • **kAreaMetric (variant A, the default).** Keep the PLIC polygon's own FOOTPRINT — the projection of the reconstructed polygon on the transverse plane, A_PLIC |n_d|, which is what makes the cells of a column tile — and replace only the SLOPE:
    A = A_PLIC(m_MYC, alpha) * |n_MYC . e_d| * sqrt(1 + h_x^2 + h_y^2) .
    
    The metric is second order on a sphere (it is the height function's own, the same central differences hfPatchKappa differentiates once more), and the footprint carries the exact cell volume through alpha. Equivalently: A = A_PLIC * |n_MYC.e_d| / |n*.e_d|, i.e. the MYC area with its slope factor divided out and the accurate one put back.
  • **kAreaNormal (variant B).** Rebuild the whole plane with the accurate normal and the cell's own colour: A = plicArea(n*, plicAlpha(n*, C)). No footprint bookkeeping at all — the plane shift does it — but the polygon is the one a plane with THAT normal cuts, which on a curved interface is not the one the transport actually carries.

Both are exact on a plane (there h_x, h_y are exact and n* = n_MYC), which is what keeps rungs P0-P2 where they were.

Definition in file interface_area.hpp.