Unfitted discontinuous Galerkin method for surface hyperbolic PDEs

Maxime Bouyges, Ghislain Blanchard, Lokman Bennani

ONERA - Multiphysics Department for Energetics

Foreword : icing simulation recipe

Four steps needed to perform an icing simulation around an obstacle:

  • compute the aerodynamic flow-field (body-fitted or Immersed Boundary Method (IBM) for instance)
  • compute droplets trajectories (body-fitted or IBM for instance)
  • compute a surface equation to obtain the ice height and/or the liquid film height (body-fitted)
  • let the ice grow (re-mesh if body-fitted)

Streamlines of the air around an iced airfoil (in red) using IBM

Droplets trajectories around the iced airfoil (in red) using IBM

Would it be possible the solve the surface equation without re-meshing, with an IBM-type method?

Back on surface equations

Consider the following advection-diffusion equation on a surface \(\Gamma \in \Omega\) \[ \partial_t q + \nabla_\Gamma \cdot F = S_\Gamma \]

Tangential gradient definition \[ \begin{align} \nabla_\Gamma u & = P \nabla u \\ & = \nabla u - (\nabla u \cdot \nu) \nu\\ \end{align} \]

One possible definition for the normal vector \[ \nu = \nabla d \] where \(d\) is the signed distance to \(\Gamma\)

The brilliant idea (Bertalmío 2001)

The distance to \(\Gamma\) (level-set) is defined in the entire1 volume \(\Omega\), hence we can define \(\nu_d = \nabla d\) everywhere in \(\Omega\)…

We can define a tangential projector \(P_d = I - \nu_d \otimes \nu_d\) everywhere in \(\Omega\)!

We can define a tangential gradient2 \(\nabla_d u = P_d\nabla u\) everywhere in \(\Omega\)!!

Interpretation

What’s the meaning of the solution of the following (volume) heat equation? \[ \partial_t T - \nabla_d \cdot \left( k \nabla_d T \right) = 0 \]

on each isocontour \(d = cte\), the solution is identical to the (surface) solution of the surface heat equation associated to this isocontour

the \(\Delta_d\) operator prevents any diffusion in the normal direction of the isocontour

The volume PDE can be solved with any classical numerical method (FDM, FVM, FEM, DG etc)…


Caution with the weak form, a special “integration by parts” formula must be applied: \[ \int_\Omega \varphi \color{red}{\nabla_d} \cdot F \, \mathrm{d}\Omega = \int_{\partial \Omega} \varphi F \cdot \color{red}{P_d}\nu_{\partial \Omega} \,\mathrm{d}S - \int_\Omega F \cdot ( \color{red}{\nabla_d} \varphi - \color{red}{\varphi}\underbrace{\color{red}{H}}_{=\nabla \cdot \nu_d} \color{red}{\nu_d} ) \, \mathrm{d} \Omega \]

Surface displacement and deformation

Taking into account surface displacement and/or deformation over time is easy because the mesh is fixed. Let \(v = -\partial_t d\) the displacement-deformation velocity, the previous heat equation reads \[ \partial_t T + \underbrace{v \cdot \nabla T}_{\text{displacement}} + T \underbrace{\nabla_d \cdot v}_{\text{expansion}} - \nabla_d \cdot (k\nabla_d T) = 0 \]

Heat equation around a growing circle in movement: \(\tiny{d(x,t) = |x - (\underbrace{x_0 + v_{trans}t}_{\substack{\text{center}\\ \text{translation}}})| - (R_0 + \underbrace{\alpha t}_{\text{expansion}})}\)

Can we solve hyperbolic equations as well?

Hyperbolic system

Consider the following shallow water non-linear hyperbolic system: \[ \begin{cases} & \partial_t h + \nabla_\Gamma \cdot hu = 0\\ & \partial_t hu + \nabla_\Gamma \cdot (hu\otimes u + 0.5gh^2 \mathcal{I}) = 0 \end{cases} \]

Apply a discontinuous Galerkin discretization (in \(\Omega\)):

