Issues Solving Pure Neumann Problems in dolfinx

dolfinx
1

Is there any currently recommended way to solve pure neumann problems in Dolfinx?

This seems to be a major limitation for me using dolfinx. There are many cases were we would need to enforce that some average condition for the solution (contraining rigid body modes in mechanics, solving the pressure-poisson equation with incompressible fluid flows, etc), but currently I don’t know of a recommended way to solve these types of problems in dolfinx. The normal way (involving the use of LaGrange multipliers) works fine in legacy dolfin:

https://fenicsproject.org/olddocs/dolfin/latest/python/demos/neumann-poisson/demo_neumann-poisson.py.html

but it is well know that this is not possible in dolfinx because of the Real Type element does not exist:

Does anyone have a work around or a best practice until real type elements are fully implemented?

2

You can remove the constant nullspace with PETSc. See for instance

1 Like
3

Thanks @dokken.

I was able to replicate the Singular Poisson problem in dolfinx using your help. The dolfin example is here.

It might be nice to include this example in the dolfinx tutorial. Either way, I’ll just post the code I was able to get working here (apologies if this already exists elsewhere, it just took me a bit to put the pieces together):

from mpi4py import MPI
from dolfinx import mesh,plot
from dolfinx.fem import FunctionSpace,Function,form
from dolfinx.fem.petsc import create_vector,assemble_matrix,assemble_vector
from ufl import inner,TrialFunction,TestFunction,ds,dx,grad,exp,sin,SpatialCoordinate
from petsc4py import PETSc
import pyvista

# Create mesh and define function space
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20,20, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)

x = SpatialCoordinate(domain)
f = 10*exp(-(((x[0] - 0.5)**2) + (x[1] - 0.5)**2) / 0.02)
g = -sin(5*x[0])

a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

# Assemble system
A = assemble_matrix(form(a))
A.assemble()
b=create_vector(form(L))
with b.localForm() as b_loc:
            b_loc.set(0)
assemble_vector(b,form(L))

# Solution Function
uh = Function(V)

# Create Krylov solver
solver = PETSc.KSP().create(A.getComm())
solver.setOperators(A)

# Create vector that spans the null space
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
A.setNullSpace(nullspace)

# orthogonalize b with respect to the nullspace ensures that 
# b does not contain any component in the nullspace
nullspace.remove(b)

# Finally we are able to solve our linear system ::
solver.solve(b,uh.vector)

#plotting stuff
tdim = domain.topology.dim
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True,opacity=0.25)
u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh.x.array.real
u_grid.set_active_scalars("u")
plotter.add_mesh(u_grid.warp_by_scalar("u",factor=.5), show_edges=True)
plotter.view_vector((-0.25,-1,0.5))
if not pyvista.OFF_SCREEN:
    plotter.show()
1 Like
4

If you want to add it to the tutorial, feel free to make a Pull request at:
https://github.com/jorgensd/dolfinx-tutorial/pulls
and I will review it.