How to Impliment a Time Dependent Non-Linear PDE Solver

1

Hi everyone,

I am very new to Fenics and I am currently trying to model the heating of a metal block in 2D via a laser including conduction and convection in order to find a steady state temperature solution for the whole block. As a result, my resultant weak form equation looks like:

F = T*v*dx + dot(grad(T), grad(v))*dx*dt + dt*(1/kappa)*Q*v*ds(3) + (1/kappa)*sum(integrals_R)-T_n*v*dx

where integrals_R is defined via:

boundary_conditions = {0: {'Robin':     (es, T_a4)},
                       1: {'Robin':     (es, T_a4)},
                       2: {'Robin':     (es, T_a4)},
                       3: {'Robin':     (es, T_a4)}}


ds = Measure('ds', domain=mesh, subdomain_data = boundary_markers)
integrals_R = []
T = TrialFunction(V)
v = TestFunction(V)
T_ = Function(V)
T_n = interpolate(T_am, V)


#Account for Radiation
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        integrals_R.append(es*dt*(T**4 - T_a4)*v*ds(i))

This dictionary section may not be necessary I was originally using this for Robin boundary conditions, which is not what I have here. The first issue I have is with the multiplication of ds(i) and dt which gives me an error:

  File "test2D.py", line 110, in <module>
    F = T*v*dx + dot(grad(T), grad(v))*dx*dt + dt*(1/kappa)*Q*v*ds(3) + (1/kappa)*sum(integrals_R)-T_n*v*dx
TypeError: unsupported operand type(s) for *: 'Form' and 'float'

I am unsure of how I can rectify this error.

In addition, I am aware that I need to solve this problem in a Non-linear solver, but I am unsure of how to do this. I have written some code down but I am unsure if this is on the right lines or not.

my full code follows:

rom __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import math

time = 100.0            # final time
num_steps = 100     # number of time steps
dt = time / num_steps # time step size


#build rectangular mesh to model 2D slice
nx = 32
ny = 16
x0 = 0
y0 = 0
x1 = 4
y1 = 2
pi = 3.141592653589793
T_am = Constant("300")
T_a4 = T_am**4
epsilon = 1  # material emissivity
sigma = 5.67E-8 # W/(m2.K4)
es = epsilon*sigma
Q = Expression('2*A*t*Pmax/(pi*w*w)*exp(-2*pow(x[0]-2, 2)/(w*w))',degree=2, A=1, Pmax=200,w=1,x0=0, y0=0,t=0)
#----------THERMAL--PROPERTIES----------
kappa = 94         #conductivity [W/m K] 
c = 0.444          #Specific Heat Capacity [J/gK]
rho = 7            #Density [g/cm^3]
#---------------------------------------
 
sw = Point (x0, y0)
ne = Point (x1,y1)

mesh = RectangleMesh(sw, ne, nx, ny, diagonal= "right")
V = FunctionSpace(mesh, 'P', 2)



#-----BOUNDARY-CONDITIONS--------------------------


#-----Define Markers for the different parts of the boundary-------
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
tol = 1E-14
class BoundaryX0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x0, tol)
    
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0)

class BoundaryX1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], x1, tol)
    
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1)

class BoundaryY0(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y0, tol)
    
by0 = BoundaryY0()
by0.mark(boundary_markers, 2)

class BoundaryY1(SubDomain):
    tol = 1E-14
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], y1, tol)
    
by1 = BoundaryY1()
by1.mark(boundary_markers, 3)


# For the implimentation of these boundary conditions to be general, we can let the user specify what kind of boundary condition that applies to each of the four boundaries. We set up a Python dictionary for this purpose.


#ENCODE BOUNDARY CONDITIONS (MAY NEED ALTERING OR REMOVAL)
boundary_conditions = {0: {'Robin':     (es, T_a4)},
                       1: {'Robin':     (es, T_a4)},
                       2: {'Robin':     (es, T_a4)},
                       3: {'Robin':     (es, T_a4)}}


ds = Measure('ds', domain=mesh, subdomain_data = boundary_markers)
integrals_R = []
T = TrialFunction(V)
v = TestFunction(V)
T_ = Function(V)
T_n = interpolate(T_am, V)


#Account for Radiation
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        integrals_R.append(es*dt*(T**4 - T_a4)*v*ds(i))

#----------------------------------------

F = T*v*dx + dot(grad(T), grad(v))*dx*dt + dt*(1/kappa)*Q*v*ds(3) + (1/kappa)*sum(integrals_R)-T_n*v*dx

T_ = Function(V)
F = action(F,T_)
J = derivative(F, T_, T)

t = dt
while t <= time:
    print 'time =', t
    problem = NonlinearVariationalProblem(F, T_, J)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()
    t += dt
    T_n.assign(T_)

Any suggestions would be most welcome. Thank you in advance!

2

Hi,

your first error comes from multiplying the form from the right with dt, i.e., you look at

form*dx*dt

which the fenics compiler does not recognize. Instead, you have to multiply with dt from the left, i.e.,

dt*form*dx

or

form*dt*dx

Additionally, you may want to use Constant(dt) instead of dt in the above two formulas.

For your second question, you might have to add

Q.t = t

into the time loop. Further, I would write down the weak form F directly in terms of a Function, i.e.,

T = Function(V)
F = T*v*dx + dot(grad(T), grad(v))*dx*dt + dt*(1/kappa)*Q*v*ds(3) + ...
J = derivative(F, T)

However, I think your approach is equivalent, I’m just not sure since I only used the one I gave you. Alternatively, you can also skip the NonlinearVariational… part and just use

# your preliminary definitions
T = Function(V)
F = T*v*dx + dot(grad(T), grad(v))*dx*dt + dt*(1/kappa)*Q*v*ds(3) + ...

t = dt
while t <= time:
    print 'time =', t
    Q.t = t
    solve(F==0)
    t += dt
    T_n.assign(T)

Which should give you the same results.

I hope this helps you get going.

3

Thank you very much from the help. That has cleared some things up for me.