How to compute the residual of a non-linear problem

38
views
0
8 days ago by
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.

L = (1+u_)/dt* usol*v*dx + 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

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 \
- v*ds  # Note: MWE uses ds here in one residual and dx in the other

#L = coeff(u_)/dt* u_*v*dx + v*dx

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

#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
Thanks for the quick and useful answer. This is exactly what I was looking for.