Lost with boundary conditions in ME spaces

My code seems to be doing what I wish, at first glance at least, but I do not understand why (and it may therefore not correspond to what I wish).
I solve 2 coupled PDE’s, the unknowns being “temp” and “volt” in a thermoelectric problem. It is question of a unit square material with anisotropic Seebeck coefficient S, and nothing depends on temperature.
The equations solved are the heat equation \nabla \cdot (\kappa \nabla T)+\rho |\vec J|^2-T\left( S_{xx}\frac{\partial J_x}{\partial x} + S_{yy}\frac{\partial J_y}{\partial x}\right)=0, and Laplace equation that comes out of the steady-state condition \nabla \cdot \vec J =0.

Physically, if this material is placed between 2 reservoirs kept at different temperatures, internal electric currents are going to appear. I am trying to simulate those, and it looks like FEniCSx is able to get those right. These currents should have no component perpendicular to any surface (or sides in 2D), as it would imply a charge transfer with the surroundings, which is forbidden (the material is supposed to be in deep space, aside the 2 thermal reservoirs).
And this is exactly what I get with FEniCSx.
However, I apply 2 Dirichlet b.c. for the temperature (the 2 reservoirs keeping 2 sides at constant temperature), and I apply no conditions on the voltage, which mean vanishing Neumann b.c. for the voltage on all surfaces.

Here is my problem:
Since \vec J=-\sigma \nabla V -\sigma S \nabla T, I am unable to understand how FEniCSx is able to determine that \vec J is parallel to the sides on the 2 sides that have the thermal Dirichlet b.c. applied. This should not be the case, because on these sides \nabla V is not supposed to have any component perpendicular to the sides, yet \nabla T does. Therefore I expect FEniCSx to get \vec J have a perpendicular component to those surfaces due to the thermal gradient. However, this is not the case, and FEniCSx is able to get the physically correct \vec J (at least it looks like so) always parallel to the sides.

I am lost on how this is possible.

The full code is:

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, 40, 40, 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)
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)
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]

# Mesh tags to apply Neumann b.c.
tdim = domain.topology.dim
fdim = tdim -1 
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)))

# 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

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])

# 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
solver.atol = 1e-8
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"
n, converged = solver.solve(TempVolt)
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())

with VTXWriter(MPI.COMM_WORLD, "results/Je_debug.bp", Je, engine="BP4") as vtx:
# Fourier heat flux.    
Q = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
q = Function(Q)
flux_calculator = Expression(-κ * grad(temp), Q.element.interpolation_points())
with VTXWriter(MPI.COMM_WORLD, "results/Fourier_heat_flux_debug.bp", q, engine="BP4") as vtx:

# Grad V.
Q2 = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
q2 = Function(Q2)
flux_calculator = Expression(grad(volt), Q2.element.interpolation_points())
with VTXWriter(MPI.COMM_WORLD, "results/grad_volt_debug.bp", q2, engine="BP4") as vtx:

Examining the 2 output files with Paraview reveals that \nabla T looks normal, it is parallel to all sides except those 2 where the temperature is kept fixed with Dirichlet b.c., grad T is perpendicular to those surfaces as it should.

The gradient of the voltage looks wrong, in that it is not compatible with vanishing Neumann b.c.s. It is parallel to the sides (as it should) except for the 2 sides where the thermal Dirichlet b.c.s are applied, where it isn’t entirely perpendicular either, but has a changing direction.

It really looks like “magic” in that it correctly gets the physical \vec J but it shouldn’t, since I haven’t specified any b.c.s for volt.
I attach the picture of grad V. If someone can tell me why it isn’t always parallel to the sides, this would be a big first step with my understanding of the situation.

Hi people, let me know if my question is unclear or if there is something that prevents an answer.

The question can be recast to ‘‘How to apply vanishing Neumann b.c.s to J_vector?’’, taking into account the Dirichlet b.c.s for ‘‘temp’’.

The only explanation I can think of, other than magic, is that one of the PDE I solve is div (J) equals 0. And so the natural Neumann b.c. would be related to J because of this equation.

It means natural Neumann b.c.s aren’t directly related to the fluxes of the unknowns that are solved for, but to a flux that appears in the PDE(s). That looks very strange to me, but I don’t have other ideas.

I’ve given more thoughts about this observation. My current conclusion is that the claim that vanishing Neumann boundary conditions are applied by default if no boundary conditions are explicitly specified, for a given finite element function, is not correct. And the above MWE is indeed a case where this happens.

If you go through the hand derivation of the weak form of the PDE \nabla \cdot \vec J (equivalently, if you integrate this equation over the whole volume and then apply Gauss theorem), you’ll find a different expression than the one I posted in my MWE. In fact, in my MWE, it is exactly the same problem except that I enforce a given electrical current. By doing so, \nabla V has to adjust such that the input and output currents match what I prescribe on the boundaries. And \nabla V does not vanish on all boundaries, in particular those where Dirichlet boundary conditions are specified for the other unknown (temperature in this case).

This also responds to Implementing a potential problem with a current condition - #21 by Daniel.Baker. I do not believe prescribing the current on the 2 boundaries (in and out current) is insufficient to fully specify the voltage (scalar field). I believe it is perfectly fine, keeping in mind that voltage (or electrostatics potential) is defined up to a constant.

I would appreciate if @kamensky could comment.