Blocked solver for mixed elements

General
1

Hi!
I’m trying to use block preconditioners in mixed finite element, but I guess my my mapping is wrong. I’m get the following error:

petsc4py.PETSc.Error: error code 60
[15] SNESSolve() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/snes/interface/snes.c:4692
[15] SNESSolve_NEWTONLS() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/snes/impls/ls/ls.c:210
[15] KSPSolve() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/ksp/ksp/interface/itfunc.c:1071
[15] KSPSolve_Private() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/ksp/ksp/interface/itfunc.c:825
[15] KSPSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/ksp/ksp/interface/itfunc.c:406
[15] PCSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/ksp/pc/interface/precon.c:994
[15] PCSetUp_FieldSplit() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/ksp/pc/impls/fieldsplit/fieldsplit.c:667
[15] MatCreateSubMatrix() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/mat/interface/matrix.c:8382
[15] MatCreateSubMatrix_MPIAIJ() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/mat/impls/aij/mpi/mpiaij.c:3427
[15] MatCreateSubMatrix_MPIAIJ_nonscalable() at /home/conda/feedstock_root/build_artifacts/petsc_1677582998761/work/src/mat/impls/aij/mpi/mpiaij.c:3809
[15] Nonconforming object sizes
[15] Local column sizes 132618 do not add up to total number of columns 8484

Below follows a mwe of my current code:

