A problem about Boundary condition setting

General
1

Hi，I am learning to use fenics to solve 3D heat conduction problems.When I was replicating the code，I met some difficulties. The simulation object is a simple model consisting of an air region, a wire region and a substrate. The question is how to set the boundary conditions so that the air dissipates the heat from the substrate and the heated area. When I set the Constant temperature boundary at the top of the air area and the bottom of the substrate,I get the following error

Discontinuous type Coefficient must be restricted.
Traceback (most recent call last):
  File "/mnt/write_air.py", line 156, in <module>
    problem = NonlinearVariationalProblem(weak_form, TsteadyVsteady, bcs, J = Jacsteady)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py", line 157, in __init__
    F = Form(F, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 44, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 50, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 100, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 66, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 143, in compile_form
    prefix, parameters, jit)
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 90, in analyze_ufl_objects
    for form in forms)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 90, in <genexpr>
    for form in forms)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 174, in _analyze_form
    do_apply_restrictions=True)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 321, in compute_form_data
    form = apply_restrictions(form)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py", line 165, in apply_restrictions
    only_integral_type=integral_types)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 47, in map_integrand_dags
    form, only_integral_type)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 28, in map_integrands
    for itg in form.integrals()]
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 28, in <listcomp>
    for itg in form.integrals()]
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 35, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
  File "/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py", line 46, in <lambda>
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 26, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)
  File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 73, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py", line 88, in reference_value
    g = self(f)
  File "/usr/lib/python3/dist-packages/ufl/corealg/multifunction.py", line 89, in __call__
    return self._handlers[o._ufl_typecode_](o, *args)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py", line 134, in coefficient
    return self._require_restriction(o)
  File "/usr/lib/python3/dist-packages/ufl/algorithms/apply_restrictions.py", line 47, in _require_restriction
    error("Discontinuous type %s must be restricted." % o._ufl_class_.__name__)
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Discontinuous type Coefficient must be restricted.

The code is shown below

import numpy as np
import meshio
from scipy import interpolate as sci_interpolate

msh = meshio.read("Sample.msh")

print("Available cell types in msh.cells:")
for cell in msh.cells:
    print(f"Cell type: {cell.type}, Number of cells: {len(cell.data)}")

tetra_cells = [cells.data for cells in msh.cells if cells.type == "tetra"]
triangle_cells = [cells.data for cells in msh.cells if cells.type == "triangle"]

if tetra_cells:
    cells = np.vstack(np.array(tetra_cells))
    tetra_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
    tetra_mesh = meshio.Mesh(points=msh.points, cells=[("tetra", cells)], cell_data={"name_to_read": [tetra_data]})
    meshio.xdmf.write("mesh_test0.xdmf", tetra_mesh)
else:
    raise ValueError("No tetra cells found in the mesh!")

if triangle_cells:
    facet_cells = np.vstack(np.array(triangle_cells))
    facet_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
    facet_mesh = meshio.Mesh(points=msh.points, cells=[("triangle", facet_cells)], cell_data={"name_to_read": [facet_data]})
    meshio.xdmf.write("facet_mesh_test0.xdmf", facet_mesh)
else:
    raise ValueError("No triangle cells found in the mesh!")

# Use dolfin below.
from dolfin import *
mesh = Mesh()

