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[1], 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], 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 = (
    - inner( k * grad(T), grad(s)) * dx 
    - 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[0] = 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[0] = 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[0] = 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[0] = 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[0] = 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[1], 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], 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 = (
    - inner( k * grad(T), grad(s)) * dx 
    - 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
***
***     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.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[1]+x[0]", 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[1], 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], 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[1]", degree=1, T0=Tref))

bcs = []

F = (
    - inner( k * grad(T), grad(s)) * dx 
    - 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