Symmetric DirichletBC?

1

Hi

I am solving nonlinear solid mechanics computations. So I overrides dolfin.NonlinearProblem to apply assembly functionalities as well as applying DirichletBC using bc.apply

When I check if the system Jacobian matrix is symmetric (A.mat().isSymmetric(...)), it says no. Yet the formulation tells the symmetry. It happens to the case of solving nonlinear Poisson’s equation.

Does anyone have an idea how to get DirichletBC applied without breaking symmetry? Or is it just me missing something here?

Thanks,
Victor

2

One way is to use assemble_system, as in the following example:

from dolfin import *
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh,"CG",1)
bc = DirichletBC(V,Constant(0.0),"on_boundary")
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u),grad(v))*dx
L = Constant(1.0)*v*dx

# The `assemble_system` function maintains symmetry:
A,B = assemble_system(a,L,bcs=[bc])
print(as_backend_type(A).mat().isSymmetric())

# Using `assemble` followed by `apply` does not:
A = assemble(a)
bc.apply(A)
print(as_backend_type(A).mat().isSymmetric())
3

Thanks, will try this out to see how i can incorporate it in subclass of NonlinearProblem

4

Didn’t get any luck when using assemble_system along with NonlinearProblem for nonlinear problems. A very sloppy trial for idea testing is:

class MyNP(NonlinearProblem):
    def __init__(self, f, j, bcs):
        self.A = PETScMatrix()
        self.b = PETScVector()
        self.f = f
        self.j = j
        self.bcs = bcs
    
    def F(self, b, x):
        assemble_system(A_form=self.j, b_form=self.f, bcs=self.bcs, A_tensor=self.A, b_tensor=b)
    
    def J(self, A, x):
        assemble_system(A_form=self.j, b_form=self.f, bcs=self.bcs, A_tensor=A, b_tensor=self.b)

Unfortunately, it doesn’t work that residual just goes up iteration by iteration until nan come up and stops the simulation.

Any thing wrong?
Thanks
Victor

5

You can use the optional x0 argument of assemble_system to set BCs in a way consistent with the system for a Newton step:

from dolfin import *

mesh = UnitSquareMesh(100,100)
V = FunctionSpace(mesh,"Lagrange",1)
u = Function(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
f = x[0]*x[1]
res = (1.0+u*u)*inner(grad(u),grad(v))*dx - inner(f,v)*dx
Dres = derivative(res,u)

expr = Expression("sin(10.0*x[0]) + x[1]*x[0]",degree=3)
bc = DirichletBC(V,expr,"on_boundary")

class CustomNonlinearProblem(NonlinearProblem):
    def form(self,A,P,b,x):
        assemble_system(Dres,res,A_tensor=A,b_tensor=b,bcs=[bc,],x0=x)
    def F(self,b,x):
        pass
    def J(self,A,x):
        pass

problem = CustomNonlinearProblem()
solver = PETScSNESSolver()
solver.solve(problem,u.vector())

from matplotlib import pyplot as plt
plot(u)
plt.show()
1 Like
7

After clearing the misunderstanding of the symbol x that you used at different spots for different meanings, this works as a charm!

Thanks!
Victor

8

One question tho: when using the CustomNonlinearProblem along with dolfin.PETScSNESSolver, it seems to be called multiple times. To see that, I modified form:

def form(self, A, P, b, x):
    print('assembly starts')
    assemble_system(...)
    print('assembly ends')

I saw multiple times of these two lines coming out when using PETScSNESSolver, but not with NewtonSolver. I suspect that assembly of jacobian then happens multiple times per iteration, you know a resolution to avoid this?

Thanks
Victor

9

Hmm… Looking into the code here

https://bitbucket.org/fenics-project/dolfin/src/master/dolfin/nls/

it’s clear that form is being called before each of J and F (and not just once, before both), although I don’t see any obvious workaround. (Trying to control calling of assemble_system using a Boolean that is switched off by J and on again by F doesn’t seem to work.)

10

What I was thinking was for the first time form being called, I can assemble both A and b.

Case 1
If the first call is for A, then I assemble with self.b, and when calling for assembling b, I can directly assign b with self.b.

Case 2
If the first time is for b, then I assemble with self.A and the second call, I assign A with self.A. The reason is assembly of Jacobian for my case seems not to be trivial.

But I don’t see a way of doing for the Case 2. If this idea sounds reasonable, do you know how?

Thanks
Victor

11

Hi Victor,

I’m not sure that you’re still working on this problem, but something that may be of use to you (or future readers), when trying to avoid reassembly of the A matrix, is to switch from using the assemble_system function, and instead use a SystemAssembler. The SystemAssembler class is designed to apply boundary conditions at assembly (thus maintaining the symmetry in A), but can assemble b or A independently or simultaneously. You initialize an assembler as:

assembler = SystemAssembler(Dres,res,bcs=[bc,])

Then, if you want to assemble both A and b simultaneously, you can use:

assembler.assemble(A,b,x0)

or if you just want to reassemble the right-hand side (without having to re-assemble A):

assembler.assemble(b,x0)

Cheers,
Daniel

12

Hi Daniel

Ah, I c. Will try when having a chance. Thanks!