Hello,
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)
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]
# 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)
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])
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
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"
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)
# Fourier heat flux.
Q = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
q = Function(Q)
flux_calculator = Expression(-κ * grad(temp), Q.element.interpolation_points())
q.interpolate(flux_calculator)
with VTXWriter(MPI.COMM_WORLD, "results/Fourier_heat_flux_debug.bp", q, engine="BP4") as vtx:
vtx.write(0.0)
# Grad V.
Q2 = functionspace(domain, ("DG", deg-1, (domain.geometry.dim,)))
q2 = Function(Q2)
flux_calculator = Expression(grad(volt), Q2.element.interpolation_points())
q2.interpolate(flux_calculator)
with VTXWriter(MPI.COMM_WORLD, "results/grad_volt_debug.bp", q2, engine="BP4") as vtx:
vtx.write(0.0)
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.