V = VectorElement("CG", domain.ufl_cell(), 1)
Q = FiniteElement("CG", domain.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(domain, solution_space)
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()

w = Function(W)
(u, p) = w.split()
(du, dp) = TrialFunctions(W)
(v, q) = TestFunctions(W)

F_momentum = (
    dot(dot(u,nabla_grad(u)),v)*dx
    + (1/Re)*inner(grad(u),grad(v))*dx
    - inner(p,div(v))*dx
    + inner(tau_SUPG*grad(v)*u,dot(u,grad(u)) - (1/Re)*div(grad(u)) + grad(p))*dx    
)
F_continuity = (
    (q*div(u))*dx
    + inner(tau_PSPG*grad(q),  dot(u,grad(u)) - (1/Re)*div(grad(u)) + grad(p))*dx
)
F = form([F_momentum, F_continuity])
jacobian = form(
    [
        [derivative(F_momentum, u, du)   , derivative(F_momentum, p, dp)],
        [derivative(F_continuity, u, du) , derivative(F_continuity, p, dp)],
    ]
)
P = form([[jacobian[0][0], None], [None, inner(dp,q)*dx]])

# Assemble Jacobian and mass matrices and residual vector
jacobian_matrix = create_matrix_block(jacobian)
P_matrix        = create_matrix_block(P)
residual_vector = create_vector_block(F)

# Create nonlinear solver
nonlinear_solver = PETSc.SNES().create(comm)
nonlinear_solver.getKSP().setType("fgmres")
nonlinear_solver.getKSP().getPC().setType("fieldsplit")

V_map = W0.dofmap.index_map
Q_map = W1.dofmap.index_map
offset_u = V_map.local_range[0]*W0.dofmap.index_map_bs + Q_map.local_range[0]
offset_p = offset_u + V_map.size_local*W0.dofmap.index_map_bs
is_u = PETSc.IS().createStride(V_map.size_local*W0.dofmap.index_map_bs, offset_u, 1, comm=PETSc.COMM_SELF)
is_p = PETSc.IS().createStride(Q_map.size_local, offset_p, 1, comm=PETSc.COMM_SELF)
nonlinear_solver.getKSP().getPC().setFieldSplitIS(("u", is_u),("p", is_p))

ksp_u, ksp_p = nonlinear_solver.getKSP().getPC().getFieldSplitSubKSP()
# Velocity block
ksp_u.setType("preonly")
ksp_u.getPC().setType("hypre")
# Pressure block
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

# Solve nonlinear problem
problem = Nonlinear_PDE_SNES_Problem(F, jacobian, [u, p], boundary_conditions_navier_stokes, P)
nonlinear_solver.setFunction(problem.F_block, residual_vector)
nonlinear_solver.setJacobian(problem.J_block, jacobian_matrix, P_matrix)
# Create solution vector
x = create_vector_block(F)
u = Function(W0)
p = Function(W1)
with u.vector.localForm() as _u, p.vector.localForm() as _p:
    scatter_local_vectors(
        x,
        [_u.array_r, _p.array_r],
        [
            (u.function_space.dofmap.index_map, u.function_space.dofmap.index_map_bs),
            (p.function_space.dofmap.index_map, p.function_space.dofmap.index_map_bs),
        ],
    )
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
nonlinear_solver.solve(None, x)

Can anyone figure out what’s going on? Thank you.

2

Without having the full context, ie imports and a mesh, it is quite a high threshold to help you debug this.
Please try to set up a simpler form, say a projection in each sub space, and see if you can get that to work or throw an error similar to the one you have atm.

3

Hi Dokken,

thank you for your quick response. After many hours unfortunately I wasn’t able to figure out a solution exactly the way you suggested. However, I created the code below, which is simpler than what I posted previously. This code should be fully reproducible and generates the same issue I first reported. Can you help me understanding what’s going on?

from mpi4py import MPI
from dolfinx.cpp.la.petsc import scatter_local_vectors
from dolfinx.fem import (
    Function, FunctionSpace, dirichletbc, form, locate_dofs_topological
)
from dolfinx.fem.petsc import (
    assemble_matrix_block, assemble_vector_block, create_matrix_block,
    create_vector_block
)
from dolfinx.mesh import (
    CellType, locate_entities_boundary, GhostMode, create_rectangle
)
from ufl import (
    dx, grad, inner, derivative, TrialFunctions, dot, VectorElement,
    FiniteElement, MixedElement, TestFunctions, nabla_grad, div
)
from petsc4py.PETSc import ScalarType
import numpy as np
from dolfinx import mesh
import dolfinx
from petsc4py import PETSc
import warnings
warnings.filterwarnings("ignore")

mu = 0.0089                 
rho = 2.109502311
domain = create_rectangle(MPI.COMM_WORLD,
                     [np.array([0., 0.]), np.array([3., 1.])],
                     [45, 15],
                     CellType.triangle, GhostMode.none)

V = VectorElement("CG", domain.ufl_cell(), 2)
Q = FiniteElement("CG", domain.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(domain, solution_space)
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()

tolerance = 1e-6
bcs = []
def no_slip_boundary_front(x):
    return np.logical_and(x[0] < tolerance, (x[1]-0.5*1.0)**2 > 0.25**2)
def no_slip_boundary_bottom(x):
    return (x[1] < tolerance)
def no_slip_boundary_top(x):
    return (x[1] > (1.0 - tolerance))
def no_slip_boundary_back(x):
    return np.logical_and(x[0] > (3 - tolerance), (x[1]-0.5*1.0)**2 > 0.25**2)
def inlet_front_1(x):
    return np.logical_and(x[0] < tolerance, (x[1]-0.5*1.0)**2 <= 0.25**2)
def outlet_back_1(x):
    return np.logical_and(x[0] > (3 - tolerance), (x[1]-0.5*1.0)**2 <= 0.25**2)
def inlet_velocity_1(x):
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = (1/0.25)*(0.25**2 - (x[1]-0.5*1.0)**2)**0.5
    return values 

facets_no_slip_boundary_front  = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_front)
facets_no_slip_boundary_bottom = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_bottom)
facets_no_slip_boundary_top    = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_top)
facets_no_slip_boundary_back   = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_back)
facets_inlet_front_1           = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, inlet_front_1)
facets_outlet_back_1           = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, outlet_back_1)

