# Unable to compute the curl of a vector field

Hi, I work on a 2D problem and my use-case code is very long, I know I must create a MWE but I greatly struggle, as I am not used to employ a mesh that doesn’t come from Gmsh. Currently my attempt is unfinished and I will post it at the end (so maybe we can finish the MWE and I might use it for my next questions).

My problem lies in errors when computing the curl of a vector field. The relevant part of my code causing me problem is:

(...)
deg = 4
# Define ME function space
el = FiniteElement("CG", mesh.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)
(...)
S_xx = 3e-3
S_yy = S_xx * 10
rho = 1
sigma = 1.0 / rho
κ = 1.8
Seebeck_tensor = as_tensor([[S_xx, 0], [0, S_yy]])
(...)
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J = Jac)
(...)
n, converged = solver.solve(TempVolt)

# Current density.
J = VectorFunctionSpace(mesh, ("DG", deg-1))
Je = Function(J)
Je_flux_calculator = Expression(-sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp), J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)

# Curl of current density.
Jcurl = functionspace(mesh, ("DG", deg-2, (mesh.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(as_vector([J_vector[0],J_vector[1],0])), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)


with traceback:

Traceback (most recent call last):
File "/root/shared/TE_simulation/Bridgman_Nye_example.py", line 434, in <module>
post_proc.post_proc()
File "/root/shared/TE_simulation/Bridgman_Nye_example.py", line 363, in post_proc
J_curl = build_vector_field(curl(as_vector([J_vector[0],J_vector[1],0])), mesh, deg)
File "/root/shared/TE_simulation/Bridgman_Nye_example.py", line 51, in build_vector_field
flux_calculator = Expression(expression, Q.element.interpolation_points())
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 149, in __init__
self._ufcx_expression, module, self._code = jit.ffcx_jit(comm, (e, _X),
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 210, in ffcx_jit
r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 247, in compile_expressions
raise e
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 238, in compile_expressions
impl = _compile_objects(decl, expressions, expr_names, module_name, p, cache_dir,
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 95, in analyze_ufl_objects
processed_expression = _analyze_expression(original_expression, options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 117, in _analyze_expression
expression = ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(expression)
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 148, in apply_algebra_lowering
return map_integrand_dags(LowerCompoundAlgebra(), expr)
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in map_integrand_dags
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 66, in map_integrands
return function(integrand)
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in <lambda>
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 34, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress,
File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 100, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 142, in curl
return as_vector((c(1, 2), c(2, 0), c(0, 1)))
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 135, in c
return a[j].dx(i) - a[i].dx(j)
File "/usr/local/lib/python3.10/dist-packages/ufl/exproperators.py", line 435, in _dx
File "/usr/local/lib/python3.10/dist-packages/ufl/differentiation.py", line 252, in __new__
dim = find_geometric_dimension(f)
File "/usr/local/lib/python3.10/dist-packages/ufl/domain.py", line 396, in find_geometric_dimension
raise ValueError("Cannot determine geometric dimension from expression.")
ValueError: Cannot determine geometric dimension from expression.


I also tried the Expression: J_vector[1].dx(0)-J_vector[0].dx(1) and as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)]), without any success.

I am sorry for not posting my whole code (I can, but it’s hard to parse, I’d rather make a MWE).

And here’s my attempt to the MWE:

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector)

from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8
T_cold = 100

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

def boundary_D_left(x):
return np.isclose(x[0], 0)

def boundary_D_right(x):
return np.isclose(x[0], 1)

# Define ME function space
deg = 3
el = FiniteElement("CG", mesh.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

# The Seebeck coefficient has to be an anisotropic tensor for the Bridgman heat to take place.
s_xx, s_yy = seebeck_components
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

# Define the boundary conditions
#left_facets = facet_markers.find(inj_current_curve)
#right_facets = facet_markers.find(out_current_curve)
#left_dofs = locate_dofs_topological(ME.sub(1), mesh.topology.dim-1, left_facets)

#left_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, left_facets)
#right_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, right_facets)
left_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_left)
right_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_right)

bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(T_cold + ΔT), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right]

x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))

