Lumped mass scheme and Picard iteration--fails with EPS=0.0

1

Hi, everyone. I’m trying to solve the Richards’ equation with lumped mass scheme and Picard iteration. The variational problem can be expressed as:

(u-u0)*v*dx-dt*D*dot(grad(u),grad(v))*dx=0

and following is my core code with reference to Mass lumping in time dependent problem and Picard Iteration Not Updating Properly

from dolfin import *
import numpy

#define time step
dt=60
time_steps=120

#define mesh and function space
mesh=UnitSquareMesh.create(1,100,CellType.Type.quadrilateral)
V=FunctionSpace(mesh,'CG',1)

#define boundary condition
u_D=Expression('0.41',degree=0)
def boundary(x):
    return 1-x[1]<DOLFIN_EPS
bc=DirichletBC(V,u_D,boundary)

#define initial value u0 and known uk
u_init=Expression('0.065',degree=0)
u0=project(u_init,V)
uk_ex=Expression('0.2',degree=0)
uk=project(uk_ex,V)

#define lumped mass matrix
u=TrialFunction(V)
v=TestFunction(V)
mass_form = u*v*dx
mass_action_form = action(mass_form, Constant(1))
M_lumped = assemble(mass_form)
M_lumped.zero()
M_lumped.set_diagonal(assemble(mass_action_form))

#define variational problem
def D(x,tol):     #a function to calculate D
    Se=(x-0.065)/0.345+tol
    K=1.22801e-5*Se**0.5*(1-(1-Se**2.123)**0.471)**2
    B=0.345/(x-0.065+tol)+tol
    return K*0.052/(x-0.065+tol)**2*B**1.123*(B**2.123-1)**0.89 

D=D(uk,1e-8)
A=M_lumped-assemble(dt*D*inner(grad(u),grad(v))*dx)
b=M_lumped*u0.vector()
bc.apply(A,b)

#Picard iteration
u=Function(V)
output=File('Richards/solution.pvd')
for t in range(0,time_steps):
    iter=0
    eps=1
    tol=1e-6
    converged=False
    while not converged:
        solve(A,u.vector(),b)
        diff=u.vector().get_local()-uk.vector().get_local()
        eps=numpy.linalg.norm(diff,ord=numpy.Inf)
        uk.assign(u)
        print(eps)
        converged=eps<tol
        iter+=1
    output << u
    print('EPS: '+str(eps)+'; Itrations: '+str(iter))
    u0.assign(u)

However, this code finally leads to a constant EPS of 0.0 for each iteration with iteration steps staying for 1. It seems that the Picard iteration fails to iterate but I can’t figure out where went wrong. Any help would be appreciated, thank you!

2

@dokken dear Dokken, do you have any suggestions?

3

It is the Matrix A and b should be updated in time stepping but not only uk and u0. Thus the definition of A and b should be included into the iteration loop in this problem for solution. Thanks anyway.