\[ \begin{align} & \forall \varphi_h,~\int_\Omega \varphi_h \partial_t h - \int_\Omega \mathcal{F}_h \cdot (\nabla_d \varphi_h - \varphi_h H \nu_d ) + \int_{\partial \Omega} [\varphi_h] \mathcal{F}_h^* \cdot (P_d\nu_{\partial \Omega})=0\\ & \forall \varphi_{hu},~\int_\Omega \varphi_{hu} \cdot \partial_t hu - \int_\Omega \mathcal{F}_{hu} : (\nabla_d \varphi_{hu} - H\varphi_{hu} \otimes \nu_d) + \int_{\partial \Omega} \mathcal{F}_{hu}^* : ([\varphi_{hu}] \otimes (P_d\nu_{\partial \Omega})) =0 \end{align} \]

\[ \begin{align} & \forall \varphi_h,~\int_\Omega \varphi_h \partial_t h - \int_\Omega \mathcal{F}_h \cdot (\color{red}{\nabla_d} \varphi_h - \color{red}{\varphi_h H \nu_d} ) + \int_{\partial \Omega} [\varphi_h] \mathcal{F}_h^* \cdot (\color{red}{P_d}\nu_{\partial \Omega})=0\\ & \forall \varphi_{hu},~\int_\Omega \varphi_{hu} \cdot \partial_t hu - \int_\Omega \mathcal{F}_{hu} : (\color{red}{\nabla_d} \varphi_{hu} - \color{red}{H\varphi_{hu} \otimes \nu_d}) + \int_{\partial \Omega} \mathcal{F}_{hu}^* : ([\varphi_{hu}] \otimes (\color{red}{P_d}\nu_{\partial \Omega})) =0 \end{align} \]


Spot the 6 Differences : what differences with the usual weak form of these equations?

Dambreak… around a circle

We study the evolution of the previous shallow-water equations around a circle with the following initial condition: \[ h(\theta, t=0) = \begin{cases} 3\text{ if }\theta \leq \pi/3\\2\text{ otherwise}\end{cases},~~~~~u(\theta,t=0) = 0 \]

Promising! But observe how the discontinuities progressively differ from circle radii…

Origin of the problem

Consider the following scalar transport equation with constant velocity \(c\) \[ \partial_t q + c \cdot \nabla_\Gamma q = 0 \]

