A Solver for Problems of the p-Laplacian PDE Operator.

This package provides a solver for problems of the scalar or vector-valued p-Laplacian with finite p including source terms and mixed Dirichlet-Neumann boundary conditions. The solver works iteratively based on a piece-wise linear finite element discretization and interior-point methods. Theory can be found in the publications [1] and [3] or summarized in the documentation.

The implementation is based on the finite element library MinFEM. Thus, it is able to import meshes in GMSH v1, v2 and v4 format and outputs VTK format for Paraview as well as statistics in plain txt.

Start by adding the PLaplace package to our julia installation and test it. Note that the test might take some time depending on your machine since it solves a full validation problem provided by the method of manufactured solutions [3]. Thus, open the julia REPL, hit the ] key and type

add PLaplace
test PLaplace

An Example Problem

Lets go through a code for the p-Laplace equation on a unit square once with inhomogeneous Dirichlet boundary conditions and once with mixed homogenous Dirichlet and inhomogeneous Neumann boundary conditions.

First we have to load the package PLaplace. Here we also load the package MinFEM to gain access to it's mesh generation functions. With these we then generate a uniform, triangular 30x30 mesh for the unit square.

using PLaplace, MinFEM

mesh = unit_square(30)

As an alternative: Download the package from github to obtain the examples and meshes and navigate, within the julia console, to the examples folder. Then import one of the mesh files generated with GMSH. Here, the additional include of MinFEM is not necessary.

using PLaplace

mesh = import_mesh("meshes/square.msh")

As preparation, it is often useful to create or specify a directory for all the output files.

output_path::String = "results/readme-example/"
mkpath(output_path)

In particular, we can also prepare a file where we can later write statistics. In this command, the guarded flag means that the file is not overwritten if it already exist. This way, we can keep this command in the script even if we already used it before and just want to append new results, e.g., for different parameters to the file.

statistics_file = output_path * "statistics.txt"
write_statistics_header(statistics_file, guarded = true)

In the next step we specify the boundary conditions. As the default procedure we simply specify julia functions as shown here for a non-homogenous Dirichlet boundary condition.

g(x) = x[1]^2

The same has to be done for mixed boundary conditions denoting one function for the Dirichlet part and one for the Neumann part. Note that also homogenous conditions need to be specified.

g(x) = 0
h(x) = x[2]^2 - x[1]^2

Alternatively, the package also support providing all boundary conditions and the source term as discrete values. Either as FEM coefficients on the nodes or values on quadrature points of the boundary elements. Note that this procedure is only recommended for experienced users or if an analytical description of the function is not available.

g = evaluate_mesh_function(mesh, x -> x[1]^2)

Next, we need to specify the boundary sets for the respective conditions. For pure Dirichlet problems it is easiest to select all physical boundaries of the mesh.

dirichlet_boundary = select_boundaries(mesh)

For mixed problems it is required to select the boundaries manually. Note that the mesh in this example is designed to have four physical boundaries identified by the indices 1001-1004.

dirichlet_boundary = select_boundaries(mesh, 1001, 1004)
neumann_boundary = select_boundaries(mesh, 1002, 1003)

Finally, we solve the problem and then write the solution in a file for visualization with Paraview. The object data of the custom type PLaplaceData contains further relevant data generated during the iteration. The most important parameters, such as iteration counts and obtained accuracy, can now be stored to the previously created statistics file.

p = 3.0
data = solve_plaplace(p, mesh, g, dirichlet_boundary)

write_statistics(statistics_file, data)

write_result_to_vtk(output_path * "result_p=$p", data)

While the solution routine always needs the parameter p, the mesh and the Dirichlet boundary information, it can be provided with a number of optional keyword arguments. The most prominent ones for Neumann boundary handling are shown here. For a full list including solver parameters see the documentation.

data = solve_plaplace(
    p,
    mesh,
    g,
    dirichlet_boundary,
    h = h,
    neumann_boundary = neumann_boundary
)

Both of the presented example and many more are available in the examples directory for experimenting.

References

[1] S. Loisel. “Efficient algorithms for solving the p-Laplacian in polynomial time”. Numerische Mathematik 146(2). 2020.
[2] K. Salari and P. Knupp. “Code Verification by the Method of Manufactured Solutions”. Sandia Report. 2000.
[3] H. Wyschka and M. Siebenborn.“Towards computing high-order p-harmonic descent directions and their limits in shape optimization”. In: "Modeling, Simulation and Optimization of Fluid Dynamic Applications". 2023.