Something is definitely wrong. I would really appreciate some help/guidance in debugging.
Here are 2 screenshots showing the problem. This represents a heating (called Bridgman heat), it is a scalar field computed as -T\left( S_{xx}\frac{\partial J_x}{\partial x} + S_{yy}\frac{\partial J_y}{\partial y}\right) where the current density is given in my OP. I think this heating term is computed correctly almost everywhere, except in some elements of the meshing. This messes up my computations, I think (for example when I integrate this quantity over the whole volume).
import numpy as np
from dolfinx import log
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile, gmshio
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)
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio, VTXWriter, XDMFFile
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/half_annulus.msh", MPI.COMM_WORLD, gdim=2)
inj_current_curve = 12
out_current_curve = 13
reading_voltage_surface_0 = 12
reading_voltage_surface_1 = 13
# Define ME function space
el = FiniteElement("CG", mesh.ufl_cell(), 3)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)
u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)
rho = 1
sigma = 1.0 / rho
κ = 1.8
S_xx = 3000e-6
S_yy = S_xx * 10#* 3#* 300
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)
T_cold = 300.0
ΔT = 2000
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)
the_current = 1e-15
J = the_current / assemble_scalar(form(1 * ds(inj_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 * ds(out_current_curve) + v * J * ds(inj_current_curve)
weak_form = F_T + F_V
# Solve the PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
Jac = derivative(weak_form,TempVolt,dTV)
# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J = Jac)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-10
solver.report = True
solver.max_it = 10000
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "ksp"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
#opts[f"{option_prefix}ksp_max_it"]= 10000
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(TempVolt)
assert(converged)
print(f'''------- Processing --------
Number of interations: {n:d}''')
# Compute fluxes on boundaries
n = FacetNormal(mesh)
down_heat_flux = form(dot(-κ * grad(temp) + Seebeck_tensor * temp * J_vector, n)*ds(out_current_curve))
print(f'''--------- Post processing -------
Computed generated Joule heat through the whole volume: {assemble_scalar(form(dot(rho * J_vector, J_vector) * dx))} W.
Make sure div(J_e) = 0 throughought the whole volume: {assemble_scalar(form((J_vector[0].dx(0) + J_vector[1].dx(1)) * dx))}
Computed generated Bridgman heat through the whole volume: {assemble_scalar(form(- temp * (S_xx * J_vector[0].dx(0) + S_yy * J_vector[1].dx(1)) * dx))} W.
''')
# Current density.
J = VectorFunctionSpace(mesh, ("DG", 2))
Je = Function(J)
Je_flux_calculator = Expression(-sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp), J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
# Bridgman heat.
Q_B = FunctionSpace(mesh, ("DG", 2))
QB = Function(Q_B)
QB_flux_calculator = Expression(- temp * (S_xx * J_vector[0].dx(0) + S_yy * J_vector[1].dx(1)), Q_B.element.interpolation_points())
QB.interpolate(QB_flux_calculator)
Output_0 = TempVolt.sub(0).collapse()
Output_1 = TempVolt.sub(1).collapse()
output_4 = Je
output_5 = QB
with VTXWriter(MPI.COMM_WORLD, "results/temperature.bp", [Output_0], engine="BP4") as vtx:
vtx.write(0.0)
with VTXWriter(MPI.COMM_WORLD, "results/voltage.bp", [Output_1], engine="BP4") as vtx:
vtx.write(0.0)
with VTXWriter(MPI.COMM_WORLD, "results/current_density.bp", [output_4], engine="BP4") as vtx:
vtx.write(0.0)
with VTXWriter(MPI.COMM_WORLD, "results/Bridgman_heat.bp", [output_5], engine="BP4") as vtx:
vtx.write(0.0)
The geo file used to produce the mesh:
h = 0.00915;
r1 = 0.3;
r2 = 0.1;
Point(0) = {r1,0,0,h};
Point(1) = {0,0,0,h};
Point(2) = {r2,0,0,h};
Point(3) = {2*r1-r2,0,0,h};
Point(4) = {2*r1,0,0,h};
Line(1) = {1,2};
Circle(2) = {2,0,3};
Line(3) = {3,4};
Circle(4) = {4,0,1};
Line Loop(1) = {1,2,3,4};
Plane Surface(1) = {1};
Physical Curve("bottom_left", 12) = {1};
Physical Curve("bottom_right", 13) = {3};
Physical Surface("material", 11) = {1};
(discourse website being too limited to accept the content of the file.msh due to a line number limit.)