# Robin bcs with local gradient dependence

The code below should solve a heat transfer problem with Robin boundary conditions, but with the heat transfer coefficient being a function of the local heat flux:

``````import numpy as np
from dolfin import *

mesh = UnitSquareMesh(32, 32)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

class Base_face(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0)
base_face = Base_face()
base_face.mark(boundaries, 1)

class Top_face(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 1)
top_face = Top_face()
top_face.mark(boundaries, 2)

# Build function space
V = FunctionSpace(mesh, 'CG', 1)
T = TrialFunction(V)
s = TestFunction(V)
n = FacetNormal(mesh)

k = Constant(1.0)
Tref = Constant(0.0)
f = Constant(1.0)
#HTC = Constant(1.0)
HTC = Expression('20. * pow(abs(dot(n, k*grad(T))), 0.67 )')

bcs = []

F = (
- HTC * (T - Tref) * s * ds(1)
+ f * s * ds(2)
)

a, L = lhs(F), rhs(F)

# Create VTK files for visualization output
vtkfile_T = File('Temperature.pvd')

T = Function(V, name='Temperature')
solve(a == L, T, bcs, solver_parameters={'linear_solver':'mumps'})

# Save solution to file (VTK)
vtkfile_T << T
``````

However because of `HTC = Expression('20. * pow(abs(dot(n, k*grad(T))), 0.67 )')`
I get the error

``````------------------- Start compiler output ------------------------
/tmp/tmp42fb9sup/dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45.cpp: In member function 'virtual void dolfin::dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const':
/tmp/tmp42fb9sup/dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45.cpp:61:41: error: 'n' was not declared in this scope; did you mean 'yn'?
61 |           values = 20. * pow(abs(dot(n, k*grad(T))), 0.67 );
|                                         ^
|                                         yn
/tmp/tmp42fb9sup/dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45.cpp:61:44: error: 'k' was not declared in this scope
61 |           values = 20. * pow(abs(dot(n, k*grad(T))), 0.67 );
|                                            ^
/tmp/tmp42fb9sup/dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45.cpp:61:51: error: 'T' was not declared in this scope
61 |           values = 20. * pow(abs(dot(n, k*grad(T))), 0.67 );
|                                                   ^
/tmp/tmp42fb9sup/dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45.cpp:61:46: error: 'grad' was not declared in this scope
61 |           values = 20. * pow(abs(dot(n, k*grad(T))), 0.67 );
|                                              ^~~~
/tmp/tmp42fb9sup/dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45.cpp:61:37: error: 'dot' was not declared in this scope
61 |           values = 20. * pow(abs(dot(n, k*grad(T))), 0.67 );
|                                     ^~~

