Manual implementation of preconditioning and proper BC application

Question

Hi, I implemented a DG upwinding scheme for solution of incompressible NS equations in 2D/3D. The method leads to a saddle-point problem which I want to solve using a simple preconditioner that I would like to write myself.

If I would use in-build tools I would use among others the following (see Newton method programmed manually)

A,b = assemble_system(Jac,-T,bcs_hom)

in every Newton iteration. Basically, this line I want to rewrite manually, i.e. I have to assemble the system manually and also apply the homogenise BC.

You can see in the following code (sorry for the length, you can ignore most of the code if you want, it is there just in case you would like to run it) in the solve_SP() procedure how I assemble and solve the system

A*x + B*y = R1
BT*x = R2.

As I said, I'm basically interested just in the application of homogenise BC on this saddle-point system - I apply the homogenise BC on the matrix A and vector R1, but if I'm not wrong, I should apply somehow the homogenise BC even on BT and R2, or not? If so, how can I do it? Thanks a lot for any suggestions.

from dolfin import *
import petsc4py
petsc4py.init()
from petsc4py import PETSc
import numpy as np

dim = 2

# Mesh
######
Right   = 2.5
Top = 1.0
mesh    = RectangleMesh(0.0,0.0,Right,Top, 25, 10,"crossed")
bndry   = FacetFunction("size_t", mesh)

# Boundary description
######################
class LeftBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[0] < DOLFIN_EPS

class RightBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[0] > Right-DOLFIN_EPS

class TopBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[1] > Top-DOLFIN_EPS

class BottomBoundary(SubDomain):
    def inside(self, x, bndry):
        return x[1] < DOLFIN_EPS

# Function spaces
#################
V = VectorFunctionSpace(mesh, 'DG', 1)
Q = FunctionSpace(mesh, 'DG', 0)

# Boundary conditions
#####################
noslip    = Expression(("0.0", "0.0"), t = 0.0)
u_IN      = Expression(("x[1]*(1.0-x[1])", "0.0"), t = 0.0)

bc_top    = DirichletBC(V, noslip, TopBoundary(), 'pointwise')
bc_bottom = DirichletBC(V, noslip, BottomBoundary(), 'pointwise')
bc_left   = DirichletBC(V, u_IN, LeftBoundary(), 'pointwise')

bcs     = [bc_top, bc_bottom, bc_left]
bcs_hom = homogenize(bcs) 

# Unknown and test functions
############################
u   = Function(V)
u_inc   = Function(V) # incremental function for Newton process
p   = Function(Q)
p_inc   = Function(Q) # incremental function for Newton process

v   = TestFunction(V)
q   = TestFunction(Q)

u_1 = Function(V) # velocity from previous time step
du  = TrialFunction(V)
dp  = TrialFunction(Q)

# Problem parameters
####################
T_end       = 0.02
dt      = 0.01
t       = dt # time variable
nu      = 1.36*1.e-4 # air viscosity

u0      = Expression(("0.0", "0.0"))
f_rhs       = Expression(("0.0", "0.0"), t = 0.0)

sigma       = 10.0 # penalty parameter in J() form

atol, rtol = 1e-7, 1e-10      # abs/rel tolerances
lmbda      = 1.0              # relaxation parameter
maxIter    = 10 

# Variational formulation
#########################
n   = FacetNormal(mesh)
ds  = Measure("ds", subdomain_data = bndry)
I   = Identity(u.geometric_dimension())
D   = 0.5*nu*(grad(u) + grad(u).T)

def a(u,v): 
    M = inner(D,grad(v))*dx - inner(avg(D)*n('+'),jump(v))*dS 
    return M

def J(u,v):
    M = sigma*inner(jump(u),jump(v))*dS
    return M

def b(p,v):
    M = -p*div(v)*dx  + avg(p)*dot(jump(v),n('+'))*dS
    return M

def c(u,w,v):   
    P = avg(dot(w,n))
    H = conditional(P < 0.0, dot(u('+'),w('+')), dot(u('-'),w('-')))
    M = -0.5*inner(grad(v)*u,w)*dx + inner(0.5*H*n('+'),jump(v))*dS
    return M

def L(v):
    M = inner(f_rhs,v)*dx 
    return M

r1  = (1/dt)*inner(u - u_1,v)*dx + a(u,v) + J(u,v) + c(u,u_1,v) - L(v) + b(p,v)
r2  = b(q,u)

a_d     = derivative(r1,u,du)
b_d = derivative(r1,p,dp)
b_dT    = derivative(r2,u,du)

# Solution
##########

plt = plot(u, mesh = mesh, mode="color", interactive = True, wireframe = False)
pltP = plot(p, mesh = mesh, mode="color", interactive = True, wireframe = False)

u_1.assign(interpolate(u0,V)) # apply initial condition
for bc in bcs:
    bc.apply(u.vector())  # apply boundary condition


