Comment my code + many questions

General
1

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.

  1. 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.

  2. As it is currently, Newton’s method fails to converge if I modify Seebeck coefficient to values higher than about 30e-6. However in real materials this number can reach 1000e-6, 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.

  3. 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.

  4. 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.dim-1)
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 = 1e-9
sigma = 1.0 / rho

kappa = 144

Seebeck_coeff = 0#30e-6

# 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.dim-1, 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.dim-1, left_facets_temperature)
right_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, 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 = 1e-2*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 = 1e-6
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_S30e-6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(E_field)

with io.XDMFFile(mesh.comm, "solution/Je_S30e-6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(Je)

with io.XDMFFile(mesh.comm, "solution/JQ_S30e-6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(JQ)

with io.XDMFFile(mesh.comm, "solution/gradT_S30e-6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(gradT)

with io.XDMFFile(mesh.comm, "solution/voltage_S30e-6_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_S30e-6_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_S30e-6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(gradV)
2

I would suggest using class structures.
Then you can have different functions applying appropriate boundary conditions.
Also keep in mind that the only thing that should differ in the core of your variational formulation is the mesh/function space, the UFL forms should be the same in 1D/2D/3D.

Oasisx (disclaimer Im the author) GitHub - ComputationalPhysiology/oasisx: Next generation Open Source Navier Stokes solver using FEniCSx
is an example of how to re-use the variational formulation for different problems.

See for instance: Default absolute tolerance and relative tolerance - #4 by nate

I would strongly suggest you try to simplify your problem to illustrate this issue.
If you won’t simplify it, you should really be careful about the integration by parts of your problem, making sure that you are not missing a sign in the transformations. Additionally, when you say the code seem to find some “zero” voltage by itself, you should look at the equation you are now solving.

Again, using classes, and boolean inputs to either:

  1. Add terms to your variational form
  2. Add a c=dolfinx.fem.Constant(mesh, value) to your variational form (i.e. F += c*joule_heating_terms), where value is 0 or 1 depending on if you want to add the effect or not.

I do not have time to go through your code in detail, but I would strongly suggest you to start small, and take out simple parts of your code that you can re-use, and make up a general structure as you go.

1 Like
3

Thank you very much for all the tips and your reply.
I will comment to this first:

I would strongly suggest you try to simplify your problem to illustrate this issue.
If you won’t simplify it, you should really be careful about the integration by parts of your problem, making sure that you are not missing a sign in the transformations. Additionally, when you say the code seem to find some “zero” voltage by itself, you should look at the equation you are now solving.

For the voltage equation, \nabla \cdot \vec{J_e}=0 \Rightarrow \nabla \cdot (-\sigma \nabla V) - \nabla \cdot (\sigma S \nabla T)=0. Multiplying by the test function v and integrating over the domain, I get \int_{\Omega} \nabla \cdot (-\sigma \nabla V)vd\Omega-\int_{\Omega} \nabla \cdot (-S \sigma \nabla T)vd\Omega=0. Integrating by parts, I get -\sigma \nabla V v|_{\partial \Omega}-\int_\Omega -\sigma \nabla V \cdot \nabla v d\Omega-\sigma S \nabla T v |_{\partial \Omega}-\int_\Omega \sigma S \nabla T \cdot \nabla v d\Omega=0 which is the weak form I should implement in FEniCSx. I think this is what I implemented. In my case, I made a correspondance between \sigma \nabla V v|_{\partial \Omega} and v \vec J \cdot d\vec A. Essentially I am assuming that \sigma \nabla V gives the current density I inject. I am not sure this is correct… since \vec J_e =-\sigma \nabla V -\sigma S \nabla T. This seems this would hold only in the case where S=0, i.e. when the Seebeck effect is neglected. Maybe I should just let the term \sigma \nabla V v|_{\partial \Omega} as is, in the weak form. Do you have an idea if that would make sense?

Prior to that, I would use a Dirichlet boundary condition for one of the 2 surfaces, and Neumann’s boundary condition for the other surface. The current I would get through the “wire” used to look correct. (It also looks correct now with my conversion to pure Neumann b.c’s).

4

I have refactored my code into something that looks much simpler, but I am still not using classes (I’m not a well versed programmer but I’ll do my best).

I face an error with my latest change. Here is my full code:

from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio, XDMFFile
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological, VectorFunctionSpace, Expression)
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner, as_tensor, inv, MixedElement, FiniteElement, split, Measure, dot, derivative)
import numpy as np
from petsc4py.PETSc import ScalarType
from ufl import Measure
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
from dolfinx import log, io

def thermoelectric_solver(mesh_file_name, physical_param, physical_surfaces_ID_numbers, temperature_of_surfaces):
    '''
    mesh_file_name: string of the mesh file name.
    physical_param: list or tuple containing values for rho, kappa and S (they can be tensors or floats), and the current amplitude in amperes (float).
    physical_surfaces_ID_numbers: List of IDs of the physical surfaces defined by Gmsh. (int)
    temperature_of_surfaces: List containing the values of T_cold and T_hot (floats).   
    '''
    
    # Read the mesh file created by Gmsh.
    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.dim-1)
    with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
        ft = xdmf.read_meshtags(mesh, name="Grid")
    # Define ME function space
    el = FiniteElement("CG", mesh.ufl_cell(), 2)
    mel = MixedElement([el, el])
    ME = FunctionSpace(mesh, mel)

    u, v = TestFunction(ME)
    TempVolt = Function(ME)
    temp, volt = split(TempVolt)
    dTV = TrialFunction(ME)

    rho = physical_param[0]
    sigma = 1.0 / rho
    kappa = physical_param[1]
    Seebeck_coeff = physical_param[2]
    current_amplitude = physical_param[3]
    T_cold = temperature_of_surfaces[0]
    T_hot = temperature_of_surfaces[1]
    
    inj_current_surface = physical_surfaces_ID_numbers[0]
    out_current_surface = physical_surfaces_ID_numbers[1]
    reading_voltage_surface_0 = physical_surfaces_ID_numbers[2]
    reading_voltage_surface_1 = physical_surfaces_ID_numbers[3]
    hot_surface = physical_surfaces_ID_numbers[4]
    cold_surface = physical_surfaces_ID_numbers[5]
    
    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.dim-1, left_facets_temperature)
    right_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, 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]
    
    x = SpatialCoordinate(mesh)
    dx = Measure("dx", domain=mesh,subdomain_data=ct)
    ds = Measure("ds", domain=mesh, subdomain_data=ft)
    the_current = current_amplitude#1 # Current, in amperes.
    J = the_current / assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh)))
    print('Current density at entry surface:', str(J) + ' A / m^2.')
    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(out_current_surface, domain=mesh))))
    L = 1e-2*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.')
    

    # 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(out_current_surface) - v * J * ds(inj_current_surface) 
    print('This should correspond to the current injected: ', assemble_scalar(form(J*ds(out_current_surface))))

    weak_form = F_T + F_V

    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-6
    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()
    

    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_surface_integral_2 = assemble_scalar(form(temp * ds(inj_current_surface, domain=mesh))) / assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh)))
    print(avg_temperature_surface_integral_2,'Temperature of current injection contact (K)')
    avg_temperature_surface_integral_3 = assemble_scalar(form(temp * ds(out_current_surface, domain=mesh))) / assemble_scalar(form(1 * ds(out_current_surface, domain=mesh)))
    print(avg_temperature_surface_integral_3,'Temperature of current outjection 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.
    W = VectorFunctionSpace(mesh, ("DG", 1))
    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)
    

    with io.XDMFFile(mesh.comm, "solution/Pack_of_solutions_S30e-6_I0_newWF_test_noVoltbc.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(E_field)
        xdmf.write_function(Je)
        xdmf.write_function(JQ)
        xdmf.write_function(gradT)
        xdmf.write_function(TempVolt.sub(1))
        xdmf.write_function(TempVolt.sub(0))
        xdmf.write_function(gradV)
    print('Problem solved.')

thermoelectric_solver("meshes/long_wire_with_contacts.msh", [1e-9, 144, 30e-6, 1.0], [28, 48, 28, 48, 47, 49], [300.0, 320.0])

The error occurs since I try to write all the solutions into a single file, for Paraview. I just followed your tutorial @dokken (specifically the bottom of Implementation — FEniCSx tutorial). The error I get is:

HDF5-DIAG: Error detected in HDF5 (1.12.2) MPI-process 0:
  #000: /tmp/hdf5-1.12.2/src/H5D.c line 141 in H5Dcreate2(): unable to create dataset
    major: Dataset
    minor: Unable to create file
  #001: /tmp/hdf5-1.12.2/src/H5VLcallback.c line 1806 in H5VL_dataset_create(): dataset create failed
    major: Virtual Object Layer
    minor: Unable to create file
  #002: /tmp/hdf5-1.12.2/src/H5VLcallback.c line 1771 in H5VL__dataset_create(): dataset create failed
    major: Virtual Object Layer
    minor: Unable to create file
  #003: /tmp/hdf5-1.12.2/src/H5VLnative_dataset.c line 73 in H5VL__native_dataset_create(): unable to create dataset
    major: Dataset
    minor: Unable to initialize object
  #004: /tmp/hdf5-1.12.2/src/H5Dint.c line 396 in H5D__create_named(): unable to create and link to dataset
    major: Dataset
    minor: Unable to initialize object
  #005: /tmp/hdf5-1.12.2/src/H5L.c line 1767 in H5L_link_object(): unable to create new link to object
    major: Links
    minor: Unable to initialize object
  #006: /tmp/hdf5-1.12.2/src/H5L.c line 2009 in H5L__create_real(): can't insert link
    major: Links
    minor: Unable to insert object
  #007: /tmp/hdf5-1.12.2/src/H5Gtraverse.c line 837 in H5G_traverse(): internal path traversal failed
    major: Symbol table
    minor: Object not found
  #008: /tmp/hdf5-1.12.2/src/H5Gtraverse.c line 613 in H5G__traverse_real(): traversal operator failed
    major: Symbol table
    minor: Callback failed
  #009: /tmp/hdf5-1.12.2/src/H5L.c line 1802 in H5L__link_cb(): name already exists
    major: Links
    minor: Object already exists

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In [6], line 1
----> 1 thermoelectric_solver("meshes/long_wire_with_contacts.msh", [1e-9, 144, 30e-6, 1.0], [28, 48, 28, 48, 47, 49], [300.0, 320.0])

Cell In [5], line 150, in thermoelectric_solver(mesh_file_name, physical_param, physical_surfaces_ID_numbers, temperature_of_surfaces)
    148 xdmf.write_mesh(mesh)
    149 xdmf.write_function(E_field)
--> 150 xdmf.write_function(Je)
    151 xdmf.write_function(JQ)
    152 xdmf.write_function(gradT)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/utils.py:150, in XDMFFile.write_function(self, u, t, mesh_xpath)
    149 def write_function(self, u, t: float = 0.0, mesh_xpath="/Xdmf/Domain/Grid[@GridType='Uniform'][1]"):
--> 150     super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)