B, _      = W.sub(0).collapse()
u_no_slip = Function(B)
u_no_slip.x.set(0)
bcu_no_slip_front  = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_front), W.sub(0))
bcu_no_slip_bottom = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_bottom), W.sub(0))
bcu_no_slip_top    = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_top), W.sub(0))
bcu_no_slip_back   = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_back), W.sub(0))

u_inlet_profile_1 = Function(B)
u_inlet_profile_1.interpolate(inlet_velocity_1)
bcu_inlet_1 = dirichletbc(u_inlet_profile_1, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front_1), W.sub(0))

G, _          = W.sub(1).collapse()
zero_pressure = Function(G)
zero_pressure.x.set(0)
bcu_outlet_1  = dirichletbc(zero_pressure, locate_dofs_topological((W.sub(1), G), domain.topology.dim - 1, facets_outlet_back_1), W.sub(1))

bcs = [bcu_no_slip_front, bcu_no_slip_bottom, bcu_no_slip_top, bcu_no_slip_back,\
       bcu_inlet_1, bcu_outlet_1]

class Nonlinear_PDE_SNES_Problem:
    def __init__(self, F, J, soln_vars, bcs, P=None):
        self.L = F
        self.a = J
        self.a_precon = P
        self.bcs = bcs
        self.soln_vars = soln_vars

    def F_block(self, snes, x, F):
        assert x.getType() != "nest"
        assert F.getType() != "nest"
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        with F.localForm() as f_local:
            f_local.set(0.0)

        offset = 0
        x_array = x.getArray(readonly=True)
        for var in self.soln_vars:
            size_local = var.vector.getLocalSize()
            var.vector.array[:] = x_array[offset : offset + size_local]
            var.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
            offset += size_local

        assemble_vector_block(F, self.L, self.a, bcs=self.bcs, x0=x, scale=-1.0)

    def J_block(self, snes, x, J, P):
        assert x.getType() != "nest" and J.getType() != "nest" and P.getType() != "nest"
        J.zeroEntries()
        assemble_matrix_block(J, self.a, bcs=self.bcs, diagonal=1.0)
        J.assemble()
        if self.a_precon is not None:
            P.zeroEntries()
            assemble_matrix_block(P, self.a_precon, bcs=self.bcs, diagonal=1.0)
            P.assemble()

w = Function(W)
(u, p)   = w.split()
(v, q)   = TestFunctions(W)
(du, dp) = TrialFunctions(W)

F_momentum = (
    rho*dot(dot(u,nabla_grad(u)),v)*dx
    + mu*inner(grad(u),grad(v))*dx
    - inner(p,div(v))*dx
)
F_continuity = (
    (q*div(u))*dx
)
F = form([F_momentum, F_continuity])
jacobian = form(
    [
        [derivative(F_momentum, u, du)  , derivative(F_momentum, p, dp)],
        [derivative(F_continuity, u, du), derivative(F_continuity, p, dp)],
    ]
)
P = form([[jacobian[0][0], None], [None, inner(dp,q)*dx]])

jacobian_matrix = create_matrix_block(jacobian)
P_matrix        = create_matrix_block(P)
residual_vector = create_vector_block(F)

options = PETSc.Options()
nonlinear_solver = PETSc.SNES().create(MPI.COMM_WORLD)
nonlinear_solver.getKSP().setType("fgmres")
nonlinear_solver.getKSP().getPC().setType("fieldsplit")
nonlinear_solver.getKSP().getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
nonlinear_solver.getKSP().getPC().setFieldSplitSchurFactType(PETSc.PC.SchurFactType.FULL)
nonlinear_solver.getKSP().getPC().setFieldSplitSchurPreType(PETSc.PC.SchurPreType.SELFP)