# Weak form.
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve)

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
Current density injected: {J_inj}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-13
solver.atol = 1e-13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(mesh)

# Current density.
J = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)

Jcurl = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)

J_curl = build_vector_field(as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)]), mesh, deg)
J_curl = build_vector_field(J_vector[1].dx(0)-J_vector[0].dx(1), mesh, deg)
J_curl = build_vector_field(curl(as_vector([J_vector[0],J_vector[1],0])), mesh, deg)


But AttributeError: module 'dolfinx.mesh' has no attribute 'ufl_cell' so I’ll have to debug myself.

Here is my new attempt to make a working MWE:

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector)

from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8
T_cold = 100

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

def boundary_D_left(x):
return np.isclose(x[0], 0)

def boundary_D_right(x):
return np.isclose(x[0], 1)

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

left_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_left)
right_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_right)

bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(T_cold + ΔT), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right]

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=facet_markers)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))

# Weak form.
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve)

weak_form = F_T + F_V

Jac = derivative(weak_form,TempVolt,dTV)

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-13
solver.atol = 1e-13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)

# Current density.
J = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)

Jcurl = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)


The traceback:

Traceback (most recent call last):
File "/root/shared/debug_discourse/curl_MWE.py", line 56, in <module>
left_dofs_temp = locate_dofs_geometrical(ME.sub(0), boundary_D_left)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/bcs.py", line 51, in locate_dofs_geometrical
return _cpp.fem.locate_dofs_geometrical(V._cpp_object, marker)  # type: ignore
RuntimeError: Cannot tabulate coordinates for a FunctionSpace that is a subspace.



Try to follow one of the demos that uses mixed_element, for instance dolfinx/python/demo/demo_stokes.py at v0.8.0 · FEniCS/dolfinx · GitHub . (the fact that locate_dofs_topological is used there instead of locate_dofs_geometrical shouldn’t make any difference)

1 Like

Grazie Francesco, I am able to make some progress thanks to you.

I am a bit stuck. My latest attempt:

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix,
)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8
T_cold = 100

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
def boundary_D_left(x):
return np.isclose(x[0], 0)

def boundary_D_right(x):
return np.isclose(x[0], 1)

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
facets = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T_cold = Function(Q)
facets = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T0, Q), 1, facets)
bc_temp_right = dirichletbc(T_cold, dofs, T0)

bcs = [bc_temp_left, bc_temp_right]

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain)#, subdomain_data=facet_markers)

the_current = 1.0
J_inj = the_current / assemble_scalar(form(1 * ds))
J_out = the_current / assemble_scalar(form(1 * ds))


1st question: I do not understand in the example where is the value of “noslip” specified. I had T_cold and T_hot as floats before (or Constants), but now they are “Function(Q)” without any value. How do I specify their values? I hope not by interpolating a constant into Q?

2nd question: How do I specify ds(1) and ds(2) for example for the 2 surfaces or curves where I apply the Dirichlet b.c.? I would need to mark the curves? This is not done in the example? I am confused on this.

A dolfinx.fem.Function has zero values by default, which is indeed the value for no split boundary conditions. Interpolating the constant value into a CG space is simply setting all entries of the dolfinx.fem.Function.x.array to that value. If someone else reads your code, you may want to be more expressive and interpolate.

Since you are already locating the entities, you may want to use a meshtags defined on those entities as in dolfinx/python/test/unit/fem/test_assemble_domains.py at v0.8.0 · FEniCS/dolfinx · GitHub
In principle you could even avoid the meshtags object if adopting the strategy in dolfinx/python/test/unit/fem/test_assemble_domains.py at v0.8.0 · FEniCS/dolfinx · GitHub . I’ll leave it to you to decide what you prefer.

1 Like

Thanks once more Francesco.
I could proceed further with my MWE, but my code fails to solve the problem (Newton method fails to converge, and in fact diverges at the very first step).
I tried to simplify the equation (made S tensor equal to 0, removed Fourier and Bridgman terms in the weak form), instead of -nan I sometimes got “inf” according to the logger.
The traceback is:

