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]])
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
(...)
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
d = Grad(d)
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.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
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.
