How can I use Block Forms, Multiphenicsx, and Lagrange Multipliers to Solve Poisson's Equation with Pure Neumann Boundary Conditions?

dolfinx
1

Hi!

I am using dolfinx version 0.6.0. It was installed with conda-forge. I am also using multiphenicsx. My aim is to solve a complicated set of linear elasticity equations with Robin boundary conditions using multiphenicsx and block forms. Also, I eventually want to impose additional constrains with Lagrange multipliers.

As a starting point, I am trying to solve Poisson’s equation with Pure Neumann boundary conditions. To do so, I am trying to translate this code using multiphenicsx and block forms. Below is a minimal working code example that shows how I am trying to do it. The output that this code produces shows that I am doing something terribly wrong.

How should I write Poisson’s equation with Pure Neumann boundary conditions and Lagrange multipliers in block form? Why am I getting such wrong output? What should I fix in the code below to solve the problem correctly?

 import numpy as np
 from mpi4py import MPI
 import dolfinx.fem
 import dolfinx.fem.petsc
 from dolfinx import mesh, plot
 from petsc4py import PETSc
 import ufl
 import pyvista
 import multiphenicsx.fem
 import multiphenicsx.fem.petsc

 def solve_poisson_demo():

    # MESH

    # mpi vars
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()
    size = comm.Get_size()

    # rectangular mesh and facets
    domain = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                                           points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                                           cell_type=mesh.CellType.triangle,)

    facets = mesh.locate_entities_boundary(domain, dim=1,
                                           marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                           np.isclose(x[0], 2.0)))

    # FUNCTION SPACES / BLOCK FUNCTION SPACE

    V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2))
    R = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2))
    (u, c) = (ufl.TrialFunction(V), ufl.TrialFunction(R))
    (v, d) = (ufl.TestFunction(V), ufl.TestFunction(R))


    # IMPORTANT EXPRESSIONS

    f = dolfinx.fem.Function(V)
    f.interpolate(lambda x:  10 * np.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02))
    g = dolfinx.fem.Function(V)
    g.interpolate(lambda x: np.sin(5 * x[0]))
    _zero = dolfinx.fem.Function(V)
    _zero.interpolate(lambda x: np.zeros(len(x[0])))
    dx = ufl.Measure("dx",domain=domain)
    ds = ufl.Measure("ds",domain=domain)


    # LINEAR AND BILINEAR FORMS

    aa = [[ufl.inner(ufl.grad(u),ufl.grad(v))*dx,v*c*dx],
             [u*d*dx,None]]
    ff = [ufl.inner(f,v)*dx + ufl.inner(g,v)*ds,_zero*d*dx]
    aa_cpp = dolfinx.fem.form(aa)
    ff_cpp = dolfinx.fem.form(ff)


    # ASSEMBLE BLOCK EQUATIONS

    AA = dolfinx.fem.petsc.assemble_matrix_block(aa_cpp, bcs=[])
    AA.assemble()
    FF = dolfinx.fem.petsc.assemble_vector_block(ff_cpp, aa_cpp, bcs=[])


    # SOLVE

    uu = dolfinx.fem.petsc.create_vector_block(ff_cpp)

    # solver
    solver = PETSc.KSP().create(domain.comm)
    solver.setOperators(AA)
    solver.setType("preonly")
    solver.getPC().setType("lu")
    solver.getPC().setFactorSolverType("superlu_dist")

    # solve
    solver.solve(FF,uu)
    uu.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    solver.destroy()


    # SPLIT

    # Split the block solution in components
    (uh,c) = (dolfinx.fem.Function(V), dolfinx.fem.Function(R))
        with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(uu, [V.dofmap, R.dofmap]) as 
        uu_wrapper:
            for u_wrapper_local, component in zip(uu_wrapper, (uh,c)):
                with component.vector.localForm() as component_local:
                    component_local[:] = u_wrapper_local


     # POST-PROCESS

     if comm.Get_size() == 1:

        # plot displacement
        topology, cell_types, geometry = plot.create_vtk_mesh(V)
        u = uh.x.array
        grid = pyvista.UnstructuredGrid(topology,cell_types,geometry)
        grid.point_data["u (cm)"] = u
        warped = grid.warp_by_scalar("u (cm)", factor=1.0)
        plotter = pyvista.Plotter()
        plotter.add_mesh(warped, show_scalar_bar=True, show_edges=False, scalars="u (cm)")
        plotter.show(screenshot='poisson_demo_visualize.png')
        plotter.clear()

     return None

 if __name__=='__main__':
     solve_poisson_demo()
2

Can you discuss where for which part of the code do you plan to use multiphenicsx? (I guess the Lagrange multipliers that are not there yet, but I want to be sure)
As the author of that, I cannot see it being actually used right now [1].

On pure Neumann boundary conditions, you may want to have a look at Singular poisson equation in DOLFINx

Together with @hherlyng in the next few weeks we would like to add a tutorial on pure Neumann BCs in multiphenicsx. Please get in touch with me by email to remind me that you are interested, and to post a link here when we are done.

[1] I do understand that you are indeed using multiphenicsx.fem.petsc.BlockVecSubVectorWrapper. I like that approach, otherwise I wouldn’t have prepared it ;), but just keep in mind that until you actually have restrictions you can do a similar operation using plain dolfinx, see e.g.

As soon as you will have multiphenicsx restrictions, then multiphenicsx.fem.petsc.BlockVecSubVectorWrapper will be the only way to go.

1 Like
3

Dear Francesco,

Thanks for your quick response.

In the non-MWE version of this code I use a disk. Eventually, I will need to implement an additional boundary condition on the top of the disk for the linear elasticity problem. That is where multiphenicsx and restrictions are handy.

In the MWE code I provided, the functions c and d are supposed to be Lagrange multipliers (just as in the code provided in this tutorial. I am trying to assemble the equations provided in that tutorial into block forms. Do you know how I can assemble those equations into block forms correctly?

I am very interested in that tutorial. I will get in touch with you by email.

4

The final implementation, that will resemble the one in the post I linked above, will be different from the old tutorial.

In particular

R = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 2))

is not an equivalent of the former

R = FunctionSpace(mesh, "R", 0)

(i.e., a FE space with a single global degree of freedom). The former "R" space has not been reimplemented yet in dolfinx, and that is one of the reasons why the new implementation will need to be different.

1 Like
5

The new tutorial in multiphenicsx about this topic is now available at

(joint work with @hherlyng )

Please keep in mind that the link above points to the notebook on the main branch. I do not backport new tutorials to branches related to older dolfinx versions. However, as of now (December 2023), I think that only a couple of minor changes are required for you to run that tutorial with dolfinx 0.6.0. Just download that tutorial and:

  1. replace every dolfinx.fem.functionspace with dolfinx.fem.FunctionSpace,
  2. drop the import of viskex, and every cell containing it. (I started working on viskex after the dolfinx 0.6.0 release).

(Future readers should be aware that items 1 and 2 may not necessarily apply in future, as the tutorial on the main branch may keep changing as dolfinx/multiphenicsx development continues)

3 Likes