### How to compute the residual of a non-linear problem

38

views

0

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:

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

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):

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

Please login to add an answer/comment or follow this question.