# Singular poisson equation in DOLFINx

Hello, I want to ask a question regarding the solution to singular systerm.

It is introduced here Issues Solving Pure Neumann Problems in dolfinx that we can solve a singular system Ax=b by removing the null space of A from b if the matrix A is singular.

And I was able to solver the Pure Neuman Poisson discussed there

But I find that sometimes it doesn’t give the correct solution

If we simply assume
-\Delta u = f at [0,1] where u = sin(2\pi x), then we can get f = 4\pi^2sin(2\pi x)
the weak form should be
a = inner(grad(u),grad(v))*dx
L = f * v * dx + g * v * ds,
where f = 4\pi^2sin(2\pi x), g = 2\pi cos(2\pi x). (Maybe there is something wrong here?)

After I solve it, the solution is very different from the correct one, does anyone know what is the reason?
I also want to ask, is this method can be applied together with DOLFINX_MPC, say if I want to solve a Poisson equation with pure periodic BC, is it possible in DOLFINX? In the old fenics, I usually solved it by using Lagrange multiplyer, but it seems like the Real element has not been implement in dolfinx yet.

Thanks!

As you have not supplied your code here. it is hard to give any guidance.
I would have a look at: Add singular Poisson · Issue #139 · jorgensd/dolfinx-tutorial · GitHub

I think it would work with DOLFINx_MPC.

Hello, thanks for replying! I have read the page that you suggested, and I was able to get some solution of singular system with PBC, but the soltion seems to be very wird,

Is this a bug? or maytbe I have implement something wrong.
I attached my codes here:

int main(int argc, char *argv[])
{
dolfinx::init_logging(argc, argv);
PetscInitialize(&argc, &argv, nullptr, nullptr);
{
// Define basic parameters
auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
auto mesh = std::make_shared<mesh::Mesh>(mesh::create_interval(MPI_COMM_WORLD, N, {xmin, xmax}, part));

// Define function space
auto V = std::make_shared<fem::FunctionSpace>(fem::create_functionspace(functionspace_form_poisson_a, "u", mesh));

// Prepare and set Constants for the bilinear form
auto f = std::make_shared<fem::Function<T>>(V);
auto u = std::make_shared<fem::Function<T>>(V);

// Define variational forms
std::map<std::string, std::shared_ptr<const fem::Function<T>>> empty_map;
auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(*form_poisson_a, {V, V}, empty_map, {}, {}));
auto L = std::make_shared<fem::Form<T>>(fem::create_form<T>(*form_poisson_L, {V}, {{"f", f}}, {}, {}));

f->interpolate(source);

io::VTKFile file1(MPI_COMM_WORLD, "results/f.pvd", "w");
file1.write<T>({*f}, 0.0);

auto pbc = dolfinx_mpc::create_periodic_condition_geometrical(V, periodic_condition_x, periodic_relation_x, {}, 1, false);
auto mpc = std::make_shared<dolfinx_mpc::MultiPointConstraint<T>>(dolfinx_mpc::MultiPointConstraint<double>(V, {pbc.slaves}, {pbc.masters}, {pbc.coeffs}, {pbc.owners}, {pbc.offsets}));

LinearProblem Poisson_eq(a, L, mpc, {});
Poisson_eq.assemble_Mat();
Poisson_eq.assemble_Vec();
Poisson_eq.solve(*u);

io::VTKFile file2(MPI_COMM_WORLD, "results/u.pvd", "w");
file2.write<T>({*u}, 0.0);
}

PetscFinalize();

return 0;
}


Here are my source codes:

using namespace dolfinx;
using T = PetscScalar;

