How to compute the residual of a non-linear problem


38
views
0
8 days ago by
ac  
Hi everybody,

I'm trying to compute transient non-linear heat equation. However, to simplify implementation and spare cpu time, I evaluate the specific heat on the previous time step, so that at each time step i solve a linear problem. I'm trying to print the residual to check the error I made by using the command residual. But I face an error I'm not understanding. Can you help me fix this ? Is it possible to call residual outside a nonliearproblem, as I'm trying to do ?
The code looks this:

from fenics import *
file = File("poisson.pvd")
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
u_=Function(V) # temperature field at the previous time step
usol = Function(V) # temperature field at the current time step
u_.assign(usol)
v = TestFunction(V)

t=0.
tf=1.
dt=tf/10.

a = (1+u_)/dt* u *v*dx + inner(grad(u), grad(v))*dx
L = (1+u_)/dt* usol*v*dx + v*dx

a_res = (1+usol)/dt* u *v*dx + inner(grad(u), grad(v))*dx
L_res = (1+usol)/dt* usol*v*dx + v*ds

# Compute solution
while t<tf:
	t+=dt
	solve(a == L, usol, bc)
	print(residual(a_res,usol,L_res))
	u_.assign(usol)
# Save solution in VTK format
	file << (u_,t)

​

The error I get:

TypeError: in method 'residual', argument 1 of type 'dolfin::GenericLinearOperator const &'

Thanks for the help. kind regards,

Community: FEniCS Project

2 Answers


0
8 days ago by
The residual() function is for evaluating residuals of linear algebraic systems, rather than nonlinear variational problems.  

I think you may be looking for something like the following modified version of your code (although I have done some interpretation, so read carefully):

from fenics import *
file = File("poisson.pvd")
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
u_=Function(V)
usol = Function(V)
u_.assign(usol)
v = TestFunction(V)

t=0.
tf=1.
dt=tf/10.

def specificHeat(u):
    return 1.0 + u

def resFunc(u, u_, u_specificHeat,v):
    return specificHeat(u_specificHeat)*((u-u_)/dt)*v*dx \
        + inner(grad(u),grad(v))*dx \
        - v*ds  # Note: MWE uses ds here in one residual and dx in the other

#a = coeff(u_)/dt* u *v*dx + inner(grad(u), grad(v))*dx
#L = coeff(u_)/dt* u_*v*dx + v*dx

resLinear = resFunc(u,u_,usol,v)
a = lhs(resLinear)
L = rhs(resLinear)

#a_res = (1+usol)/dt* u *v*dx + inner(grad(u), grad(v))*dx
#L_res = (1+usol)/dt* usol*v*dx + v*ds

# Compute solution
while t<tf:
    t+=dt
    solve(a == L, usol, bc)
    
    #print(residual(a_res,usol,L_res))
    
    # Simplest code:
    resVec = assemble(resFunc(usol,u_,usol,v))
    bc.apply(resVec)
    print(norm(resVec))
    
    # Alternative using residual(), producing the same results to
    # machine precision:
    A_res = assemble(a) # (Assembled with updated usol in specific heat)
    B_res = assemble(L)
    bc.apply(A_res,B_res)
    print(residual(A_res,usol.vector(),B_res))
    
    u_.assign(usol)
    # Save solution in VTK format
    file << (u_,t)
0
8 days ago by
ac  
Thanks for the quick and useful answer. This is exactly what I was looking for.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »