I have a code that seems to be doing what I seek: it solves for the temperature and voltage in a material where thermoelectricity is taken into account. It returns correct results for some cases I could check with analytical solution, so I tend to trust it in general. The equations to be solved are \nabla \cdot \vec J_e =0 where \vec J_e = \sigma \nabla V \sigma S \nabla T for the voltage V and the equation \nabla \cdot (\kappa \nabla T) +\rho \vec{J_e}^2=0 for the temperature T. Where \rho=1/\sigma and S is the Seebeck coefficient.

One of my goals now is to use the code on different meshes for example, without tweaking the code. What would be a nice way to accomplish this task? Make a function out of my whole code, where as input I would manually specify the physical boundaries where boundary conditions are applied? It would become pretty complex if it could deal with 1D, 2D and 3D meshes without crashing.

As it is currently, Newtonâ€™s method fails to converge if I modify Seebeck coefficient to values higher than about 30e6. However in real materials this number can reach 1000e6, so I would like to see what would happen in that case, but I am unable to make the Newtonâ€™s method converge. Any pointer is welcome.

I wanted to specify an electric current magnitude (say 1 ampere), which I used to do by applying a 0 voltage Dirichlet boundary condition on one surface (where the current is injected), and by adding a term in the weak form for the surface where the current exits the material (Neumannâ€™s boundary condition). It used to work well, however it doesnâ€™t seem to work here. I use mixed elements, so I suspect this is the cause, but still, this sounds strange to me. To fix the problem, I considered pure Neumannâ€™s boundary conditions, so I donâ€™t use any Dirichlet boundary conditions for the voltage. As a result, the code seems to find some â€śzeroâ€ť of voltage by itself. Voltage is usually defined up to a constant. However the solution I get if I set a boundary to a constant is wrong. But when I let FEniCSx decide which constant it is, it works fine. I do not understand why it is the case. Any pointer is welcome. It is as if FEniCSx thinks the voltage isnâ€™t defined up to a constant.

I would like to reuse the code with some modularity. For example, I would like to have the choice to ignore thermoelectricity (so by setting Seebeck coefficient to 0). I also want to be able to ignore Joule heating. I can easily do this if I define a function of my whole code and where I have a boolean variable for Joule_heat and where if Joule_heat ==True, then I would include this term in the weak form, otherwise set it to 0 or do nothing. But, as a non programmer, would this be a nice way to achieve this? I realize I do not use any class nor dictionary in my whole code. I am surely missing improvements.
My whole code is here:
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
# Read the mesh files produced above.
from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(mesh, name="Grid")
# Read the mesh file created by Gmsh.
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/long_wire_with_contacts.msh", MPI.COMM_WORLD, gdim=3)
# Convert the mesh file to XDMF using meshio.
proc = MPI.COMM_WORLD.rank
#print(proc)
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
if proc == 0:
# Read in mesh
msh = meshio.read("meshes/long_wire_with_contacts.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
tetrahedron_mesh = create_mesh(msh, "tetra", prune_z=False)
meshio.write("mt.xdmf", triangle_mesh)
meshio.write("mesh.xdmf", tetrahedron_mesh)
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
dx, grad, inner, as_tensor, inv, MixedElement, FiniteElement, split)
import numpy as np
# Define ME function space
el = FiniteElement("CG", mesh.ufl_cell(), 1)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)
u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)
rho = 1e9
sigma = 1.0 / rho
kappa = 144
Seebeck_coeff = 0#30e6
# Reminder of the mesh:
#Physical Volume("the_wire", 13) = {1};
#Physical Surface("left_contact_2", 28) = {17};
#Physical Surface("right_end", 47) = {26};
#Physical Surface("right_contact_2", 48) = {23};
#Physical Surface("left_end", 49) = {25};
inj_current_surface = 28
vanish_voltage_surface = 48 # This corresponds to the curve the current leaves the material.
V_current_contact_out = 0.0 # Voltage value of the curve where the current leaves the material.
#V_current_contact_in = 1.6
reading_voltage_surface_0 = 28
reading_voltage_surface_1 = 48
hot_surface = 47
cold_surface = 49
T_cold = 300.0
T_hot = 320.0
# Define the boundary conditions. For this, we need to find the facets where the b.c. will be applied.
# I use Dirichlet b.c. where I inject the current as well as where the current exists the material,
# except in the 1st iteration.
from petsc4py.PETSc import ScalarType
left_facets_current = ft.find(inj_current_surface)
right_facets_current = ft.find(vanish_voltage_surface)
left_dofs_current = locate_dofs_topological(ME.sub(1), mesh.topology.dim1, right_facets_current)
bc_volt = dirichletbc(ScalarType(V_current_contact_out), left_dofs_current, ME.sub(1))
left_facets_temperature = ft.find(hot_surface)
right_facets_temperature = ft.find(cold_surface)
left_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim1, left_facets_temperature)
right_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim1, right_facets_temperature)
bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(T_hot), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right]#, bc_volt]
from dolfinx.fem import assemble_scalar, form
from ufl import Measure
x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh,subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)
the_current = 1 # Current, in amperes.
J = the_current / assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh)))
print('Current density:', J)
print('Surface area where current is injected', assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh))))
print('Surface area where current leaves the wire', assemble_scalar(form(1 * ds(vanish_voltage_surface, domain=mesh))))
L = 1e2*0.9
estimated_Joule_heat = the_current ** 2 * rho * L / assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh)))
print('Estimated Joule heat generated through the whole volume: ', str(estimated_Joule_heat) + ' W.')
from ufl import dot, derivative
from petsc4py import PETSc
def temp_init(x):
values = np.full(x.shape[1], T_cold, dtype=PETSc.ScalarType)
return values
def volt_init(x):
values = np.full(x.shape[1], V_current_contact_out, dtype=PETSc.ScalarType)
return values
# Initialize values for the temperature and voltage.
TempVolt.sub(0).interpolate(temp_init)
TempVolt.sub(1).interpolate(volt_init)
# Weak form.
J_vector = sigma * grad(volt)  sigma * Seebeck_coeff * grad(temp)
Fourier_term = dot(kappa * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
F_T = Fourier_term + Joule_term
F_V = dot(grad(v), sigma * grad(volt))*dx dot(grad(v), sigma * Seebeck_coeff * grad(temp))*dx + v * J * ds(vanish_voltage_surface)  v * J * ds(inj_current_surface)
print('This should correspond to the current injected: ', assemble_scalar(form(J*ds(vanish_voltage_surface))))
weak_form = F_T + F_V
Jac = derivative(weak_form,TempVolt,dTV)
# Solve the PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J = Jac)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e6
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
#opts[f"{option_prefix}ksp_type"] = "cg"
#opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"]= 10000
ksp.setFromOptions()
from dolfinx import log
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(TempVolt)
assert(converged)
print(f"Number of interations: {n:d}")
avg_voltage_curve_integral_0 = assemble_scalar(form(volt * ds(reading_voltage_surface_0, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_surface_0, domain=mesh)))
avg_voltage_curve_integral_1 = assemble_scalar(form(volt * ds(reading_voltage_surface_1, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_surface_1, domain=mesh)))
voltage_diff = avg_voltage_curve_integral_1  avg_voltage_curve_integral_0
print(voltage_diff,'voltage measured (V)')
print('voltage left', avg_voltage_curve_integral_0)
print('voltage right', avg_voltage_curve_integral_1)
avg_temperature_curve_integral_1 = assemble_scalar(form(temp * ds(inj_current_surface, domain=mesh))) / assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh)))
print(avg_temperature_curve_integral_1,'Temperature of current injection contact (K)')
print('Computed generated Joule heat through the whole volume: ', str(assemble_scalar(form(dot(rho * J_vector, J_vector) * dx))) + ' W.')
# Need to interpolate E_field on some vector space. Because it wouldn't be continuous otherwise.
from dolfinx.fem import VectorFunctionSpace, Expression
W = VectorFunctionSpace(mesh, ("DG", 0))
E_field = Function(W)
E_field_expr = Expression(grad(TempVolt.sub(1)), W.element.interpolation_points())
E_field.interpolate(E_field_expr)
# Check the current density.
W_current = VectorFunctionSpace(mesh, ("DG", 1))
Je = Function(W_current)
Je_expr = Expression(sigma * grad(TempVolt.sub(1))  sigma * Seebeck_coeff * grad(TempVolt.sub(0)), W_current.element.interpolation_points())
Je.interpolate(Je_expr)
# Check the heat flux.
W_heat_flux = VectorFunctionSpace(mesh, ("DG", 1))
JQ = Function(W_heat_flux)
JQ_expr = Expression(kappa * grad(TempVolt.sub(0)) + Seebeck_coeff * TempVolt.sub(0) * Je, W_heat_flux.element.interpolation_points())
JQ.interpolate(JQ_expr)
# Check the temperature gradient, which is the "hot to cold" direction.
W_gradT = VectorFunctionSpace(mesh, ("DG", 1))
gradT = Function(W_gradT)
gradT_expr = Expression(grad(TempVolt.sub(0)), W_gradT.element.interpolation_points())
gradT.interpolate(gradT_expr)
# Check grad(volt), which is proportional to grad mu
W_gradV = VectorFunctionSpace(mesh, ("DG", 1))
gradV = Function(W_gradV)
gradV_expr = Expression(grad(TempVolt.sub(1)), W_gradV.element.interpolation_points())
gradV.interpolate(gradV_expr)
from dolfinx import io
with io.XDMFFile(mesh.comm, "solution/Efield_S30e6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(E_field)
with io.XDMFFile(mesh.comm, "solution/Je_S30e6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(Je)
with io.XDMFFile(mesh.comm, "solution/JQ_S30e6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(JQ)
with io.XDMFFile(mesh.comm, "solution/gradT_S30e6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(gradT)
with io.XDMFFile(mesh.comm, "solution/voltage_S30e6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(TempVolt.sub(1))
with io.XDMFFile(mesh.comm, "solution/temperature_S30e6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(TempVolt.sub(0))
with io.XDMFFile(mesh.comm, "solution/gradV_S30e6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(gradV)