Error when using adaptive linear solver and subdomains

Question

Hi, I am trying to use the adaptive linear solver to adapt the mesh for a heat conduction problem. For the problem I am using a python class that uses subdomains from mshr to have non uniform material properties (thermal conductivities).

When I run the problem with an adaptive tolerance of 1.0e-3 the code runs but when I run it with a lower tolerance like 1.0e-5 and the mesh is being subdivided I get the following error:

Generating forms required for error control, this may take some time...
Solving variational problem adaptively
Adaptive iteration 0
  Solving linear variational problem.
  Value of goal functional is -0.759127.
  Solving linear variational problem.
  Number of cells increased from 1100 to 1187 (7.9% increase).
  Interpolate from parent to child
Adaptive iteration 1
  Solving linear variational problem.
Traceback (most recent call last):
  File "exp.py", line 102, in <module>
    solver.solve(tol)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/adaptivesolving.py", line 81, in solve
    cpp.AdaptiveLinearVariationalSolver.solve(self, tol)
  File "exp.py", line 76, in eval_cell
    if markers[cell.index] == 0:
  File "/usr/lib/python2.7/dist-packages/dolfin/cpp/mesh.py", line 4472, in __getitem__
    raise IndexError("index out of range")
IndexError: index out of range
Aborted (core dumped)

The code is shown below:

from fenics import *
from mshr import * #library for meshing
import math as m   #Import like this to prevent overiding of fenics' UFL math functions

# Use compiler optimizations
parameters["form_compiler"]["cpp_optimize"] = True
# Allow approximating values for points that may be generated outside
# of domain (because of numerical inaccuracies)
parameters["allow_extrapolation"]  = True
# Choose refinement algorithm (needed for solver to choose right refinement algorithm)
parameters["refinement_algorithm"] = "plaza_with_parent_facets"

R_inner  =  0.009525
R_2      =  0.025
R_outer  =  0.06
h_p   =  0.50
h_2   =  0.30
T_g   = -70.0 #degrees C
T_inf =  20.0 #degrees C
#Conduction values
k_1  =  16.3
k_2  =  0.25
#Convection values
h_a    =  10.0  #Free convection
h_g    =  500.0 #Forced convection

#Generate mesh with sub-domains
domain = Rectangle(Point(R_inner,0.0), Point(R_outer,h_p/2.0))
h = h_2/2.0
pipe_vertices = [Point(R_inner,0.0), Point(R_2,0.0), Point(R_2,h), Point(R_inner,h)]
pipe   = Polygon(pipe_vertices)
domain.set_subdomain(1, pipe)
meshResolution = 30
mesh = generate_mesh(domain, meshResolution)

# Define boundaries
tol = 1.0e-2
class insideEdge(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_inner)

class outsideEdge(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], R_outer)

insideEdge  = insideEdge()
outsideEdge = outsideEdge()

boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
insideEdge.mark(boundaries, 1)
outsideEdge.mark(boundaries, 2)
ds = Measure('ds')[boundaries]

# Define function space
V = FunctionSpace(mesh, 'P', 1)

# Define subdomain markers
markers = MeshFunction('size_t', mesh, 2, mesh.domains())

# Define thermal conductivity
class Thermal_Conductivity(Expression):
    def __init__(self, mesh, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if markers[cell.index] == 0:
            values[0] = k_1
        else:
            values[0] = k_2

k = Thermal_Conductivity(mesh, degree=1)

# Define variational problem
h_g = Constant(h_g)
T_g = Constant(T_g)
h_a = Constant(h_a)
T_inf = Constant(T_inf)
T = TrialFunction(V)
q = TestFunction(V)
a = k*inner(grad(T), grad(q))*dx() + h_g*T*q*ds(1) + h_a*T*q*ds(2)
L = h_g*T_g*q*ds(1) + h_a*T_inf*q*ds(2)


# Compute solution (adaptive)
T   = Function(V)
M   = T*dx()
tol = 1.0e-5
problem = LinearVariationalProblem(a,L,T)
solver  = AdaptiveLinearVariationalSolver(problem,M)
solver.solve(tol)
solver.summary()
plot(mesh.root_node(), title='Original coarse mesh')
plot(mesh.leaf_node(), title='Fine mesh')
plot(T.root_node(), title='Temperature distribution (coarse)')
plot(T.leaf_node(), title='Temperature distribution (fine)')

# Plot heat flow
plot(-k*grad(T.leaf_node()), scale=1.5, title='Heat flow vector field')
interactive()

Answer 1

When running the Adaptive linear solver, the markers used in Thermal_Conductivity has to be updated at every solve call. In stead consider doing something like this:

V0 = FunctionSpace(mesh, "DG", 0)            
k = interpolate(Thermal_Conductivity(mesh, degree=1), V0)  

By the way Thermal_Conductivity takes and argument mesh which is not used, but instead it sets self.markers which is defined in the line above this class. I guess what you meant to do was something like the following:

class Thermal_Conductivity(Expression):
    def __init__(self, markers, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 0:
            values[0] = k_1
        else:
            values[0] = k_2

V0 = FunctionSpace(mesh, "DG", 0)            
k = interpolate(Thermal_Conductivity(markers, degree=1), V0)

To make the plot work in the final line you can do e.g

mesh_leaf = mesh.leaf_node()
V_leaf = FunctionSpace(mesh_leaf, "DG", 0)
markers_leaf = MeshFunction("size_t", mesh_leaf, 2, mesh_leaf.domains())
k_leaf = interpolate(Thermal_Conductivity(markers_leaf, degree=1), V_leaf)
plot(-k_leaf*grad(T.leaf_node()), scale=1.5, title='Heat flow vector field')

Comments

Excellent that's what I was looking for. I was actually also trying to plot the material properties and ended up using DG elements for that as well. The code can use different domains now. For the final plot however I get the following error:

Traceback (most recent call last):
  File "insulatedConnector.py", line 174, in <module>
    k_leaf = interpolate(Thermal_Conductivity(markers_leaf, degree=1), V_leaf)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/interpolation.py", line 66, in interpolate
    Pv.interpolate(v)
  File "insulatedConnector.py", line 120, in eval_cell
    if markers[cell.index] == 0:
  File "/usr/lib/python2.7/dist-packages/dolfin/cpp/mesh.py", line 4472, in __getitem__
    raise IndexError("index out of range")
IndexError: index out of range
Aborted (core dumped)

Do you know what this this might be caused by?

---

Did you make the changes I proposed in the definition of Thermal_Conductivity? In eval_cell there should be self.markers and not markers.

---

I didn't change the second markers to self.markers it works now thank you