# Read mesh.
with XDMFFile("mesh_test0.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("facet_mesh_test0.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh_test0.xdmf") as infile:
    infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = Measure("dx", domain=mesh,subdomain_data=cf)
# External surface
ds = Measure("ds", domain=mesh, subdomain_data=mf)
# Internal surface
dS = Measure("dS", domain=mesh, subdomain_data=mf)
subdomains = cf
boundary_parts = mf

# Volumes
# 1 thin film
# 2 heater wire
# 3 air

# Surfaces
# 4 left end heater wire.
# 5 right end heater wire.
# 6 heater_front
# 7 "heater_behind"
# 8 "heater_top"
# 9 Interface between wire and thin film.
# 10 "thin_film_top"
# 11 Thin film bottom
# 12 "air_top"
# 13 "air_left"
# 14 "air_front"
# 15 "air_behind"
# 16 "air_right"

dx_air=dx(subdomain_id=3)
dx_thin_film = dx(subdomain_id=1)
dx_wire = dx(subdomain_id=2)
ds_bottom_thin_film = ds(subdomain_id=11)
dS_interface=dS(subdomain_id=9)

P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

ME = FunctionSpace(mesh, P1*P1)

h_air_film = Constant(5.0)  #Heat exchange coefficient W/(m²·K)
h_air_wire = Constant(5.0)  #Heat exchange coefficient W/(m²·K)
T_ref = 300.0
kappa_air = Constant(0.025)  # Thermal conductivity W/(m·K)
Dens_air = Constant(1.225)  # density kg/m^3
cp_air = Constant(1005)  #  Heat Capacity J/(kg·K)
kappa_wire = Constant(75.0)
kappa_thin_film = Constant(7.8)
rho_interpolated = sci_interpolate.InterpolatedUnivariateSpline([10,400], [2.25e-9, 1.14e-7], k=1) # linear interpolation
current = Constant(90e-1)
J_0 = current / assemble(1 * dS(4))

dTsteadyVsteady = TrialFunction(ME)
usteady, vsteady = TestFunctions(ME)
TsteadyVsteady = Function(ME)
Tsteady, Vsteady = split(TsteadyVsteady)

V = FunctionSpace(mesh, 'CG', 1)
rho_steady = project(float(rho_interpolated(T_ref)),V)
rho_thin_film = project(5.0e-3,V)

TsteadyVsteady_init=Expression((T_ref,'0.225'),degree=1)
TsteadyVsteady.interpolate(TsteadyVsteady_init)
J_injected = J_0
J_local = -(1.0/rho_steady)*grad(Vsteady)
J_thin_film = -(1.0/rho_thin_film)*grad(Vsteady)

convection_term_film = h_air_film * (Tsteady - T_ref) * usteady * dS(subdomain_id=10)
convection_term_wire = h_air_wire * (Tsteady - T_ref) * usteady * dS(subdomain_id=8)


Fourier_terms = -dot(kappa_wire * grad(Tsteady), grad(usteady)) * dx_wire - dot(kappa_thin_film * grad(Tsteady), grad(usteady)) * dx_thin_film  - dot(kappa_air * grad(Tsteady), grad(usteady)) * dx_air
Joule_term = dot(rho_steady*J_local, J_local) * usteady*dx_wire + dot(rho_thin_film*J_thin_film, J_thin_film) * usteady*dx_thin_film# Add Joule heat in thin film.

FsteadyT = Fourier_terms + Joule_term  + convection_term_wire + convection_term_film
FsteadyV = dot(1.0/rho_steady * grad(Vsteady), grad(vsteady))*dx_wire + J_injected * vsteady * ds(4) + dot(1.0/rho_thin_film * grad(Vsteady), grad(vsteady))*dx_thin_film

weak_form = FsteadyT + FsteadyV
Jacsteady = derivative(weak_form,TsteadyVsteady,dTsteadyVsteady)

V_left = Constant(0.0)
bc0 = DirichletBC(ME.sub(0), T_ref, boundary_parts, 11)
bc1 = DirichletBC(ME.sub(1), V_left, boundary_parts, 4)
bc_air = DirichletBC(ME.sub(0), T_ref, boundary_parts, 12)
#bc_air_2 = DirichletBC(ME.sub(0), T_ref, boundary_parts, 13)
#bc_air_3 = DirichletBC(ME.sub(0), T_ref, boundary_parts, 14)
#bc_air_4 = DirichletBC(ME.sub(0), T_ref, boundary_parts, 15)
#bc_air_5 = DirichletBC(ME.sub(0), T_ref, boundary_parts, 16)
bcs = [bc0, bc1, bc_air]

problem = NonlinearVariationalProblem(weak_form, TsteadyVsteady, bcs, J = Jacsteady)
solver = NonlinearVariationalSolver(problem)
solver.solve()

Thank you very much.

2

Searching the forum / google for the error message, you’ll find that it relates to specifying which side of an interface should be evaluated, if you want to integrate a double-valued function.

In your case, I assume it goes wrong in the lines

convection_term_film = h_air_film * (Tsteady - T_ref) * usteady * dS(subdomain_id=10)
convection_term_wire = h_air_wire * (Tsteady - T_ref) * usteady * dS(subdomain_id=8)

which probably should look something like:

convection_term_film = h_air_film * (Tsteady("+") - T_ref) * usteady("+") * dS(subdomain_id=10)
convection_term_wire = h_air_wire * (Tsteady("-") - T_ref) * usteady("-") * dS(subdomain_id=8)

I’m sure that’s not entirely correct yet, but it may give you a handle on figuring out exactly how to fix your formulation.

1 Like
3

Thank you for your reply！
I change the code and it works, but there is a problem that the computation does not converge. I think the boundary conditions are still wrong. The error message is as follows：

Traceback (most recent call last):
  File "/mnt/write_air.py", line 157, in <module>
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------