def solve_SP():
    # Solution vectors
    x = Function(V)
    y = Function(Q)

    # Assembling Saddle-point blocks
    A = PETScMatrix()
    assemble(a_d, tensor=A)
    for bc in bcs_hom:
        bc.apply(A) 

    B = PETScMatrix()
    assemble(b_d, tensor=B)

    BT = PETScMatrix()
    assemble(b_dT, tensor=BT)   

    # Assembling RHS for Newton
    R1 = assemble(-r1)
    R2 = assemble(-r2)
    for bc in bcs_hom:
        bc.apply(R1)

    # Solver
    solver = LUSolver(A)
    solver.parameters['reuse_factorization'] = True

    # solving A*x_ = R1 
    x_ = Function(V)
    solver.solve(x_.vector(),R1)

    # solving A*xv = b, where xv are solution vectors and b is i-th column of B
    n = V.dim()
    m = Q.dim()
    xv, b = [Vector(None, n) for i in range(2)]
    ei = Vector(None,m) 
    # Matrix of solutions x as columns
    K = PETSc.Mat().createDense([n, m])
    K.setUp()
    rows = range(n)
    for i in range(m):
        # This is i-th R^m basis vector
        ei.zero()
        ei[i] = 1
        # B.ei gives the i-th colum
        B.mult(ei, b)
        solver.solve(xv, b)
        # Set the column
        col = [i]
        K.setValues(rows, col, xv.array())
    K.assemblyBegin()
    K.assemblyEnd()
    K = PETScMatrix(K)

    # Solving BT*K*y = BT*x_ - R2
    BTK = PETSc.Mat().createDense([m, m])
    BTK.setUp()
    iCol = Vector(None,n)
    iBTK = Vector(None,m)
    rows = range(m)
    for i in range(m):
        # This is i-th R^m basis vector
        ei.zero()
        ei[i] = 1
        # K.ei gives the i-th column iCol
        K.mult(ei, iCol)
        # BT.iCol gives i-th column of BT*K
        BT.mult(iCol,iBTK)
        # Set the column
        col = [i]
        BTK.setValues(rows, col, iBTK.array())
    BTK.assemblyBegin()
    BTK.assemblyEnd()
    BTK = PETScMatrix(BTK)  
    # Finalizing BT*K*y = BT*x_ - R2
    rhs = Vector(None,m)
    BT.mult(x_.vector(),rhs)
    rhs = rhs - R2
    solve(BTK,y.vector(),rhs,'gmres','amg') 

    # Computing K*y in order to set x = x_ - K*y    
    Ky = Function(V)
    K.mult(y.vector(),Ky.vector())
    x.vector()[:] = x_.vector() - Ky.vector()

    # COMPUTE ERROR
    err = np.sqrt(x_.vector().norm('l2')**2 + rhs.norm('l2')**2)
    return [x,y,err]

while t<=T_end: 
    print 'time step: ', t
    nIter = 0
    res = 1
    rel_res = 1
    while res > atol and rel_res > rtol and nIter < maxIter:     # Newton iterations
        nIter += 1
        [u_inc,p_inc,res] = solve_SP()
        if (nIter == 1):
            rel_res = 1
        else:
            rel_res = res/res_1
        res_1 = res
        string = "Newton iteration %d: r (abs) = %.3e (tol = %.3e) r (rel) = %.3e (tol = %.3e)"
        print string % (nIter, res, atol, rel_res, rtol)                            
            u.vector()[:] += lmbda*u_inc.vector()    # New u vector
            p.vector()[:] += lmbda*p_inc.vector()    # New p vector
    u_1.assign(u)
    t += float(dt)

plt.plot(u)
pltP.plot(p)

PS: I'm comparing the results with results that I get using the in-build NewtonSolver (with GMRES and ILU precond.) and I'm getting different results even though convergence is nice.

Comments

Hi, for applying bcs to matrices, have a look at how the DirichletBC class is implemented. In particular methods zero, zero_columns and apply are relevant. If the symmetry is not important you now only need to modify the B matrix. Since all your bcs are on the V space you have to work with V row of the block matrix. You have already changed A and R1.

---

Thanks a lot for a hint. I'll look at it more thoroughly in a few days. Let me in a meantime add a few remarks/questions:
1) Is it possible to "mix" python and C++ code or should I make a library in C++ and just import it?
2) I should apply bcs even on BT and R2, right? Since I have zero BC on V, it should affect columns in BT and R2 (i.e. make them zero), or not?
3) It seems to me now that main difficulty is to identify rows/columns that are associated with basis functions whose coefficients should be zero.
4) I'm playing along with cdc.block library for an alternative solution.

---

Hi, 2) if the problem is symmetric it is desirable to keep that property also after the bcs are applied. In that case you also modify the BT and R2. This is what cbc.block does if you pass symmetric=True to block_bc constructor. For 3) have a look at DirichletBC.get_boundary_values. 4) cbc.block will make your life a lot easier. 1) is a matter of personal preference.

---

Hi, I wasn't able to find much time last month to solve the problem. Nevertheless, I recently encounter a few problems and would be grateful for a help.
1) It turns out, that cbc.block's routine 'block_assemble' works only with forms made of TrialFunctions and fails for forms made of Functions. Since I need to work with classical Functions, it sucks ...
2) Can you give some example on how to correctly import a self-made c++ library or even how to combine a python and c++ code in one .py file?
Thanks a lot for any help.