Level class

The Level class solves the diffusion-transport equation:

\[ \begin{aligned} \partial_t \phi+\theta\cdot\nabla\phi &= h^2\Delta\varphi && x\in \mathcal{D},\, t>0, \\ \partial_{n} \phi &= 0 && x\in \partial\mathcal{D},\, t>0, \\ \phi(0, x) &= \phi^{i}(x) && x\in \mathcal{D}, \end{aligned} \]

where \(\theta\) is the velocity field, \(h\) the mesh diameter, and \(\phi^i\) the level set function at iteration \(i\). A time step \(\Delta t\) and a final time \(T\) must be provided. We choose the next level set function as \(\phi^{i+1}(x) = \phi(T,x)\).

---

Details

This class implements the Petrov Galerkin and Crank-Nicolson methods to solve the transport equation corresponding to the level set function.

Attributes:

Name	Type	Description
dt	Constant	Time step.
domain	Mesh	Problem domain.
a	Form	Bilinear part of the iterative scheme.
L	Form	Linear part of the iterative scheme.

Methods:

Name	Description
run	Run the iterative method over a fixed number of time steps.

Source code in code/formopt.py

class Level:
    """
    This class implements the Petrov Galerkin and
    Crank-Nicolson methods to solve
    the transport equation corresponding
    to the level set function.

    Attributes
    ----------
    dt : Constant
        Time step.
    domain : Mesh
        Problem domain.
    a : Form
        Bilinear part of the iterative scheme.
    L : Form
        Linear part of the iterative scheme.

    Methods
    -------
    run(phi: Function, steps: int, tend: float) -> None
        Run the iterative method over a fixed number of time steps.
    """

    def __init__(
        self,
        domain: Mesh,
        space: FunctionSpace,
        phi: Function,
        tht: Function,
        diam2: float,
        smooth: bool,
    ) -> None:
        """
        Set up the iterative method to solve the level set equation.
        The bilinear and linear parts are pre-compiled.

        Parameters
        ----------
        domain : Mesh
            Problem domain.
        space : FunctionSpace
            Space of functions.
        phi : Function
            Level set function.
        tht : Function
            Velocity function.
        diam2 : float
            Square of the mesh diameter.
        smooth : bool
            Flag to add diffusion.
        """

        self.domain = domain
        self.dt = const(domain, 0.0)
        dx = Measure("dx", domain=domain)

        u = TrialFunction(space)
        v = TestFunction(space)

        # # ----
        # sign_phi = conditional(
        #     lt(phi, -0.25), -0.25, conditional(gt(phi, 0.25), 0.25, phi)
        # )
        # # ----

        # Petrov-Galerkin test function
        tau = 2.0 * sqrt(1.0 / self.dt**2 + dot(tht, tht) / diam2)
        new_v = v + dot(grad(v), tht) / tau

        # Crank-Nicolson weak formulation
        a = (u + (self.dt / 2.0) * dot(tht, grad(u))) * new_v * dx
        L = (phi - (self.dt / 2.0) * dot(tht, grad(phi))) * new_v * dx

        # Add diffusion
        if smooth:
            a += (self.dt / 2.0) * diam2 * dot(grad(u), grad(v)) * dx
            L -= (self.dt / 2.0) * diam2 * dot(grad(phi), grad(v)) * dx

        # Pre-compilation
        self.a = form(a)
        self.L = form(L)

    def run(self, phi: Function, steps: int, tend: float) -> None:
        """
        Run the iterative method over a fixed number of time steps.
        Since the bilinear part remains unchanged across iterations,
        it is compiled only once. Therefore, the linear part must be
        correctly updated.

        Parameters
        ----------
        phi : Function
            Level set function.
        steps : int
            Number of time steps.
        tend : float
            Final time.
        """

        self.dt.value = tend / steps

        solver = build_solver(self.domain, self.a, [])
        b = create_vector(self.L)

        for _ in range(steps):
            with b.localForm() as loc_b:
                loc_b.set(0)
            assemble_vector(b, self.L)
            b.ghostUpdate(
                addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE
            )
            solver.solve(b, phi.x.petsc_vec)
            phi.x.scatter_forward()

__init__(domain, space, phi, tht, diam2, smooth)

Set up the iterative method to solve the level set equation. The bilinear and linear parts are pre-compiled.

Parameters:

Name	Type	Description	Default
domain	Mesh	Problem domain.	required
space	FunctionSpace	Space of functions.	required
phi	Function	Level set function.	required
tht	Function	Velocity function.	required
diam2	float	Square of the mesh diameter.	required
smooth	bool	Flag to add diffusion.	required

Source code in code/formopt.py

def __init__(
    self,
    domain: Mesh,
    space: FunctionSpace,
    phi: Function,
    tht: Function,
    diam2: float,
    smooth: bool,
) -> None:
    """
    Set up the iterative method to solve the level set equation.
    The bilinear and linear parts are pre-compiled.

    Parameters
    ----------
    domain : Mesh
        Problem domain.
    space : FunctionSpace
        Space of functions.
    phi : Function
        Level set function.
    tht : Function
        Velocity function.
    diam2 : float
        Square of the mesh diameter.
    smooth : bool
        Flag to add diffusion.
    """

    self.domain = domain
    self.dt = const(domain, 0.0)
    dx = Measure("dx", domain=domain)

    u = TrialFunction(space)
    v = TestFunction(space)

    # # ----
    # sign_phi = conditional(
    #     lt(phi, -0.25), -0.25, conditional(gt(phi, 0.25), 0.25, phi)
    # )
    # # ----

    # Petrov-Galerkin test function
    tau = 2.0 * sqrt(1.0 / self.dt**2 + dot(tht, tht) / diam2)
    new_v = v + dot(grad(v), tht) / tau

    # Crank-Nicolson weak formulation
    a = (u + (self.dt / 2.0) * dot(tht, grad(u))) * new_v * dx
    L = (phi - (self.dt / 2.0) * dot(tht, grad(phi))) * new_v * dx

    # Add diffusion
    if smooth:
        a += (self.dt / 2.0) * diam2 * dot(grad(u), grad(v)) * dx
        L -= (self.dt / 2.0) * diam2 * dot(grad(phi), grad(v)) * dx

    # Pre-compilation
    self.a = form(a)
    self.L = form(L)

run(phi, steps, tend)

Run the iterative method over a fixed number of time steps. Since the bilinear part remains unchanged across iterations, it is compiled only once. Therefore, the linear part must be correctly updated.

Parameters:

Name	Type	Description	Default
phi	Function	Level set function.	required
steps	int	Number of time steps.	required
tend	float	Final time.	required

Source code in code/formopt.py

def run(self, phi: Function, steps: int, tend: float) -> None:
    """
    Run the iterative method over a fixed number of time steps.
    Since the bilinear part remains unchanged across iterations,
    it is compiled only once. Therefore, the linear part must be
    correctly updated.

    Parameters
    ----------
    phi : Function
        Level set function.
    steps : int
        Number of time steps.
    tend : float
        Final time.
    """

    self.dt.value = tend / steps

    solver = build_solver(self.domain, self.a, [])
    b = create_vector(self.L)

    for _ in range(steps):
        with b.localForm() as loc_b:
            loc_b.set(0)
        assemble_vector(b, self.L)
        b.ghostUpdate(
            addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE
        )
        solver.solve(b, phi.x.petsc_vec)
        phi.x.scatter_forward()