area of material: 1.0
1.0000000000000002
1.0000000000000002
------- Pre-processing --------
Current density injected: 1.0000000000000002
Length of the side where current is injected: 0.9999999999999999
Length of the side where current leaves the wire: 0.9999999999999999
This should correspond to the current injected: 0.9999999999999999

2024-06-22 21:22:24.037 (   0.740s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 0: r (abs) = -nan (tol = 1e-06) r (rel) = -nan(tol = 1e-06)
2024-06-22 21:22:24.039 (   0.743s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.045 (   0.749s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.005737, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.046 (   0.749s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 1: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.048 (   0.751s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.051 (   0.754s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.003069, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.052 (   0.755s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.054 (   0.757s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.054 (   0.757s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000033, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.055 (   0.758s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 3: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.057 (   0.760s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.057 (   0.760s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000047, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.058 (   0.761s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 4: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.060 (   0.763s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.060 (   0.763s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000021, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.060 (   0.764s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 5: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.062 (   0.766s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.062 (   0.766s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.063 (   0.766s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 6: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.065 (   0.768s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.065 (   0.768s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.066 (   0.769s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 7: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.068 (   0.771s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.068 (   0.771s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000021, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.068 (   0.772s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 8: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.070 (   0.774s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.070 (   0.774s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.071 (   0.774s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 9: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.073 (   0.777s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.073 (   0.777s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000029, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.074 (   0.777s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 10: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.076 (   0.779s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.076 (   0.780s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000038, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.077 (   0.780s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 11: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.079 (   0.782s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.079 (   0.782s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.079 (   0.783s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 12: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.081 (   0.785s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.081 (   0.785s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000017, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.082 (   0.785s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 13: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.084 (   0.787s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.084 (   0.788s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000018, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.085 (   0.788s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 14: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.087 (   0.790s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.087 (   0.790s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000017, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.087 (   0.791s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 15: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.089 (   0.793s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.089 (   0.793s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000039, 0.000000, 0.000000 (PETSc Krylov solver)
(...)
2024-06-22 21:22:24.317 (   1.020s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000021, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.317 (   1.021s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 99: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
2024-06-22 21:22:24.319 (   1.023s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2024-06-22 21:22:24.319 (   1.023s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.000022, 0.000000, 0.000000 (PETSc Krylov solver)
2024-06-22 21:22:24.320 (   1.023s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 100: r (abs) = nan (tol = 1e-06) r (rel) = nan(tol = 1e-06)
Traceback (most recent call last):
File "/root/shared/debug_discourse/curl_MWE.py", line 153, in <module>
n, converged = solver.solve(TempVolt)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py", line 46, in solve
n, converged = super().solve(u.vector)
RuntimeError: Newton solver did not converge because maximum number of iterations reached
2024-06-22 21:22:24.366 (   1.069s) [main            ]             loguru.cpp:526   INFO| atexit


My full code:

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix,
)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
return np.isclose(x[0], 0)

def boundary_D_right(x):
return np.isclose(x[0], 1)

def T_cold_func(x):
return eval(str(np.asarray([100.0])))
def T_hot_func(x):
return eval(str(np.asarray([200.0])))

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 0,0#1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T_cold = Function(Q)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T0, Q), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T0)

bcs = [bc_temp_left, bc_temp_right]

# Mesh tags to apply Neumann b.c.
tdim = domain.topology.dim
domain.topology.create_entities(tdim - 1)
facet_map = domain.topology.index_map(tdim - 1)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

# T_cold surface markers.
inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, tdim - 1, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

the_current = 1.0
J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

print(f'area of material: {assemble_scalar(form(1 * dx))}')
print(J_inj)
print(J_out)
# Weak form.
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve)

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
Current density injected: {J_inj}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))}
This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-6#13
solver.atol = 1e-6#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)

# Current density.
J = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)

Jcurl = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)

J_curl = build_vector_field(as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)]), domain, deg)
J_curl = build_vector_field(J_vector[1].dx(0)-J_vector[0].dx(1), domain, deg)
J_curl = build_vector_field(curl(as_vector([J_vector[0],J_vector[1],0])), domain, deg)


I’m very close to have it working, which would be awesome… But not sure how to make it converge now.

An initial residual with value nan or inf typically means that the initial guess you provided to the Newton solver, i.e. the value of TempVolt, is incorrect. Unless you interpolate some expression in TempVolt is zero, which may lead to inf/nan as soon as your residual has expressions like 1/TempVolt (or similar algebraic operations that do not make sense when the solution is zero).

1 Like

I see…
I’m having a hard time to do that, i.e. interpolate an intial guess for temp and volt.
Here is my attempt (exact same whole code as above, with this addition):

def initial_guess_temp():
return np.asarray([100.0])

def initial_guess_volt():
return np.asarray([0.0])
TempVolt.sub(0).interpolate(initial_guess_temp, facet_indices)
TempVolt.sub(1).interpolate(initial_guess_volt, facet_indices)
# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)


This yields:

area of material: 1.0
0.0
0.0
------- Pre-processing --------
Current density injected: 0.0
Length of the side where current is injected: 0.9999999999999999
Length of the side where current leaves the wire: 0.9999999999999999
This should correspond to the current injected: 0.0

Traceback (most recent call last):
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 397, in interpolate
_interpolate(u, cells)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 373, in _interpolate
self._cpp_object.interpolate(u, cells, nmm_interpolation_data)
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. (self: dolfinx.cpp.fem.Function_float64, f: numpy.ndarray[numpy.float64], cells: numpy.ndarray[numpy.int32]) -> None
2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32], nmm_interpolation_data: Tuple[List[int], List[int], List[float], List[int]]) -> None
3. (self: dolfinx.cpp.fem.Function_float64, expr: dolfinx::fem::Expression<double, double>, cells: numpy.ndarray[numpy.int32]) -> None

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7da93e7f9770>, <function initial_guess_temp at 0x7da93e4ebf40>, array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
143]), ((), (), (), ())

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/root/shared/debug_discourse/curl_MWE.py", line 133, in <module>
TempVolt.sub(0).interpolate(initial_guess_temp, facet_indices)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 401, in interpolate
x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh.geometry, cells)
TypeError: interpolation_coords(): incompatible function arguments. The following argument types are supported:
1. (element: dolfinx.cpp.fem.FiniteElement_float32, V: dolfinx.cpp.mesh.Geometry_float32, cells: numpy.ndarray[numpy.int32]) -> numpy.ndarray[numpy.float32]
2. (element: dolfinx.cpp.fem.FiniteElement_float64, V: dolfinx.cpp.mesh.Geometry_float64, cells: numpy.ndarray[numpy.int32]) -> numpy.ndarray[numpy.float64]

Invoked with: <dolfinx.cpp.fem.FiniteElement_float64 object at 0x7da93e2ed3b0>, <dolfinx.cpp.mesh.Geometry_float64 object at 0x7da93e2eea70>, array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
143])


I also have a hard time to find ME examples with interpolation.

You cannot interpolate on facets only.

What i would do is using compute_incident_entities to find all cells connected to facet_indices and send this into interpolate.

I see, thanks a lot for your help.
What is wrong with my new attempt though?

def initial_guess_temp():
return np.asarray([100.0])

def initial_guess_volt():
return np.asarray([0.0])

the_cells = compute_incident_entities(domain.topology, facet_indices, tdim, tdim-1)

TempVolt.sub(0).interpolate(initial_guess_temp, the_cells)
TempVolt.sub(1).interpolate(initial_guess_volt, the_cells)


The output is:

area of material: 1.0
0.0
0.0
------- Pre-processing --------
Current density injected: 0.0
Length of the side where current is injected: 0.9999999999999999
Length of the side where current leaves the wire: 0.9999999999999999
This should correspond to the current injected: 0.0

Traceback (most recent call last):
File "/root/shared/debug_discourse/curl_MWE.py", line 134, in <module>
the_cells = compute_incident_entities(domain.topology, facet_indices, tdim, tdim-1)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/mesh.py", line 182, in compute_incident_entities
return _cpp.mesh.compute_incident_entities(topology, entities, d0, d1)
TypeError: compute_incident_entities(): incompatible function arguments. The following argument types are supported:
1. (mesh: dolfinx.cpp.mesh.Topology, entities: numpy.ndarray[numpy.int32], d0: int, d1: int) -> numpy.ndarray[numpy.int32]

Invoked with: <dolfinx.cpp.mesh.Topology object at 0x7c21e66ac220>, array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
143]), 2, 1



Please note that it is hard to give you guidance when you provide small code snippets that are not reproducible.
Say, for the problem above, you should start by testing the method on a unit square.

from mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
v = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u = dolfinx.fem.Function(v)

def top_facets(x):
return np.isclose(x[1], 1.0)

tdim = mesh.topology.dim
fdim = tdim -1

mesh.topology.create_connectivity(fdim, tdim)
facets = dolfinx.mesh.locate_entities(mesh, fdim, top_facets)

cells = dolfinx.mesh.compute_incident_entities(mesh.topology, facets, fdim, tdim)

u.interpolate(lambda x: np.sin(x[0]), cells=cells)
u.x.scatter_forward()

with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u], engine="BP4") as bp:
bp.write(0.0)


returning

I’m sorry dokken, even with your code, I am having a hard time. I do not want to do exactly what you’ve done, I believe. I want to set “temp” and “volt” to constant values over the whole mesh, not on particular submeshes.
I am not sure, maybe I am missing the domain.topology.create_connectivity(tdim-1, tdim) part (I tried to add it but this didn’t fix my error), but I am not sure because I do not see it in the code linked by Francesco. Also, you seem to have altered the order of tdim and fdim in the arguments of compute_incident_entities. Which order would be correct?

Here is my full code so far:

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix,
)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags, compute_incident_entities
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
return np.isclose(x[0], 0)

def boundary_D_right(x):
return np.isclose(x[0], 1)

def T_cold_func(x):
return eval(str(np.asarray([100.0])))
def T_hot_func(x):
return eval(str(np.asarray([100.0])))

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 0,0
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])
the_current = 0.0

T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T_cold = Function(Q)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T0, Q), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T0)

bcs = [bc_temp_left, bc_temp_right]

# Mesh tags to apply Neumann b.c.
tdim = domain.topology.dim
domain.topology.create_entities(tdim - 1)
facet_map = domain.topology.index_map(tdim - 1)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

# T_cold surface markers.
inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, tdim - 1, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

print(f'area of material: {assemble_scalar(form(1 * dx))}')
print(J_inj)
print(J_out)
# Weak form.
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve)

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
Current density injected: {J_inj}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))}
This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

#print(facet_indices)
# Create initial guess.
def initial_guess_temp():
return np.asarray([100.0])

def initial_guess_volt():
return np.asarray([0.0])

tdim = mesh.topology.dim
fdim = tdim -1

domain.topology.create_connectivity(tdim-1, tdim)
the_cells = compute_incident_entities(domain.topology, facet_indices, tdim-1, tdim)

TempVolt.sub(0).interpolate(initial_guess_temp, the_cells)
TempVolt.sub(1).interpolate(initial_guess_volt, the_cells)
# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"#"incremental"#"residual"
solver.rtol = 1e-6#13
solver.atol = 1e-6#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)


But why dont you the just do interpolation over the whole mesh.

I suggested compute_incident_entities because you had

in your code, which is simply not correct.

What is wrong with doing
TempVolt.sub(0).interpolate(initial_guess_temp)?

The order I have in my code is correct. First dimension is the dimension of the input entities, second dimension is output dimension.
See for instance: What is the difference between compute_incident_entities and entities_to_geometry?

I guess that’s what I tried and what I thought I was doing. I did not realize I wasn’t doing this.

Answer: TypeError: initial_guess_temp() takes 0 positional arguments but 1 was given. This was the very first thing I tried.
The full traceback is:

Traceback (most recent call last):
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 397, in interpolate
_interpolate(u, cells)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 373, in _interpolate
self._cpp_object.interpolate(u, cells, nmm_interpolation_data)
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
1. (self: dolfinx.cpp.fem.Function_float64, f: numpy.ndarray[numpy.float64], cells: numpy.ndarray[numpy.int32]) -> None
2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32], nmm_interpolation_data: Tuple[List[int], List[int], List[float], List[int]]) -> None
3. (self: dolfinx.cpp.fem.Function_float64, expr: dolfinx::fem::Expression<double, double>, cells: numpy.ndarray[numpy.int32]) -> None

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7b0a0c6ce770>, <function initial_guess_temp at 0x7b0a0c21fa30>, array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63], dtype=int32), ((), (), (), ())

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/root/shared/debug_discourse/curl_MWE.py", line 127, in <module>
TempVolt.sub(0).interpolate(initial_guess_temp)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 402, in interpolate
self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
TypeError: initial_guess_temp() takes 0 positional arguments but 1 was given


Great. Note that this differs from your Github issue (I’m sure you’re aware… but just in case): [BUG]: Packing of codim-0 form not working for mixture of measures · Issue #3256 · FEniCS/dolfinx · GitHub.

Thats because there i am mapping cells to facets, while you are mapping facets to cells.

Change them to

def initial_guess_temp(x):
return np.full(x.shape[1], 100.)

def initial_guess_volt(x):
return np.zeros(x.shape[1])


The interpolate function takes in an x that is (3, num_points) and returns an array of length num_points mapping each entry x[:,i] to a corresponding output value

1 Like

Thanks a lot dokken. I had found an older code of mine and the syntax was a bit different, I couldn’t make it work, but your code did the trick.

I still get the same problem of instant nan values with Newton method. I am guessing something else is wrong in my code.

I might try to impose a Dirichlet b.c. for “volt” although I do not have it in my original (non MWE) code.

I do have a question regarding the spaces of the boundary condition functions:
Should I create different spaces “T0” and “T1” as well as “Q” and “Q1” like so, one for each of the Dirichlet b.c. for the “temp” variable, like so?:

# Define the boundary conditions
T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T1 = ME.sub(0)
Q1, _ = T1.collapse()
T_cold = Function(Q1)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T1, Q1), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T1)


Or is the code I posted previously OK regarding this part?

Either way, this does not fix Newton method failing to converge.

Here is my latest code with the implementation of a Dirichlet b.c. for volt, Newton method fails to converge:

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix,
)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags, compute_incident_entities
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
return np.isclose(x[0], 0)

def boundary_D_right(x):
return np.isclose(x[0], 1)

def T_cold_func(x):
return eval(str(np.asarray([100.0])))
def T_hot_func(x):
return eval(str(np.asarray([100.0])))

def V_left_func(x):
return eval(str(np.asarray([0.0])))

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 0,0#1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])
the_current = 0.0

# Define the boundary conditions
T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_right)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T1 = ME.sub(0)
Q1, _ = T1.collapse()
T_cold = Function(Q1)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T1, Q1), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T1)

V0 = ME.sub(1)
JV, _ = V0.collapse()
V_hot = Function(JV)
V_hot.interpolate(V_left_func)
bc_volt_left = dirichletbc(V_hot, dofs, V0)

bcs = [bc_temp_left, bc_temp_right, bc_volt_left]

# Mesh tags to apply Neumann b.c.
tdim = domain.topology.dim
fdim = tdim -1
domain.topology.create_entities(fdim)
facet_map = domain.topology.index_map(fdim)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

# T_cold surface markers.
inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, fdim, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

print(f'area of material: {assemble_scalar(form(1 * dx))}')
print(J_inj)
print(J_out)
# Weak form.
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve)

weak_form = F_T + F_V

# Solve the PDE.
Jac = derivative(weak_form,TempVolt,dTV)

print(f''' ------- Pre-processing --------
Current density injected: {J_inj}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=domain)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=domain)))}
This should correspond to the current injected: {assemble_scalar(form(J_out*ds(out_current_curve)))}
''')

def initial_guess_temp(x):
return np.full(x.shape[1], 100.0)

def initial_guess_volt(x):
return np.zeros(x.shape[1])

TempVolt.sub(0).interpolate(initial_guess_temp)
TempVolt.sub(1).interpolate(initial_guess_volt)

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"#"incremental"#"residual"
solver.rtol = 1e-6#13
solver.atol = 1e-6#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(TempVolt)
assert(converged)


Consider the following:

import dolfinx.fem
vec = dolfinx.fem.assemble_vector(dolfinx.fem.form(weak_form))
dolfinx.fem.apply_lifting(vec.array, [dolfinx.fem.form(Jac)], [bcs])


This would be part of your RHS of the newton problem (with some extra stuff).
Inspecting vec.array yields:

array([-7.54294650e+307, -2.08755243e+307, -2.02282709e+307, ...,
0.00000000e+000,  0.00000000e+000,  0.00000000e+000])


indicating that something is seriously wrong with either your jacobian or your bcs.
Further investigating, you can compute:

A = dolfinx.fem.assemble_matrix(dolfinx.fem.form(Jac))
print(A.to_scipy().min(), A.to_scipy().max())


and see that the values are bounded.

V_hot.x.array contains nan’s and infinite values, meaning that it is wrong.
I’ve now given you a full recipe on how I address convergence issues, which you should take into account for further questions.

I have no clue why you are creating interpolation functions like this

as I have explicitly stated how to create them in my last post on this thread.

1 Like

Thank you very much dokken, I was sloppy, and unaware of these techniques to debug. This will help me in the future.
I am very happy in that I reach the same point as with my original code (not the MWE).
The code solves for temp and volt and computes “Je”, the current density of electrodynamics. This vector field, at first glance, looks correct to me (it has vortices in this case).
I wish to compute its curl, and this is where I get stuck. Namely, the very last line of my script yields an error.
Here’s the full code, the traceback and a picture of the current density vector field whose curl I am looking for. Note that I removed the Dirichlet b.c. on volt, the code is able to converge with pure Neumann b.c. and no null space specified.

import numpy as np
from dolfinx import log, fem, mesh
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix,
assemble_vector, apply_lifting)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector)
from dolfinx.mesh import locate_entities_boundary, meshtags, compute_incident_entities
from dolfinx.io import VTXWriter

rho = 1
sigma = 1.0 / rho
κ = 1.8

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
x = SpatialCoordinate(domain)

def boundary_D_left(x):
return np.isclose(x[0], 0)

def boundary_D_bottom(x):
return np.isclose(x[1], 0)

def T_cold_func(x):
return np.full(x.shape[1], 100.0)
def T_hot_func(x):
return np.full(x.shape[1], 200.0)

# Define ME function space
deg = 3
el = FiniteElement("CG", domain.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(domain, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

s_xx, s_yy = 1e-3, 1e-4
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])
the_current = 0.0

# Define the boundary conditions
T0 = ME.sub(0)
Q, _ = T0.collapse()
T_hot = Function(Q)
T_hot.interpolate(T_hot_func)
facets_hot = locate_entities_boundary(domain, 1, boundary_D_bottom)
dofs = locate_dofs_topological((T0, Q), 1, facets_hot)
bc_temp_left = dirichletbc(T_hot, dofs, T0)

T1 = ME.sub(0)
Q1, _ = T1.collapse()
T_cold = Function(Q1)
T_cold.interpolate(T_cold_func)
facets_cold = locate_entities_boundary(domain, 1, boundary_D_left)
dofs = locate_dofs_topological((T1, Q1), 1, facets_cold)
bc_temp_right = dirichletbc(T_cold, dofs, T1)
bcs = [bc_temp_left, bc_temp_right]#, bc_volt_left]

tdim = domain.topology.dim
fdim = tdim -1
domain.topology.create_entities(fdim)
facet_map = domain.topology.index_map(fdim)
num_facets = facet_map.size_local + facet_map.num_ghosts
facet_indices = np.arange(0, num_facets)
facet_values = np.zeros_like(facet_indices, dtype=np.intc)

inj_current_curve = 12
out_current_curve = 13
facet_values[facets_cold] = inj_current_curve
facet_values[facets_hot] = out_current_curve
mt_facets = meshtags(domain, fdim, facet_indices, facet_values)

x = SpatialCoordinate(domain)
dx = Measure("dx", domain=domain)#, subdomain_data=cell_markers)
ds = Measure("ds", domain=domain, subdomain_data=mt_facets)

J_inj = the_current / assemble_scalar(form(1 * ds(inj_current_curve)))
J_out = the_current / assemble_scalar(form(1 * ds(out_current_curve)))

Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve)

weak_form = F_T + F_V

Jac = derivative(weak_form,TempVolt,dTV)

def initial_guess_temp(x):
return np.full(x.shape[1], 100.0)

def initial_guess_volt(x):
return np.zeros(x.shape[1])

TempVolt.sub(0).interpolate(initial_guess_temp)
TempVolt.sub(1).interpolate(initial_guess_volt)

# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, J = Jac, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"#"incremental"#"residual"
solver.rtol = 1e-8#13
solver.atol = 1e-8#13
solver.report = True
solver.max_it = 100

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(TempVolt)
assert(converged)
n = FacetNormal(domain)

# Current density.
J = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)

with VTXWriter(MPI.COMM_WORLD, "results/Je_debug.bp", Je, engine="BP4") as vtx:
vtx.write(0.0)

Jcurl = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
Jecurl = Function(Jcurl)
Jecurl_flux_calculator = Expression(curl(J_vector), Jcurl.element.interpolation_points())
Jecurl.interpolate(Jecurl_flux_calculator)

# Also tried:
#Jecurl_flux_calculator = Expression(as_vector([0,0,J_vector[1].dx(0)-J_vector[0].dx(1)]), Jcurl.element.interpolation_points())
#Jecurl_flux_calculator = Expression(J_vector[1].dx(0)-J_vector[0].dx(1), Jcurl.element.interpolation_points())
#Jecurl_flux_calculator = Expression(curl(as_vector([J_vector[0],J_vector[1],0])), Jcurl.element.interpolation_points())

Traceback (most recent call last):
File "/root/shared/debug_discourse/curl_MWE.py", line 181, in <module>
Jecurl.interpolate(Jecurl_flux_calculator)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 397, in interpolate
_interpolate(u, cells)
File "/usr/lib/python3.10/functools.py", line 889, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 388, in _
self._cpp_object.interpolate(expr._cpp_object, cells)
RuntimeError: Function value size not equal to Expression value size


The last commented lines all yield the same error except the very last, whose traceback is:

Traceback (most recent call last):
File "/root/shared/debug_discourse/curl_MWE.py", line 183, in <module>
Jecurl_flux_calculator = Expression(curl(as_vector([J_vector[0],J_vector[1],0])), Jcurl.element.interpolation_points())
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 149, in __init__
self._ufcx_expression, module, self._code = jit.ffcx_jit(comm, (e, _X),
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 210, in ffcx_jit
r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 247, in compile_expressions
raise e
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 238, in compile_expressions
impl = _compile_objects(decl, expressions, expr_names, module_name, p, cache_dir,
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 95, in analyze_ufl_objects
processed_expression = _analyze_expression(original_expression, options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/analysis.py", line 117, in _analyze_expression
expression = ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(expression)
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 148, in apply_algebra_lowering
return map_integrand_dags(LowerCompoundAlgebra(), expr)
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in map_integrand_dags
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 66, in map_integrands
return function(integrand)
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/map_integrands.py", line 73, in <lambda>
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 34, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress,
File "/usr/local/lib/python3.10/dist-packages/ufl/corealg/map_dag.py", line 100, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 142, in curl
return as_vector((c(1, 2), c(2, 0), c(0, 1)))
File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/apply_algebra_lowering.py", line 135, in c
return a[j].dx(i) - a[i].dx(j)
File "/usr/local/lib/python3.10/dist-packages/ufl/exproperators.py", line 435, in _dx
File "/usr/local/lib/python3.10/dist-packages/ufl/differentiation.py", line 252, in __new__
dim = find_geometric_dimension(f)
File "/usr/local/lib/python3.10/dist-packages/ufl/domain.py", line 396, in find_geometric_dimension
raise ValueError("Cannot determine geometric dimension from expression.")
ValueError: Cannot determine geometric dimension from expression.


You are using a 2D mesh.
The ufl curl operator is only defined in 3D, as written in the documentation:
https://docs.fenicsproject.org/ufl/2024.1.0/manual/form_language.html#curl-and-rot

See for instance

or Vector curl in 2D - #2 by nate
for solutions