-------------------  End compiler output  ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/chbrago/Documents/Fenics/Heat_Spreading/jitfailure-dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45
Traceback (most recent call last):
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/jit/jit.py", line 165, in compile_class
module, signature = dijitso_jit(cpp_data, module_name, params,
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
return local_jit(*args, **kwargs)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/jit/jit.py", line 103, in dijitso_jit
return dijitso.jit(*args, **kwargs)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dijitso/jit.py", line 216, in jit
raise DijitsoError("Dijitso JIT compilation failed, see '%s' for details"
dijitso.jit.DijitsoError: Dijitso JIT compilation failed, see '/home/chbrago/Documents/Fenics/Heat_Spreading/jitfailure-dolfin_expression_e5fd8dcbbb9e14a75c40244a34f14c45' for details

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/home/chbrago/Documents/Fenics/Heat_Spreading/test_Fenics_HTC-q.py", line 31, in <module>
HTC = Expression('20. * pow(abs(dot(n, k*grad(T))), 0.67 )')
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/function/expression.py", line 400, in __init__
self._cpp_object = jit.compile_expression(cpp_code, params)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/function/jit.py", line 158, in compile_expression
expression = compile_class(cpp_data, mpi_comm=mpi_comm)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/jit/jit.py", line 170, in compile_class
raise RuntimeError("Unable to compile C++ code with dijitso")
RuntimeError: Unable to compile C++ code with dijitso
``````

Why are you wrapping this as an expression? This should simply be
`HCT=20. * pow(abs(dot(n, k*grad(T))), 0.67 )`

Without the wrapping I get the error

``````Found Argument in <Power id=139710232517456>, this is an invalid expression.
Traceback (most recent call last):
File "/home/chbrago/Documents/Fenics/Heat_Spreading/test_Fenics_HTC-q.py", line 42, in <module>
a, L = lhs(F), rhs(F)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/formoperators.py", line 79, in lhs
return compute_form_lhs(form)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py", line 379, in compute_form_lhs
return compute_form_with_arity(form, 2)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py", line 345, in compute_form_with_arity
return map_integrands(_transform, form)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/map_integrands.py", line 38, in map_integrands
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/map_integrands.py", line 38, in <listcomp>
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/map_integrands.py", line 46, in map_integrands
return itg.reconstruct(function(itg.integrand()))
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py", line 341, in _transform
e, provides = pe.visit(e)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 110, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/transformer.py", line 114, in visit
r = h(o)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py", line 73, in expr
error("Found Argument in %s, this is an invalid expression." % ufl_err_str(x))
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Found Argument in <Power id=139710232517456>, this is an invalid expression.

``````

Your problem is non-linear due to the `pow(...,0.67)`, and thus you would need to define `T` as a `Function`, not a trial function, and solve the problem as `F = a - L`, `solve(F ==0, ...)`

The following code now runs

``````import numpy as np
from dolfin import *

mesh = UnitSquareMesh(32, 32)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

class Base_face(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0)
base_face = Base_face()
base_face.mark(boundaries, 1)

class Top_face(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 1)
top_face = Top_face()
top_face.mark(boundaries, 2)

# Build function space
V = FunctionSpace(mesh, 'CG', 1)
#T = TrialFunction(V)
T = Function(V, name='Temperature')
s = TestFunction(V)
n = FacetNormal(mesh)

k = Constant(1.0)
Tref = Constant(0.0)
f = Constant(1.0)
#HTC = Constant(1.0)
HTC = Constant(1.0) * pow(abs(dot(n, k*grad(T))), 0.67)

bcs = []

F = (
- HTC * (T - Tref) * s * ds(1)
+ f * s * ds(2)
)

a, L = lhs(F), rhs(F)

# Create VTK files for visualization output
vtkfile_T = File('Temperature.pvd')

F = a - L
solve(F==0, T, bcs)

# Save solution to file (VTK)
vtkfile_T << T
``````

but yields a convergence issue:

``````Form has no parts with arity 2.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.754e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 2: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 3: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 4: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 5: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 6: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 7: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 8: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 9: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 10: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 11: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 12: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 13: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 14: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 15: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 16: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 17: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 18: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 19: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 20: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 21: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 22: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 23: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 24: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 25: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 26: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 27: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 28: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 29: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 30: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 31: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 32: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 33: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 34: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 35: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 36: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 37: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 38: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 39: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 40: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 41: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 42: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 43: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 44: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 45: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 46: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 47: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 48: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Newton iteration 49: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton iteration 50: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
Traceback (most recent call last):
File "/home/chbrago/Documents/Fenics/Heat_Spreading/test_Fenics_HTC-q.py", line 50, in <module>
solve(F==0, T, bcs)#, solver_parameters={'linear_solver':'mumps'})
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/fem/solving.py", line 220, in solve
_solve_varproblem(*args, **kwargs)
File "/home/chbrago/.conda/envs/cadquery/lib/python3.9/site-packages/dolfin/fem/solving.py", line 266, in _solve_varproblem
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
***
***
*** 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.1.0
*** Git changeset:  d0e78c11e813d2decc337bbee614a33128aa539f
*** -------------------------------------------------------------------------
``````

I suspect the problem comes from `Form has no parts with arity 2.` but I don’t know why.

This comes from

which you do not need, as `F` is properly defined at

You can remove

In general, solving a non-linear problem with a non-integer power on the unknown is challenging, see for instance:

You need a good initial guess, which i guess (completely random) at:

``````T.interpolate(Expression("3*x+x", degree=1))
``````

which makes the problem solvable.

It was indeed because of the initial conditions. Thanks to all the contributors for the help, the code is now working, I paste it below for reference:

``````import numpy as np
from dolfin import *

mesh = UnitSquareMesh(32, 32)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

class Base_face(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 0)
base_face = Base_face()
base_face.mark(boundaries, 1)

class Top_face(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x, 1)
top_face = Top_face()
top_face.mark(boundaries, 2)

# Build function space
V = FunctionSpace(mesh, 'CG', 1)
T = Function(V, name='Temperature')
s = TestFunction(V)
n = FacetNormal(mesh)

k = Constant(0.1)
Tref = Constant(20.0)
f = Constant(10.0)
HTC = Constant(1.0) * pow(abs(dot(n, k*grad(T))), 0.67)

T.interpolate(Expression("T0+x", degree=1, T0=Tref))

bcs = []

F = (
- HTC * (T - Tref) * s * ds(1)
+ f * s * ds(2)
)

# Create VTK files for visualization output
vtkfile_T = File('Temperature.pvd')

solve(F==0, T, bcs)

"""
J = derivative(F, T) # Jacobian
rtol = 1e-5 * 1e-4
atol = 1e-2 * 1e-8
problem = NonlinearVariationalProblem(F, T, bcs, J)
solver  = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['relaxation_parameter'] = 1.
prm['newton_solver']['relative_tolerance'] = rtol
prm['newton_solver']['absolute_tolerance'] = atol
prm['newton_solver']['maximum_iterations'] = 20
prm['newton_solver']['error_on_nonconvergence'] = False
prm['newton_solver']['linear_solver'] = 'mumps'
solver.solve()
"""

# Save solution to file (VTK)
vtkfile_T << T
``````