Time-dependent nonlinear problem in 3D

Question

Hi guys,
I am trying to solve the heat equation with (nonlinear) radiation and convection from the surfaces in 3D.
In the file demo_cahn-hilliard.py they are beginning with defining two classes - the initial condition (which is more easy to me for I just want to have a constant initial temperature) and the class for interfacing with the Newton solver.

My problem is, I dont really know how to define the class for my initial condition. What I wrote is:

u0 = Constant(293)
class InitialCondition(Constant):
   def __init__(self):
      u0
   def eval(self, values, x):
      values=u0

But I can't interpolate with this conditions:

mesh = UnitCubeMesh(10,10,10)
V = FunctionSpace(mesh, 'Lagrange', 1)
ME = V*V
du = TrialFunction(ME)
v = TestFunction(ME)
u = Function(ME)
u_1 = Function(ME)
u_init = InitialCondition()
u.interpolate(u_init)
u_1.interpolate(u_init)

I get an error message now:

TypeError: in method 'Function_interpolate', argument 2 of type 'dolfin::Generic Function const &'

What should I do?
Your help would be highly appreciated.

EDIT: I was able to fix this somehow. And I reduced the functionspace to:

ME = V

instead of

ME=V*V

because I dont have a mixture of 2 fluids but a solid body of one material.

But a new error messages appeared:

ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data=''), 1, quad_scheme=None), 0, None), Argument(FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data=''), 1, quad_scheme=None), 1, None)) vs (Argument(FiniteElement('Lagrange', Domain(Cell('tetrahedron', 3), label='dolfin_mesh_with_id_0', data=''), 1, quad_scheme=None), 0, None),).

Answer 1

Hi,

The problem is that you try to interpolate a function defined on a mixed space (V*V) on a single function space (such as V for example), here's a fix that should work:

from dolfin import *
u0 = Expression('293')
class InitialCondition(Expression):
   def __init__(self):
      u0
   def eval(self, values, x):
      values[0]=293
      values[1]=0.0
   def value_shape(self):
        return (2,)

mesh = UnitCubeMesh(10,10,10)
V = FunctionSpace(mesh, 'Lagrange', 1)
ME = V*V
du = TrialFunction(ME)
v = TestFunction(ME)
u = Function(ME)
u_1 = Function(ME)
u_init = InitialCondition()
u.interpolate(u_init)
u_1.interpolate(u_init)

Here I set the values of the second space to be 0.0 but you can change it if you want.

Hoping this will help you,

Best regards.

Comments

Seems like the problem were in the radiation terms; I cant find the problem though. I tried uuu*u as well as u**4, but this error message keeps appearing.
The convection terms work in the linear solver btw.

---

Hi, sorry for the delay, but I found some time to look to the code today.

First, I think that the problem is not well posed. In the Cahn Hilliard demo the linear form (right hand side of the variational form a(u,v)=L(v) ) is set first, then the left hand side is calculated as the derivative of the right hand side so as to build the following non linear variational form : J(u,v) = L(v) where J is the jacobian of L (calculated with the derivative function). The Newton method must be built this way, so as to solve the following equation: derivative(L) = -L.

In order to obtain this equation with your specific problem, the derivative function must be applied to the linear form (with only test functions) so I think we should apply derivative to FL only (Some changes have to be made to write it correctly though...)

There are also some little things that made the code crash when I tried to run it but it may only be related to differences between our versions of FEniCS, I'm running version 1.5, I tried to make the corrections visible. Still, as far as I can see the problem should be well posed now but something doesn't work and I think it may be related to the definition of FL but I'm not much in heat equation and I don't want to write something that would not correspond to your problem. Maybe you could tell me the exact equation you are trying to solve so that we can think about how to write it in a Newton suitable way.

Here's the code with some corrections.

from dolfin import *
import numpy

class HeatEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A) 

parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"

mesh = UnitCubeMesh(10,10,10)
V = FunctionSpace(mesh, 'Lagrange', 1)
'''Define boundary parts:'''

boundary_parts = FacetFunction("size_t", mesh)
boundary_parts.set_all(0)

tol = 1E-14

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[2] - 1) < tol
top=Top()
top.mark(boundary_parts, 1)

du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
u_1 = Function(V)

u0 = Constant(293)
source = Constant(10)
del_t = 0.1
t_ges = 1
u=interpolate(u0,V)
u_1=interpolate(u0,V)

ds = Measure('ds')[boundary_parts]

# dx(mesh)=dx(), dx() is the volumic integral over th whole mesh

# del_t*u0*v*ds(1)-->FL, this part of Fa should belong to FL

#Fa = (du*v + del_t*inner(grad(du),grad(v)))*dx() - del_t*du*v*ds(1) 

FL = u0*v*dx() - del_t*source*v*ds(1) + del_t*u0*v*ds(1)

Fa = derivative(FL, u, du) #derivating FL so as to obtain Fa

problem = HeatEquation(Fa,FL)

solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

'''Time dependency'''
t = del_t
while t <= t_ges:
    print 'Zeit =', t
    u_1.vector()[:] = u.vector()
    solver.solve(problem, u.vector()) #crashes here : "All terms in form must have same rank."
    plot(u, key='u')
    #u_1.assign(u)

interactive()

---

Finally I made it work and I think it is correct, after some reformulation of the form in Fa AND FL. Please scratch the code written above, here's the working one, don't hesitate to ask if some details are needed:

from dolfin import *
import numpy
#-------------------Class for interfacing with Newton Solver----------------#
class HeatEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A) 
#-------------------------------Parameters----------------------------------#
u0 = Constant(293)
source = Constant(10)
del_t = 0.1
t_ges = 1
#-------------------------Form Compiler options-----------------------------#
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"
#----------------------Create mesh and function space-----------------------#
mesh = UnitCubeMesh(10,10,10)
V = FunctionSpace(mesh, 'Lagrange', 1)
#-------------------------Create subdomains for BCs-------------------------#
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[2] - 1) < DOLFIN_EPS

boundary_parts = FacetFunction("size_t", mesh)
boundary_parts.set_all(0)

top=Top()
top.mark(boundary_parts, 1)

ds = Measure('ds')[boundary_parts]
#-------------------------Define Trial and Test functions-------------------#
du = TrialFunction(V)
v = TestFunction(V)
#-----------------------------Define Functions------------------------------#
u = Function(V) #Current Solution
u_1 = Function(V) #Solution from previous converged step
#-------------------Interpolate for initial conditions----------------------#
u=interpolate(u0,V)
u_1=interpolate(u0,V)
#----------------------Weak statement of the equation-----------------------#
L0 = (u*v + del_t*inner(grad(u),grad(v)))*dx() - del_t*u*v*ds(1) 
L1 = u_1*v*dx() - del_t*source*v*ds(1) + del_t*u_1*v*ds(1)
L = L0 - L1
#--Compute directional derivative about u in the direction of du (Jacobian)-#
a = derivative(L, u, du)
#---------------Create nonlinear problem and Newton solver------------------#
problem = HeatEquation(a,L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
#------------------------------Time dependency------------------------------#
t = del_t
while t <= t_ges:
    print 'Zeit =', t
    u_1.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    plot(u, key='u')
    #u_1.assign(u)
interactive()

EDIT: The output looks good but the values are odd, don't know how to fix it but it is due to a lacking term in the equation I think...

---

Thank you! I am sorry for the late answer. The values must be odd because I didnt write all the material parameters here to keep the code at least a little bit clear.

---

you forgot to use "t+=del_t" in your loop. Do it or it will never finish.