Implement DirichiletBC in dolfinx (ksp solver)

Hello FEniCS family,

I am getting difficulty in implementing bc to ksp solver. Please find MWE:

    from mpi4py import MPI
    from dolfinx import mesh, fem, cpp
    mesh = mesh.create_unit_square(MPI.COMM_WORLD, 10,10, mesh.CellType.triangle)
    from dolfinx.fem import VectorFunctionSpace, form, petsc, Function
    import dolfinx 
    import ufl
    import numpy as np
    import petsc4py.PETSc
    V = VectorFunctionSpace(domain, ("Lagrange", 1),dim=3)
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    F2 = dot(u,v)*dx
    A=  petsc.assemble_matrix(form(ufl.lhs(F2)))
    F = petsc.assemble_vector(form(ufl.rhs(F2)))
    F.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
    # ksp solve
    ksp = petsc4py.PETSc.KSP()
    ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)  # detect null pivots
    ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)  # do not compute null space again
    ksp.solve(F, w.vector)
        addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
   # Apply Dirichilet BC
    uD = Function(V)
    uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
    tdim = domain.topology.dim
    fdim = tdim - 1
    mesh.topology.create_connectivity(fdim, tdim)
    boundary_facets = mesh.exterior_facet_indices(mesh.topology)
    boundary_dofs = locate_dofs_topological(V, fdim, boundary_facets)
    bc = dolfinx.fem.dirichletbc(uD, boundary_dofs)

My query is how can I implement bc to ksp solve? I choose ksp solve to implement the nullspace through operator A. The above bc is similar to dolfinx example.

Any help is greatly appreciated.

Do you want to apply the DirichletBC to the problem after solving without BCS?
It is common to apply Dirichlet conditions prior to solving, with either of the approaches described at: Application of Dirichlet boundary conditions — FEniCS Tutorial @ Sorbonne

Thanks @dokken, I want to apply bc before solving the variational form.

I am little confused about process of solving var form. Can the above MWE

be written in terms of

# In most situations, you would just pass this object to the linear problem and it would be handled for you
petsc_options = {"ksp_type": "preonly",
                 "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
problem = dolfinx.fem.petsc.LinearProblem(
    a, L, bcs=bcs, petsc_options=petsc_options)
u = problem.solve()         

including setOperator(A) (like in here). So that in solve (I can input assembled matrix A and assembled vector F instead of form a and form f.
Any help is greatly appreciated.

Yes, you can set these values through the options, or after initialization of problem access the ksp-object through problem.solver, see: dolfinx/python/dolfinx/fem/ at 387ff1279c3470fb4956e54f47e00ce81d2f9c1d · FEniCS/dolfinx · GitHub

1 Like

@dokken Thanks for your response.

Can we implement bcs to assembled matrix using apply_lifting. (mean, we assemble the matrix first and then apply bcs, like in assembling of PETSc vector through apply_lifting feature)

For case of this nullspace example from author @francesco-ballarin, I wanted to apply the nullspace constraints before setting up the bc.

However, A=assemble_matrix(a,bcs) sets bcs before applying constraints, doing so, I got the nullspace assertion error assert nullspace.test(A). Can you guide a way such that, both nullspace and bc application works correctly.

Q: Does nullspace of a matrix change when we impose dirichilet to it?

Of course it does. Given a solution u, the whole point of a null space is to find an increment n such that u + n is still a solution. If you Dirichlet BCs on some DOFs of u, then n on those DOFs must be zero, other the sum u + n will not be a solution anymore, at least on those DOFs.

1 Like

@dokken , I have a followup question.

  1. Does setting petsc_options in petsc.LinearProblem as "petsc_options = {"getFactorMatrix().setMumpsIcntl":(24,1), "getFactorMatrix().setMumpsIcntl":(25,0),"ksp_Operators": A} is correct syntax for

I guess so. I usually don’t use that way of enforcing petsc options. Please use ksp.view(), and ksp.getPC().view() to inspect the options set.