V_map = W0.dofmap.index_map
Q_map = W1.dofmap.index_map
offset_u = V_map.local_range[0]*W0.dofmap.index_map_bs + Q_map.local_range[0]
offset_p = offset_u + V_map.size_local*W0.dofmap.index_map_bs
is_u = PETSc.IS().createStride(V_map.size_local*W0.dofmap.index_map_bs, offset_u, 1, comm=PETSc.COMM_SELF)
is_p = PETSc.IS().createStride(Q_map.size_local, offset_p, 1, comm=PETSc.COMM_SELF)
nonlinear_solver.getKSP().getPC().setFieldSplitIS(("u", is_u),("p", is_p))

ksp_u, ksp_p = nonlinear_solver.getKSP().getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("hypre")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

nonlinear_solver.setFromOptions()
nonlinear_solver.getKSP().setFromOptions()

problem = Nonlinear_PDE_SNES_Problem(F, jacobian, [u, p], bcs, P)
nonlinear_solver.setFunction(problem.F_block, residual_vector)
nonlinear_solver.setJacobian(problem.J_block, jacobian_matrix, P_matrix)
x = create_vector_block(F)
u = Function(W0)
p = Function(W1)
with u.vector.localForm() as _u, p.vector.localForm() as _p:
    scatter_local_vectors(
        x,
        [_u.array_r, _p.array_r],
        [
            (u.function_space.dofmap.index_map, u.function_space.dofmap.index_map_bs),
            (p.function_space.dofmap.index_map, p.function_space.dofmap.index_map_bs),
        ],
    )
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
nonlinear_solver.solve(None, x)
1 Like
4

It seems like you are using an old version of DOLFINx.
What version are you running?
The example that you supplied, given some minor API changes does not fail on my system with DOLFINx v0.7.3, when executed in serial. It does however fail when you run in parallel.

Is this the case for you as well?

It seems to be your index sets that are wrong.
This is because you are using COMM_SELF in the createStride function, which should be

V_map = W0.dofmap.index_map
Q_map = W1.dofmap.index_map
offset_u = V_map.local_range[0]*W0.dofmap.index_map_bs + Q_map.local_range[0]
offset_p = offset_u + V_map.size_local*W0.dofmap.index_map_bs
is_u = PETSc.IS().createStride(V_map.size_local*W0.dofmap.index_map_bs,
                               offset_u, 1, comm=PETSc.COMM_WORLD)
is_p = PETSc.IS().createStride(Q_map.size_local, offset_p, 1, comm=PETSc.COMM_WORLD)
nonlinear_solver.getKSP().getPC().setFieldSplitIS(("u", is_u), ("p", is_p))

ksp_u, ksp_p = nonlinear_solver.getKSP().getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("hypre")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
5

Hi Dokken.

Thank you for your invaluable help. I was using the DOLFINx version 0.6.0, but I just updated to the same version you’re using and the code runs. I observe though that the norm of the solution is zero. I’ve checked boundary conditions and they seem correct to me. I wonder if this is due to a lack of initial guess for velocity and pressure. Do you have any idea why the solution is zero? Below follows my code updated according to version 0.7.3:

from mpi4py import MPI
from dolfinx.cpp.la.petsc import scatter_local_vectors
from dolfinx.fem import (
    Function, FunctionSpace, dirichletbc, form, locate_dofs_topological
)
from dolfinx.fem.petsc import (
    assemble_matrix_block, assemble_vector_block, create_matrix_block,
    create_vector_block
)
from dolfinx.mesh import (
    CellType, locate_entities_boundary, GhostMode, create_rectangle,
     meshtags
)
from ufl import (
    dx, grad, inner, derivative, TrialFunctions, dot, VectorElement,
    FiniteElement, MixedElement, TestFunctions, nabla_grad, div
)
from petsc4py.PETSc import ScalarType
import numpy as np
from dolfinx import mesh
import dolfinx
from petsc4py import PETSc
import warnings
warnings.filterwarnings("ignore")

mu = 0.0089                 
rho = 2.109502311
domain = create_rectangle(MPI.COMM_WORLD,
                     [np.array([0., 0.]), np.array([3., 1.])],
                     [45, 15],
                     CellType.quadrilateral)

