Nonlinear Poisson tutorial

[bhaveshshrimali]

Hi @jorgensd,

Although a pretty minor detail, but isn't scatter_forward repeated in the Nonlinear poisson tutorial at the end? NewtonSolver.solve() automatically applies scatter_forward to the solution, as I can see here?

On another note, do you feel that adding an example on, say, implementing the Newton's method by hand would add value to tutorial? This could be illustrative in use cases where one needs to implement their own nonlinear solve? Something like this for the nonlinear poisson

from ufl import TrialFunction, derivative
Jac = derivative(F, uh, TrialFunction(V))

bilinear_form = fem.form(Jac)
linear_form = fem.form(-F)

du = fem.Function(V, name="incr_u")

bilinf = fem.assemble.create_matrix(bilinear_form)
linf = fem.assemble.create_vector(linear_form)

# zero the entries
bilinf.zeroEntries()

fem.assemble_matrix(bilinf, bilinear_form, bcs=bcs)
bilinf.assemble()

with linf.localForm() as loc:
    loc.set(0.)

fem.assemble_vector(linf, linear_form)

fem.apply_lifting(linf, [bilinear_form], [bcs])
linf.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(linf, bcs)

res_0 = linf.norm()
res = res_0
print(f"Initial residual norm: {res_0}")
solver = PETSc.KSP().create(mesh.comm)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.setOperators(bilinf)
solver.solve(linf, du.vector)
du.x.scatter_forward()
uh.vector.axpy(1.0, du.vector)
iters = 0
maxIter = 30

while res/res_0 > 1.e-8 and iters < maxIter:

    
    bilinf.zeroEntries()
    fem.assemble_matrix(bilinf, bilinear_form, bcs=bcs_hom)
    bilinf.assemble()

    with linf.localForm() as loc:
        loc.set(0.)

    fem.assemble_vector(linf, linear_form)
    fem.apply_lifting(linf, [bilinear_form], [bcs_hom])
    linf.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    fem.set_bc(linf, bcs_hom)
    solver.setOperators(bilinf)
    solver.solve(linf, du.vector)
    du.x.scatter_forward()
    uh.vector.axpy(1.0, du.vector)
    res = linf.norm()
    print(f"norm: {res}")
    iters += 1

instead of the call to NewtonSolver.

[jorgensd]

I've made a patch removing the unnecessary need for ghost updates in that demo.
I'll think about adding a custom newton solver.
It would have to cover a slightly different aspect than the non-linear demo.

My idea would be:
-To introduce it as a 1D "projection" method, i.e. $ u + u^2 = f$ in $\Omega$, $u= \frac{-1 \sqrt{1+4f}}{2}$ on $$\partial\Omega$$.

• Then it would go through what a Newton solver does for a simple function of one parameter (without bcs).
• Next describe how it works for variational forms, ending up with matrix F and residual vector J.
• Next describe application of DirichletBC to the initial guess, meaning that one can use a homogenous Dirichlet for correcting this behavior
• Solve the problem as described in your algorithm above

[jorgensd]

Fixed in #85