//-----------------------------------------------------------------------------
LinearProblem::LinearProblem(
std::shared_ptr<fem::Form<T>> a,
std::shared_ptr<fem::Form<T>> L,
std::shared_ptr<dolfinx_mpc::MultiPointConstraint<T>> mpc,
std::vector<std::shared_ptr<const fem::DirichletBC<T>>> bcs) : _a(a), _L(L), _mpc(mpc), _bcs(bcs),
_b(L->function_spaces()[0]->dofmap()->index_map,
L->function_spaces()[0]->dofmap()->index_map_bs())
{
auto pattern = dolfinx_mpc::create_sparsity_pattern(*a,mpc,mpc);
pattern.assemble();
_matA = std::make_shared<la::petsc::Matrix>(MPI_COMM_WORLD,pattern);
lu = std::make_shared<la::petsc::KrylovSolver>(MPI_COMM_WORLD);
la::petsc::options::set("ksp_type", "gmres");
la::petsc::options::set("pc_type", "hypre");
lu->set_from_options();
}

//-----------------------------------------------------------------------------
LinearProblem::~LinearProblem()
{
}

//-----------------------------------------------------------------------------
void LinearProblem::assemble_Mat()
{

MatZeroEntries(_matA->mat());

dolfinx_mpc::assemble_matrix(la::petsc::Matrix::set_block_fn(_matA->mat(), ADD_VALUES), la::petsc::Matrix::set_fn(_matA->mat(), ADD_VALUES), *_a, _mpc, _mpc, {}, 1.0);
MatAssemblyBegin(_matA->mat(), MAT_FLUSH_ASSEMBLY);
MatAssemblyEnd(_matA->mat(), MAT_FLUSH_ASSEMBLY);
fem::set_diagonal<T>(la::petsc::Matrix::set_fn(_matA->mat(), INSERT_VALUES), *_L->function_spaces()[0], {_bcs});
MatAssemblyBegin(_matA->mat(), MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(_matA->mat(), MAT_FINAL_ASSEMBLY);
lu->set_operator(_matA->mat());
}

//-----------------------------------------------------------------------------
void LinearProblem::assemble_Vec()
{
_b.set(0.0);

dolfinx_mpc::assemble_vector(_b.mutable_array(), *_L, _mpc);
dolfinx_mpc::apply_lifting(_b.mutable_array(), {_a}, {{_bcs}}, {}, 1.0, _mpc);
_b.scatter_rev(std::plus<T>());
fem::set_bc(_b.mutable_array(), {_bcs});
}

//-----------------------------------------------------------------------------
void LinearProblem::solve(fem::Function<T> &u)
{
la::petsc::Vector _bv(la::petsc::create_vector_wrap(_b), false);
la::petsc::Vector _u(la::petsc::create_vector_wrap(*u.x()), false);

PetscBool has_cnst = PETSC_TRUE;
PetscInt num_vectors = 0;
Vec * vec=NULL;
MatNullSpace nullspace;
PetscErrorCode ierr = MatNullSpaceCreate(MPI_COMM_WORLD, has_cnst, num_vectors, vec, &nullspace);
MatSetNullSpace(_matA->mat(),nullspace);
MatNullSpaceRemove(nullspace, _bv.vec());
if (ierr == 0)
MatNullSpaceDestroy(&nullspace);

lu->solve(_u.vec(), _bv.vec());
_mpc->backsubstitution(u.x()->mutable_array());
}

//-----------------------------------------------------------------------------
std::pair<std::vector<T>, std::vector<std::size_t>> source(std::experimental::mdspan<const double, std::experimental::extents<std::size_t, 3, std::experimental::dynamic_extent>> x)
{
std::vector<T> u0;

for (std::size_t p = 0; p < x.extent(1); ++p)
{
double dx = x(0, p);
u0.push_back(4.0 * pi * pi * sin(2.0 * pi * dx));
}

return {u0, {u0.size()}};
}

//-----------------------------------------------------------------------------
std::vector<std::int8_t> periodic_condition_x(std::experimental::mdspan<const double, std::experimental::extents<std::size_t, 3, std::experimental::dynamic_extent>> x)
{
std::vector<std::int8_t> f;
for (std::size_t p = 0; p < x.extent(1); ++p)
{
f.push_back((std::abs(x(0, p)-xmax) < 1e-6));
}
return f;
}

//-----------------------------------------------------------------------------
std::vector<double> periodic_relation_x(std::span<const double> x)
{
std::vector<double> y(x.size(),0.0);

std::size_t len = x.size()/3;
for (std::size_t i=0; i < x.size()/3; i++)
{
y[i] = (xmin+xmax)-x[i];
y[i+len] = x[i+len];
y[i+2*len] = x[i+2*len];
}

return y;
}

import basix
from basix.ufl_wrapper import create_element
from ufl import (Coefficient, TestFunction, TrialFunction, dx, ds, inner, grad, interval)

cell = "interval"
degree = 1
element = create_element("Lagrange",cell, degree)

u = TrialFunction(element)
v = TestFunction(element)

f = Coefficient(element)

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



Thanks very much!

Hello, I have tried some times again, and the results are:

1. If solve a singular Poisson equation without removing the nulll space but with PBC, I will get a solution whose left bound is equal to right bound, but the wird thing is that I should not be able to solve it since it is singular system, and the range of sin function is not that correct ([-1,1])

2. If solve a singular Poisson equation with removing the nulll space and with PBC, the solution
the right bound will be equal to the 2nd left bound, as you can see there is a little “tail” in the solution, but the range is exactly [-1,1]

Maybe there is something wrong withe KSP options, do you have any ideas or suggestions?

_b is created with the wrong dofmap. The dofmap used by MPC can be found in function space of the MPC: https://github.com/jorgensd/dolfinx_mpc/blob/v0.7.0/cpp/MultiPointConstraint.h#L192-L195

BTW.: Well done making a C++ code for DOLFINx_mpc. If you want to add a C++ demo to DOLFINx_mpc, a PR is welcome!

Hi, thanks! Our research group is moving to the new DOLFINX, and we will work in the CPP interface. Hopefully we can ontribute some codes and demos, after we make our code clear and convincible.

I’m interested in the principle of the two methods used to solve the singular equation, or to solve the system with pure Neumann boundary conditions.

The old method available in legacy dolfinx applied a Lagrange multiplier to resolve the ill-conditioned system by imposing the condition

\int_\Omega u dx= 0

(edit: mutlphenicsx can also handle Lagrange multipliers)

The new method available in dolfinx via PETSc mentioned in this thread removes the null space of the matrix.

It’s not entirely obvious to me that the two methods must be exactly the same. Are the solutions obtained in both cases necessarily identical? I’d be interested if anyone knows any references that discuss this.

Hi, in the old dolfin, you can also solve the singular Poisson by removing the null space, see here. But for the new DOLFINx, the real element is not avaliable so far, so you can not use the Lagrange multiplyer. As for the difference, I guess maybe it depends on what KSP solver, the method which removes the null space is sensitive to the KSP solver in my case, I’m actually studying this but haven’t got the answer, maybe @dokken can make a comment on it.

Ah, true, I had missed the singular Poisson example in legacy dolfin.

Sure, the specific numerical solutions of the two methods will differ within numerical tolerance. But are they in principle the same solution?

The main difference from a numerical point of view is that a real space introduces a dense column/row in your system, meaning that the scalability/performance of a solver will be reduced for large problem sizes.

PETSc has written some notes on what nullspace removal does mathematically: MatSetNullSpace — PETSc 3.20.1 documentation

Better performance is certainly a good reason to prefer the null space method. Still, at the level of the mathematics rather than the numerics, is the solution the same? If we were to execute both methods on an ideal Türing engine, would they return the same perfect solution?

Our group has recently been working with a Neumann problem and we have been asking exactly the same question, i.e. whether the two approaches of

1. constraining the solution to have a mean equal to zero and
2. setting the nullspace as in e.g. Add singular Poisson

actually represents the same problem and thus have the same solutions. For a time-dependent problem, I myself cannot see these as being identical.

Would love to get some input on this if anyone has some!