V = VectorElement("CG", domain.ufl_cell(), 2)
Q = FiniteElement("CG", domain.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(domain, solution_space)
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()

tolerance = 1e-6
bcs = []
def no_slip_boundary_front(x):
    return np.logical_and(x[0] < tolerance, (x[1]-0.5*1.0)**2 > 0.25**2)
def no_slip_boundary_bottom(x):
    return (x[1] < tolerance)
def no_slip_boundary_top(x):
    return (x[1] > (1.0 - tolerance))
def no_slip_boundary_back(x):
    return np.logical_and(x[0] > (3 - tolerance), (x[1]-0.5*1.0)**2 > 0.25**2)
def inlet_front_1(x):
    return np.logical_and(x[0] < tolerance, (x[1]-0.5*1.0)**2 <= 0.25**2)
def outlet_back_1(x):
    return np.logical_and(x[0] > (3 - tolerance), (x[1]-0.5*1.0)**2 <= 0.25**2)
def inlet_velocity_1(x):
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = (1/0.25)*(0.25**2 - (x[1]-0.5*1.0)**2)**0.5
    return values

facets_no_slip_boundary_front  = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_front)
facets_no_slip_boundary_bottom = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_bottom)
facets_no_slip_boundary_top    = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_top)
facets_no_slip_boundary_back   = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_back)
facets_inlet_front_1           = locate_entities_boundary(domain, domain.topology.dim - 1, inlet_front_1)
facets_outlet_back_1           = locate_entities_boundary(domain, domain.topology.dim - 1, outlet_back_1)

facet_indices = [facets_no_slip_boundary_front, facets_no_slip_boundary_bottom, \
                 facets_no_slip_boundary_top, facets_no_slip_boundary_back, \
                 facets_inlet_front_1, facets_outlet_back_1]

B, _      = W.sub(0).collapse()
u_no_slip = Function(B)
u_no_slip.x.array[:] = 0
bcu_no_slip_front  = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_front), W.sub(0))
bcu_no_slip_bottom = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_bottom), W.sub(0))
bcu_no_slip_top    = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_top), W.sub(0))
bcu_no_slip_back   = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_back), W.sub(0))

u_inlet_profile_1 = Function(B)
u_inlet_profile_1.interpolate(inlet_velocity_1)
bcu_inlet_1 = dirichletbc(u_inlet_profile_1, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front_1), W.sub(0))

G, _          = W.sub(1).collapse()
zero_pressure = Function(G)
zero_pressure.x.array[:] = 0
bcu_outlet_1  = dirichletbc(zero_pressure, locate_dofs_topological((W.sub(1), G), domain.topology.dim - 1, facets_outlet_back_1), W.sub(1))

bcs = [bcu_no_slip_front, bcu_no_slip_bottom, bcu_no_slip_top, bcu_no_slip_back,\
       bcu_inlet_1, bcu_outlet_1]

class Nonlinear_PDE_SNES_Problem:
    def __init__(self, F, J, soln_vars, bcs, P=None):
        self.L = F
        self.a = J
        self.a_precon = P
        self.bcs = bcs
        self.soln_vars = soln_vars

    def F_block(self, snes, x, F):
        assert x.getType() != "nest"
        assert F.getType() != "nest"
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
        with F.localForm() as f_local:
            f_local.set(0.0)

        offset = 0
        x_array = x.getArray(readonly=True)
        for var in self.soln_vars:
            size_local = var.vector.getLocalSize()
            var.vector.array[:] = x_array[offset : offset + size_local]
            var.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
            offset += size_local

        assemble_vector_block(F, self.L, self.a, bcs=self.bcs, x0=x, scale=-1.0)

    def J_block(self, snes, x, J, P):
        assert x.getType() != "nest" and J.getType() != "nest" and P.getType() != "nest"
        J.zeroEntries()
        assemble_matrix_block(J, self.a, bcs=self.bcs, diagonal=1.0)
        J.assemble()
        if self.a_precon is not None:
            P.zeroEntries()
            assemble_matrix_block(P, self.a_precon, bcs=self.bcs, diagonal=1.0)
            P.assemble()