Around the unit circle, we obtain \[ \partial_t q + \color{#1154a6}{c} \partial_\theta q = 0 \] but the use of \(\nabla_d = P_d \nabla\) leads to (in the volume) \[ \partial_t q + \color{#1154a6}{\dfrac{c}{r}} \partial_\theta q = 0 \]

\(q(r,\theta,t=0) = sin(\theta)\)

On each isocontour \(d = R\), a transport equation with velocity \(c/R\) is solved. This is harmless on the continuous level because of the absence of normal advection-diffusion, but on a discrete level …

Solution of the problem

Initially proposed by Greer (2006). Let \(x \in \Omega\), we note \(p(x)\in\Gamma\) the orthogonal projection of \(x\) on \(\Gamma\) and \(\tilde{u} = u(p(x))\). We can show that \[ \underbrace{(\mathcal{I} - d \mathcal{H})^{-1}P_d}_{\tilde{P}_d} \nabla \tilde{u} = \left(\nabla_\Gamma u\right) \circ p \] \(\tilde{P}_d\) enables the evaluation of a gradient on the \(d=0\) contour, not “locally”!

With \(\nabla_d = \tilde{P}_d\nabla\), the previous transport equation in \(\Omega\) reads \[ \partial_t q + \color{#1154a6}{c} \partial_\theta q = 0 \]

\(q(r,\theta,t=0) = sin(\theta)\)

Corrected dambreak around a circle

With the Greer tangential gradient, the correct solution is obtained on quads and/or triangles, at first and second order.

Unfitted resolution of the hyperbolic shallow-water equation around a circle with \(Q_0\) elements

And it converges!

Convergence has been checked for \(Q_0\), \(Q_1\), \(P_0\) and \(P_1\) Lagrange elements for both the scalar transport equation and the shallow-water equations around a circle. The error is found to be independant of the tubular region radius.

\(L^2\) error for \(Q_1\) Lagrange elements for the shallow-water equations around a circle

A non-constant curvature case

Dambreak around a star, gravity always normal to the surface.

Left: volumic solution, right: comparison between unfitted and fitted solutions

A demo in 3D: the Thacker problem (and solution)

\(\Rightarrow\) Periodic oscillations of a plane surface in a paraboloid bassin. Below : “extrusion” of the surface height along the normal.

Comparison between the unfitted approach (left) and the exact solution (right)

Conclusions and perspectives

Wrap-up:

  • it is possible to solve surface hyperbolic equations using an implicit representation of the surface, ie without meshing the surface
  • 1st and 2nd order convergence has been demonstrated

Coming next:

  • how should the geometrical operators \(d\), \(\nu_d\), \(\tilde{P}_d\), \(\mathcal{H}_d\), etc be discretized when not known anatically?
  • couple a surface hyperbolic equation with a (true) volume equation, for instance shallow-water equations over a heated body
  • explore other similar methods, for instance the closest-point method

Thanks, any question?

All the codes have been written using the Bcube Julia library.

Bcube is an open-source FEM-DG discretization library, developed at ONERA. It is written in Julia.

Helmholtz eqn Phase field solidification Linear transport eqn
Heat eqn on a sphere Transport eqn on hypersurfaces Linear thermo-elasticity

Appendix

A few references

1 M. Bertalmıo, L.-T. Cheng, S. Osher, and G. Sapiro, “Variational problems and partial differential equations on implicit surfaces,” Journal of Computational Physics, (2001).
2 G. Dziuk, and C.M. Elliott, “Finite element methods for surface PDEs,” Acta Numerica, (2013).
3 J.B. Greer, “An improvement of a recent eulerian method for solving PDEs on general geometries,” Journal of Scientific Computing, (2006).
4 K. Deckelnick, G. Dziuk, C.M. Elliott, and C.-J. Heine, “An h-narrow band finite-element method for elliptic equations on implicit surfaces,” IMA Journal of Numerical Analysis, (2010).
5 S.J. Ruuth, and B. Merriman, “A simple embedding method for solving partial differential equations on surfaces,” Journal of Computational Physics, (2008).

Boundary conditions

At the embedding domain limits, the value applied in the ghost cell is obtained from the corresponding value on the iso-0 contour with orthogonal projection : \(x = p(x) + d \nu(x)\)

Integration by parts with Greer tangential gradient

Recall that \(\nabla_d = \tilde{P}_d\nabla\) with \(\tilde{P}_d = (\mathcal{I}-d\mathcal{H})^{-1}P\) where \(\mathcal{H}_{ij} = \dfrac{\partial^2 d}{ \partial x_i \partial x_j}\).

Observe that \(H = \mathrm{trace}(\mathcal{H})= \nabla \cdot \nu\) and \(\mathcal{H}^0(x) = \mathcal{H}(p(x))\).

Then, the “\(d\)-integration by parts” reads \[ \int_\Omega (\nabla_d \cdot \mathcal{F}) \varphi \, \mathrm{d}x = \int_{\partial \Omega} \varphi \mathcal{F} \cdot \tilde{P} \nu_{\partial \Omega} \, \mathrm{d}x - \int_\Omega \mathcal{F} \cdot \left[ \nabla_d \varphi - \left(H \nu - \nabla \cdot \mathcal{H}^0\right)\right] \, \mathrm{d}x \]

Errors \(L^1\), \(L^2\) and \(L^\infty\)

Errors for the shallow-water equations around a circle with \(Q_1\) Lagrange elements.

Code example

module main
using Bcube
using LinearAlgebra
using StaticArrays
using DelimitedFiles
include(joinpath(@__DIR__, "..", "common", "common.jl"))

const eps_h = 1.0e-10
const hmin₀ = 1.e-8
const hmax₀ = 1.0e10
const DMPcurv₀ = 0.0
const wall_friction = false
@warn "wall friction is $(wall_friction)"

velocity(h, hu) = (hu * 2 * h) / (h * h + max(h * h, (3e-5)^2))  #desingularization

function _flux_HLL(qL, qR, n, flux, f_λ)
    λL, λR = f_λ(qL), f_λ(qR)
    λ⁻ = min(minimum(λL), minimum(λR), zero(λL[1]))
    λ⁺ = max(maximum(λL), maximum(λR), zero(λL[1]))
    function f_HLL(qL, qR, fL, fR)
        if abs(λ⁺ - λ⁻) > 1.0e-12
            fLn, fRn = dotn(fL, n), dotn(fR, n)
            f = (λ⁺ * fLn - λ⁻ * fRn + λ⁻ * λ⁺ * (qR - qL)) / (λ⁺ - λ⁻)
        else
            f = 0.5 * (fL(qL) + fR(qR))
        end
        return f
    end
    map(f_HLL, qL, qR, flux(qL), flux(qR))
end

function shallow_water_maxeigval(q, gravity)
    h, hu = q
    u = velocity(h, hu)
    return norm(u) + √(norm(gravity) * max(h, eps_h))
end

function shallow_water_maxeigval(q, n, gravity)
    h, hu = q
    un = velocity(h, hu) ⋅ n
    return abs(un) + √(norm(gravity) * max(h, eps_h))
end

function shallow_water_eigval(q, n, gravity)
    h, hu = q
    un = velocity(h, hu) ⋅ n
    c = √(norm(gravity) * max(h, eps_h))
    return un - c, un + c
end

function compute_timestep!(q, mesh, dimcar, gn, CFL)
    h, hu = q
    degree = minimum(
        feSpace -> get_degree(Bcube.get_function_space(feSpace)),
        get_fespace.(get_fe_functions(q)),
    )

    λmax = var_on_centers(norm ∘ velocity(h, hu) + √(gn * h), mesh)
    Δt = CFL * minimum(dimcar ./ λmax)
    return Δt / (2degree + 1)
end

function flux_fitted(qi, qj, gni, gnj, Ri, Rj, φi, φj, nij, flux)
    hj, huj = qj
    _qj = (hj, Ri * huj)

    φ_hi, φ_hui = φi
    φ_hj, φ_huj = φj
    δφ = (φ_hi - φ_hj, φ_hui - transpose(Rj) * φ_huj)

    return flux(qi, _qj, gni, δφ, nij)
end

function flux_unfitted(qi, qj, gni, gnj, Pi, Pj, ℋi, ℋj, ϕi, ϕj, δφ, nij, flux)
    _nij = inv(I - ϕi * ℋi) * Pi * nij
    return flux(qi, qj, gni, δφ, _nij)
end

function flux_HLL(qi, qj, gni, δφ, nij)
    δv_h, δv_hu = δφ

    f_λ = x -> shallow_water_eigval(x, nij, gni)
    flux = _flux_HLL(qi, qj, nij, y -> flux_sw(y, gni), f_λ)

    flux_h, flux_hu = flux
    return flux_h ⋅ δv_h + flux_hu ⋅ δv_hu
end

function _flux_HLL(qL, qR, n, flux, f_λ)
    λL, λR = f_λ(qL), f_λ(qR)
    λ⁻ = min(minimum(λL), minimum(λR), zero(λL[1]))
    λ⁺ = max(maximum(λL), maximum(λR), zero(λL[1]))
    function f_HLL(qL, qR, fL, fR)
        fLn, fRn = dotn(fL, n), dotn(fR, n)
        if abs(λ⁺ - λ⁻) > 1.0e-12
            f = (λ⁺ * fLn - λ⁻ * fRn + λ⁻ * λ⁺ * (qR - qL)) / (λ⁺ - λ⁻)
        else
            f = 0.5 * (fLn + fRn)
        end
        return f
    end
    map(f_HLL, qL, qR, flux(qL), flux(qR))
end

function flux_sw(q, gn)
    h, hu = q
    u = velocity(h, hu)
    huu = hu * transpose(u)
    p_grav = 0.5 * gn * h * h
    return h .* u, huu + p_grav * I
end

function apply_limitation!(q, params, cache)
    h, hu = q
    dΩ = params.dΩ

    q_mean = Bcube.cell_mean(q, cache.cacheCellMean)

    lim_h, h_proj = linear_scaling_limiter(
        h,
        dΩ;
        bounds=(hmin₀, hmax₀),
        DMPrelax=params.DMPrelax,
        mass=cache.mass_sca,
    )
    set_dof_values!(h, get_dof_values(h_proj))

    _, hu_mean, = q_mean
    limited_var(a, a̅, lim_a) = a̅ + lim_a * (a - a̅)
    projection_l2!(hu, limited_var(hu, hu_mean, lim_h), dΩ; mass=cache.mass_vec)
end

"""
Define closest point interpolator of a function f : here we cheat because we now the closest point on Gamma,
which is assumed to be a circle of radius 1 centered in (0,0)
"""
function closest_point_interp_func(pf, f)
    return x -> begin
        θ = atan(x[2], x[1])
        Bcube.interpolate_at_point(pf, SA[cos(θ), sin(θ)], f...)
    end
end

closest_point_interp = PhysicalFunction ∘ closest_point_interp_func

"""
Warning : only valid for P1 geometrical elements
"""
function fitted_closest_point_interp_func(pf, f)
    return x -> begin
        icell = Bcube.find_cell_index(pf, x)
        cinfo = Bcube.CellInfo(pf.mesh, icell)
        ctype = Bcube.celltype(cinfo)
        @assert ctype isa Bcube.Bar2_t "not valid for elements different from P1"
        A, B = Bcube.nodes(cinfo)
        u = normalize(B.x - A.x)
        l = (x - A.x) ⋅ u
        H = A.x + l * u
        y = Bcube.CellPoint(H, cinfo, Bcube.PhysicalDomain())
        f_icell = Bcube.materialize(f, cinfo)
        return Bcube.materialize(f_icell, y)
    end
end

fitted_closest_point_interp = PhysicalFunction ∘ fitted_closest_point_interp_func


function unfitted_dg(;
    mesh,
    _ϕ,
    h0,
    u0,
    g,
    μ,
    tfinal,
    degree,
    CFL,
    flux,
    Δt_min=0.0,
    nitemax,
    use_constant_Δt=false, # if true, a constant Δt is used
    constant_Δt=0.0,
    limitation=degree > 0,
)
    (degree > 0 && !limitation) && @warn "degree > 0 but limitation disabled"

    # Mesh
    quad = Quadrature(QuadratureLobatto(), 2 * degree + 1)
    dΩ = Measure(CellDomain(mesh), quad)
    dΓ = Measure(InteriorFaceDomain(mesh), quad)
    dΛ = Measure(BoundaryFaceDomain(mesh), quad)
    nΓ = get_face_normals(dΓ)
    nΛ = get_face_normals(dΛ)
    dimcar = compute_dimcar(mesh)

    # Projection material
    ϕ = PhysicalFunction(_ϕ)
    ν = PhysicalFunction(x -> x ./ (norm(x) + 1e-20)) # analytic
    H = PhysicalFunction(x -> 1 / (norm(x) + 1e-20)) # analytic
    ℋ = PhysicalFunction(
        x -> begin
            x2 = x[1]^2
            y2 = x[2]^2
            return SA[
                y2 (-x[1]*x[2])
                (-x[1]*x[2]) x2
            ] ./ (x2 + y2)^(3 / 2)
        end,
    ) # analytic
    divℋ0 = PhysicalFunction(x -> -x ./ (norm(x) + 1e-20))
    P = I - (ν ⊗ ν)

    # FESpace
    U_h = TrialFESpace(FunctionSpace(:Lagrange, degree), mesh, :discontinuous)
    U_hu = TrialFESpace(
        FunctionSpace(:Lagrange, degree),
        mesh,
        :discontinuous;
        size=Bcube.spacedim(mesh),
    )
    V_h = TestFESpace(U_h)
    V_hu = TestFESpace(U_hu)
    U = MultiFESpace(U_h, U_hu)
    V = MultiFESpace(V_h, V_hu)

    # FEFunction
    q = FEFunction(U)
    projection_l2!(q, (h0, u0), mesh)

    # Gravity op
    gn = abs(g ⋅ ν)

    # Flux
    function flux_Γ(q, gn, P, ℋ, ϕ, φ, n)
        flux_unfitted ∘ (
            side⁻(q),
            side⁺(q),
            side⁻(gn),
            side⁺(gn),
            side⁻(P),
            side⁺(P),
            side⁻(ℋ),
            side⁺(ℋ),
            side⁻(ϕ),
            side⁺(ϕ),
            jump(φ),
            side⁻(n),
            flux,
        )
    end

    function _flux_Λ(qi, closest_q, gni, Pi, ℋi, ϕi, φi, nij)
        # Closest-point version
        qg = closest_q
        return flux_unfitted(qi, qg, gni, gni, Pi, Pi, ℋi, ℋi, ϕi, ϕi, φi, nij, flux)
    end
    function flux_Λ(q, closest_q, gn, P, ℋ, ϕ, φ, n)
        _flux_Λ ∘ (
            side⁻(q),
            side⁻(closest_q),
            side⁻(gn),
            side⁻(P),
            side⁻(ℋ),
            side⁻(ϕ),
            side⁻(φ),
            side⁻(n),
        )
    end

    function _flux_Ω(q, φ, ∇φ, gn, ν, P, ℋ, H, divℋ0, ϕ)
        h, hu = q
        u = hu / h
        φh, φhu = φ
        ∇φh, ∇φhu = ∇φ

        _q = (h, P * hu)
        f_h, f_hu = flux_sw(_q, gn)

        # Version Greer
        Pg = inv(I - ϕ * ℋ) * P
        ∇Γ_φh = Pg * ∇φh # scalar version
        ∇Γ_φhu = ∇φhu * Pg # vector version
        expr =
            f_h ⋅ (∇Γ_φh - φh * (H * ν - ϕ * divℋ0)) +
            f_hu ⊡ (∇Γ_φhu - φhu ⊗ (H * ν - ϕ * divℋ0))

        if wall_friction
            expr += -3μ * u / h ⋅ φhu
        end

        return expr
    end
    function flux_Ω(q, gn, P, ℋ, H, divℋ0, ϕ, φ)
        _flux_Ω ∘ (q, φ, map(∇, φ), gn, ν, P, ℋ, H, divℋ0, ϕ)
    end

    # Rhs function
    function compute_rhs!(rhs, qdofs)
        _q = (FEFunction(U, qdofs)...,)
        _closest_q = closest_point_interp(pf, _q)
        function l(φ)
            ∫(flux_Ω(_q, gn, P, ℋ, H, divℋ0, ϕ, φ))dΩ -
            ∫(flux_Γ(_q, gn, P, ℋ, ϕ, φ, nΓ))dΓ -
            ∫(flux_Λ(_q, _closest_q, gn, P, ℋ, ϕ, φ, nΛ))dΛ
        end
        assemble_linear!(rhs, l, V)
    end

    function compute_rhs(qdofs)
        rhs = zero(qdofs)
        compute_rhs!(rhs, qdofs)
        return rhs
    end

    # Mass matrix
    m(q, φ) = ∫(q ⋅ φ)Measure(CellDomain(mesh), 2 * degree + 1) # whole domain
    _M = assemble_bilinear(m, U, V)
    M = factorize(_M)

    # Limitation
    if degree > 0
        DMPrelax = DMPcurv₀ .* dimcar .^ 2
        params = (; dΩ, DMPrelax)
        cache = (
            mass_sca=Bcube.build_mass_matrix(U_h, V_h, dΩ),
            mass_vec=Bcube.build_mass_matrix(U_hu, V_hu, dΩ),
            cacheCellMean=Bcube.build_cell_mean_cache(q, dΩ),
        )
        apply_limitation!(q, params, cache)
    end

    # Loop
    t = 0.0
    last_ite = false
    rhs = Bcube.allocate_dofs(U)
    qdofs = copy(get_dof_values(q))
    for ite in 1:nitemax
        Δt = if !use_constant_Δt
            compute_timestep!(FEFunction(U, qdofs), mesh, dimcar, gn, CFL)
        else
            constant_Δt
        end
        Δt = max(Δt, Δt_min)
        if Δt > tfinal - t
            last_ite = true
            Δt = tfinal - t
        end

        qdofs .= get_dof_values(q)
        # udofs .= forward_euler(udofs, x -> M \ compute_rhs(x), Δt)
        qdofs .= rk3_ssp(qdofs, x -> M \ compute_rhs(x), Δt)

        t += Δt

        set_dof_values!(q, qdofs)
        limitation && apply_limitation!(q, params, cache)

        last_ite && break
    end

    return (; q, pf)
end

function run_one_case()
    # Settings
    R = 1 # radius
    CFL = 0.5
    tfinal = 0.5
    degree = 0
    outdir = joinpath(@__DIR__, "tmp")
    gn = 1.0
    ϕmax = 0.3
    n = 1024

    mkpath(outdir)

    # Initial solution / conditions
    μ = 5.0
    θc = 0.0
    θlr = π / 3
    hl0 = 3.0
    hr0 = 2.0
    _θ(x) = atan(x[2], x[1])
    u0 = PhysicalFunction(x -> @SVector zeros(Bcube.spacedim(mesh_fitted)))
    # h0 = PhysicalFunction(x -> abs(_θ(x) - θc) ≤ θlr ? hl0 : hr0)
    h0 = PhysicalFunction(
        x -> begin
            __θ = _θ(x) - θc
            if abs(__θ) ≤ θlr
                return hr0 + (hl0 - hr0) * exp(1 / θlr^2) * exp(-1 / (θlr^2 - __θ^2)) # wrong
            else
                return hr0
            end
        end,
    )
    function _g(x)
        if norm(x - xc) < 1e-9
            return SA[0.0, 0.0] # no gravity vector on center
        else
            return gn * normalize(x - xc)
        end
    end
    g = PhysicalFunction(_g)

    # Unfitted solution
    ϕ(x) = norm(x) - R # signed distance
    l = 3R

    ## Mesh 1
    println("UNFitted DG, n = $n, ϕmax = $ϕmax, degree = $degree")
    mesh_unfitted =
        rectangle_mesh(n, n; xmin=-l / 2, xmax=l / 2, ymin=-l / 2, ymax=l / 2)
    indices = Bcube.identify_cells(mesh_unfitted, x -> abs(ϕ(x)) < ϕmax)
    @assert length(indices) > 0 "no cell in clipped mesh" # the ϕmax is too harsh and there is no cell in the clipped mesh
    mesh_unfitted = Bcube.domain_to_mesh(CellDomain(mesh_unfitted, indices))
    if ncells(mesh_unfitted) <= 21 # check that all cells have at least one neighbor cell
        c2c = Bcube.connectivity_cell2cell_by_faces(mesh_unfitted)
        @assert c2c.minsize > 0 "disjoined cells"
    end

    sol_unfitted = unfitted_dg(;
        mesh=mesh_unfitted,
        _ϕ=ϕ,
        h0,
        u0,
        g,
        μ,
        tfinal,
        degree,
        flux=flux_HLL,
        CFL,
        nitemax=typemax(Int),
        Δt_min=0.0,
        limitation=false,
        use_constant_Δt=true,
        constant_Δt=CFL * l / (n - 1) / √(gn * max(hl0, hr0)) / (2degree + 1),
    )
end

run_one_case()

end