Assemble a nonlinear form

General
1

Hello Everyone. I am trying to implement a scheme for solving the obstacle problem described at this link.

The code for the elliptic problem I am interested in is available below:

from fenics import *
from obstacles import dome
from mshr import mesh
def smoothmax(r, eps=1e-4): 
	return conditional(gt(r, eps), r-eps/2, conditional(lt(r, 0), 0, r**2/(2*eps)))
mesh = generate_mesh(Circle(Point(0, 0), 2), 25)
V = FunctionSpace(mesh, "CG", 1)
w = Function(V)
v = TestFunction(V)
psi = Constant(0.)
f = Constant(-1)		
eps = Constant(pow(10, -8))
bc = DirichletBC(V, 0.8, "on_boundary")
F = dot(grad(w), grad(v))*dx - 1/eps*inner(smoothmax(-w+psi), v)*dx - f*v*dx
solve(F == 0, w, bcs=bc)
vtkel = File("output/el.pvd")
vtkel << w

Instead of using the Newton method as in the snippet above, I would like to assemble the form F using something similar to assemble() so that I can multiply this matrix by a vector (I need this because I would like to test an iterative scheme).

Would anyone know a way to solve this problem?
Thanks in advance.

2

See for instance Custom Newton solvers — FEniCSx tutorial
which goes through what you would need to create your own Newton/iterative solver.
You could also consider: Implementing the Newton Method per hand - #3 by j-507
among one of many different posts on the subject

1 Like
3

Thank you for the suggestion. I am not actually using Fenicsx yet. The scheme I tried to implement is the following

from fenics import *
from obstacles import dome
from mshr import mesh
def smoothmax(r, eps=1e-4): 
	return conditional(gt(r, eps), r-eps/2, conditional(lt(r, 0), 0, r**2/(2*eps)))
mesh = generate_mesh(Circle(Point(0, 0), 2), 25)
V = FunctionSpace(mesh, "CG", 1)
w = Function(V)
v = TestFunction(V)
u=TrialFunction(V)
psi = Constant(0.)
f = Constant(-1)		
eps = Constant(pow(10, -8))
bc = DirichletBC(V, 0.8, "on_boundary")
l=0
error=5.0
while error>1.0e-06:
    if l==0:
        a=dot(grad(u), grad(v))*dx
        L=f*v*dx
        solve(a==L,w,bc)
    else:
        wk=w.copy()
        wnew=Function(V)
        a=u*v*dx
        F = dot(grad(wk), grad(v))*dx - 1/eps*inner(smoothmax(-wk+psi), v)*dx - f*v*dx
        # Brezis-Sibony scheme
        solve(a==F, wnew, bcs=bc)
        error=errornorm(wnew,w,norm_type='L2')
        w.assign(wnew)
vtkel = File("output/el.pvd")
vtkel << w

However, when I run the solver I get an error message saying that the type of finite element I chose is not correct and the error comes from plugging w into the variational formulation.

I would like to use the function w that I found at the first iteration of the loop as a given datum in the successive iterations.
Would there be a way to do this?

EDIT. The code should now work. The trick was to set up wk=w.copy().