w = Function(W)
(u, p)   = w.split()
(v, q)   = TestFunctions(W)
(du, dp) = TrialFunctions(W)

F_momentum = (
    rho*dot(dot(u,nabla_grad(u)),v)*dx
    + mu*inner(grad(u),grad(v))*dx
    - inner(p,div(v))*dx
)
F_continuity = (
    (q*div(u))*dx
)
F = form([F_momentum, F_continuity])
jacobian = form(
    [
        [derivative(F_momentum, u, du)  , derivative(F_momentum, p, dp)],
        [derivative(F_continuity, u, du), derivative(F_continuity, p, dp)],
    ]
)
P = form([[jacobian[0][0], None], [None, inner(dp,q)*dx]])

jacobian_matrix = create_matrix_block(jacobian)
P_matrix        = create_matrix_block(P)
residual_vector = create_vector_block(F)

options = PETSc.Options()
nonlinear_solver = PETSc.SNES().create(MPI.COMM_WORLD)
nonlinear_solver.getKSP().setType("fgmres")
nonlinear_solver.getKSP().getPC().setType("fieldsplit")
nonlinear_solver.getKSP().getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
nonlinear_solver.getKSP().getPC().setFieldSplitSchurFactType(PETSc.PC.SchurFactType.FULL)
nonlinear_solver.getKSP().getPC().setFieldSplitSchurPreType(PETSc.PC.SchurPreType.SELFP)

V_map = W0.dofmap.index_map
Q_map = W1.dofmap.index_map
offset_u = V_map.local_range[0]*W0.dofmap.index_map_bs + Q_map.local_range[0]
offset_p = offset_u + V_map.size_local*W0.dofmap.index_map_bs
is_u = PETSc.IS().createStride(V_map.size_local*W0.dofmap.index_map_bs, offset_u, 1, comm=PETSc.COMM_WORLD)
is_p = PETSc.IS().createStride(Q_map.size_local, offset_p, 1, comm=PETSc.COMM_WORLD)
nonlinear_solver.getKSP().getPC().setFieldSplitIS(("u", is_u),("p", is_p))

ksp_u, ksp_p = nonlinear_solver.getKSP().getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("hypre")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

nonlinear_solver.setFromOptions()
nonlinear_solver.getKSP().setFromOptions()

problem = Nonlinear_PDE_SNES_Problem(F, jacobian, [u, p], bcs, P)
nonlinear_solver.setFunction(problem.F_block, residual_vector)
nonlinear_solver.setJacobian(problem.J_block, jacobian_matrix, P_matrix)
x = create_vector_block(F)
u = Function(W0)
p = Function(W1)
with u.vector.localForm() as _u, p.vector.localForm() as _p:
    scatter_local_vectors(
        x,
        [_u.array_r, _p.array_r],
        [
            (u.function_space.dofmap.index_map, u.function_space.dofmap.index_map_bs),
            (p.function_space.dofmap.index_map, p.function_space.dofmap.index_map_bs),
        ],
    )
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
nonlinear_solver.solve(None, x)
PETSc.Sys.Print(
    f"Norm of x: {x.norm()}\n"
    f"Norm of u: {u.vector.norm(0)}\n"
    f"Norm of p: {p.vector.norm(0)}\n"
)
6

Hello, may I ask you where can I find the detailed documentation of PESTC ? The example codes provided by petsc are mostly written in C/C++, and they are very difficult to follow (most of them lack comments). It is difficult for me to convert them to python code and call them in Fenics.

7

I usually use a mixture of
https://petsc.org/release/petsc4py/reference/petsc4py.PETSc.html
and
https://petsc.org/release/manual/

1 Like