RuntimeError: Failed to create HDF5 global dataset.

Edit: I fixed my problem by giving a name (as you did in the tutorial), to Je, JQ and so on.

5

There is something I don’t understand about the voltage solution. As my code is, I do not specify any Dirichlet boundary condition for the voltage (there are terms related to Neumann’s b.c. in the weak form though). The solution I get seems to make sense, but I do not understand how FEniCSx chooses the 0 of voltage (voltage is defined up to a constant so in theory if I added any constant to the voltage solution, it should be fine).

If I manually set a constant voltage surface (as should be the case for the 2 surfaces where Neumann’s b.c.s are applied), then I get a wrong solution, it doesn’t just shift the solution I get when I don’t specify Dirichlet’s b.c.s.

Any pointer is welcome.

6

Well, I asked chatGPT. It says FEniCS (specifically the solver) would choose an arbitrary “zero voltage” point in this case (only applied Neumann b.c. for the voltage, and no Dirichlet b.c.). As long as this zero voltage point doesn’t lie on the surfaces where I apply the Neumann b.c., it should be fine.

In FEniCSx, the reference point for the voltage is usually not specified explicitly but rather determined by the linear solver. When solving a linear system of equations, such as the one arising from the finite element discretization of the Poisson equation, the solver will typically choose a particular node or degree of freedom as the reference point for the solution. This choice is usually based on the specific implementation of the solver and may vary from one solver to another.

If you want to know which point is being used as the reference for the voltage solution in your particular case, one option is to print out the solution vector and inspect it. Since the voltage is defined up to a constant, you can simply add or subtract a constant value to the solution vector and see how the solution changes. By doing this for different values of the constant, you may be able to identify which degree of freedom is being used as the reference. Another option is to consult the documentation or source code of the linear solver you are using to see if it specifies how the reference point is chosen.

I hope it’s correct.

7

Hello @dokken , would you have an idea on whether using Lagrange elements converge towards the real solution to the PDEs in this case? One of them, div(Je), derives from Maxwell equations. Or whether a Nedelec approach is the correct one?