<div dir="ltr"><div dir="ltr">On Mon, Nov 8, 2021 at 11:00 AM barnafi <<a href="mailto:nabw91@gmail.com">nabw91@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi PETSc users,<br>
<br>
I am testing a block preconditioner for a nonlinear problem using <br>
fenics for matrix assembly, and I obtain correcet results in serial <br>
only (parallel diverges). The difficulty lies in using Nest matrices, <br>
as I want a performant code... I already see astonishing gains in <br>
serial.<br>
<br>
The code seems simple enough, where I have a UFL form from which I <br>
obtain all different blocks. Then, I create a SNESSolver as follows:<br>
<br>
<br>
--------------- SNESSolver.py -----------------<br>
class SNESProblem():<br>
def __init__(self, FFs, JJs, uu, pp, bbcs):<br>
self.Ls = FFs # UFL form for residuals<br>
self.Js = JJs # UFL form for jacobian<br>
self.bcs = bbcs<br>
self.bcs0 = [DirichletBC(_bc) for _bc in bbcs] # Homogenized <br>
boundary conditions for Newton iterations<br>
for bc in self.bcs0:<br>
bc.homogenize()<br>
self.u = uu<br>
self.p = pp<br>
<br>
def assemble_F(self, FF):<br>
Fu, Fp = FF.getNestSubVecs()<br>
assemble(self.Ls[0], tensor=PETScVector(Fu))<br>
assemble(self.Ls[1], tensor=PETScVector(Fp))<br>
for bc in self.bcs0:<br>
bc.apply(PETScVector(Fu))<br>
Fu.assemble()<br>
FF.assemble()<br>
<br>
def assemble_J(self, JJ):<br>
Juu = JJ.getNestSubMatrix(0, 0)<br>
Jup = JJ.getNestSubMatrix(0, 1)<br>
Jpu = JJ.getNestSubMatrix(1, 0)<br>
Jpp = JJ.getNestSubMatrix(1, 1)<br>
assemble(self.Js[0], tensor=PETScMatrix(Juu))<br>
assemble(self.Js[1], tensor=PETScMatrix(Jup))<br>
assemble(self.Js[2], tensor=PETScMatrix(Jpu))<br>
assemble(self.Js[3], tensor=PETScMatrix(Jpp))<br>
MM = Juu.getDiagonal().sum() / Juu.getSize()[0] # Better <br>
performance<br>
for bc in self.bcs0:<br>
indexes = [*bc.get_boundary_values().keys()]<br>
Juu.zeroRows(indexes, diag=MM)<br>
Jup.zeroRows(indexes, diag=0.0)<br>
Juu.assemble()<br>
Jup.assemble()<br>
JJ.assemble()<br>
<br>
def F(self, snes, xx, FF):<br>
<br>
xu, xp = xx.getNestSubVecs()<br>
# Update fenics vectors to keep track of ghost nodes<br>
xu.copy(self.u.vector().vec())<br>
xp.copy(self.p.vector().vec())<br>
self.u.vector().apply("")<br>
self.p.vector().apply("")<br>
self.assemble_F(FF)<br>
<br>
def J(self, snes, xx, JJ, PP):<br>
xu, xp = xx.getNestSubVecs()<br>
xu.copy(self.u.vector().vec())<br>
xp.copy(self.p.vector().vec())<br>
# Update fenics vectors to keep track of ghost nodes<br>
self.u.vector().apply("")<br>
self.p.vector().apply("")<br>
self.assemble_J(JJ)<br>
----------- SNESSolver.py end -------------<br>
<br>
I then initialize the Jacobian matrices to create the Nest (no 'isrows' <br>
nor 'iscols') and then initialize the PC from the SNES with<br>
<br>
------ PC init ------<br>
is_0, is_1 = J_mat.getNestISs()[0]<br>
ksp = snes.getKSP()<br>
pc = ksp.getPC()<br>
pc.setType("fieldsplit")<br>
pc.setFieldSplitIS(("0", is_0), ("1", is_1))<br>
---------------------<br>
<br>
The current behavior is that in serial everything works as expected, <br>
and instead in parallel not even the boundary conditions are set up <br>
correctly. The two possible problems I see regard the ordering of the <br>
dofs, or the update of ghost nodes, but I have tried modifying the code <br>
to account for this without success.<br>
<br>
If anyone has had similar issues, I would appreciate any comments. <br>
</blockquote><div><br></div><div>I am not adept enough at FEniCS to look at the code and debug. However, </div><div>I would make a very small problem for 2 processes, say two cells on each,</div><div>and then print out is_0 and is_1. That should tell you if you have the right</div><div>unknowns in each block. If so, then it seems like something is not correct in</div><div>your assembly routine, which is probably better debugged on the FEniCS list.</div><div>However, if they are not right, we can help get them right.</div><div><br></div><div> THanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Thanks in advance.<br>
<br>
Best,<br>
NB<br>
<br>
<br>
<br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>