Mixed elements, solution insensitive to function of other subspace

General
1

I am trying to convert some old (FEniCS) working code to the newer FEniCSx, but I am facing a problem. The temperature distribution I get is insensitive to the electric current, but this should not be the case, since I consider Joule heating.

In terms of code, I employ mixed elements and solve the heat equation (to get the temperature distribution), and div(J)=0, where J = sigma E. I expect a parabolic solution for the temperature distribution along z. As boundary conditions, I keep fixed the temperatures of the left and right ends of a wire (parallelepipede).
My full code starts with the file.msh:

$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
3
2 14 "left_end"
2 15 "right_end"
3 13 "the_wire"
$EndPhysicalNames
$Entities
8 12 6 1
1 0 0 1 0 
2 0 0 0 0 
3 0 0.01 1 0 
4 0 0.01 0 0 
5 0.01 0 1 0 
6 0.01 0 0 0 
7 0.01 0.01 1 0 
8 0.01 0.01 0 0 
1 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 0 2 2 -1 
2 -1e-07 -1.000000000002735e-07 0.9999999000000001 1e-07 0.0100001 1.0000001 0 2 1 -3 
3 -1e-07 0.009999900000000001 -9.999999994736442e-08 1e-07 0.0100001 1.0000001 0 2 4 -3 
4 -1e-07 -1.000000000002735e-07 -1e-07 1e-07 0.0100001 1e-07 0 2 2 -4 
5 0.009999900000000001 -1e-07 -9.999999994736442e-08 0.0100001 1e-07 1.0000001 0 2 6 -5 
6 0.009999900000000001 -1.000000000002735e-07 0.9999999000000001 0.0100001 0.0100001 1.0000001 0 2 5 -7 
7 0.009999900000000001 0.009999900000000001 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 0 2 8 -7 
8 0.009999900000000001 -1.000000000002735e-07 -1e-07 0.0100001 0.0100001 1e-07 0 2 6 -8 
9 -1.000000000002735e-07 -1e-07 -1e-07 0.0100001 1e-07 1e-07 0 2 2 -6 
10 -1.000000000002735e-07 -1e-07 0.9999999000000001 0.0100001 1e-07 1.0000001 0 2 1 -5 
11 -1.000000000002735e-07 0.009999900000000001 -1e-07 0.0100001 0.0100001 1e-07 0 2 4 -8 
12 -1.000000000002735e-07 0.009999900000000001 0.9999999000000001 0.0100001 0.0100001 1.0000001 0 2 3 -7 
1 -1e-07 -1.000000000002735e-07 -9.999999994736442e-08 1e-07 0.0100001 1.0000001 0 4 -1 4 3 -2 
2 0.009999900000000001 -1.000000000002735e-07 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 0 4 -5 8 7 -6 
3 -1.000000000002735e-07 -1e-07 -9.999999994736442e-08 0.0100001 1e-07 1.0000001 0 4 -9 1 10 -5 
4 -1.000000000002735e-07 0.009999900000000001 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 0 4 -11 3 12 -7 
5 -1.000000000002735e-07 -1.000000000002735e-07 -1e-07 0.0100001 0.0100001 1e-07 1 14 4 -4 9 8 -11 
6 -1.000000000002735e-07 -1.000000000002735e-07 0.9999999000000001 0.0100001 0.0100001 1.0000001 1 15 4 -2 10 6 -12 
1 -1.000000000002735e-07 -1.000000000002735e-07 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 1 13 6 -1 2 -3 4 -5 6 
$EndEntities
$Nodes
15 44 1 44
0 1 0 1
1
0 0 1
0 2 0 1
2
0 0 0
0 3 0 1
3
0 0.01 1
0 4 0 1
4
0 0.01 0
0 5 0 1
5
0.01 0 1
0 6 0 1
6
0.01 0 0
0 7 0 1
7
0.01 0.01 1
0 8 0 1
8
0.01 0.01 0
1 1 0 9
9
10
11
12
13
14
15
16
17
0 0 0.09999999999999987
0 0 0.1999999999999997
0 0 0.2999999999999997
0 0 0.3999999999999997
0 0 0.4999999999999998
0 0 0.5999999999999998
0 0 0.6999999999999998
0 0 0.7999999999999999
0 0 0.8999999999999999
1 3 0 9
18
19
20
21
22
23
24
25
26
0 0.01 0.09999999999999987
0 0.01 0.1999999999999997
0 0.01 0.2999999999999997
0 0.01 0.3999999999999997
0 0.01 0.4999999999999998
0 0.01 0.5999999999999998
0 0.01 0.6999999999999998
0 0.01 0.7999999999999999
0 0.01 0.8999999999999999
1 5 0 9
27
28
29
30
31
32
33
34
35
0.01 0 0.09999999999999987
0.01 0 0.1999999999999997
0.01 0 0.2999999999999997
0.01 0 0.3999999999999997
0.01 0 0.4999999999999998
0.01 0 0.5999999999999998
0.01 0 0.6999999999999998
0.01 0 0.7999999999999999
0.01 0 0.8999999999999999
1 7 0 9
36
37
38
39
40
41
42
43
44
0.01 0.01 0.09999999999999987
0.01 0.01 0.1999999999999997
0.01 0.01 0.2999999999999997
0.01 0.01 0.3999999999999997
0.01 0.01 0.4999999999999998
0.01 0.01 0.5999999999999998
0.01 0.01 0.6999999999999998
0.01 0.01 0.7999999999999999
0.01 0.01 0.8999999999999999
2 5 0 0
2 6 0 0
3 1 0 0
$EndNodes
$Elements
3 64 1 64
2 5 2 2
1 2 4 6 
2 6 4 8 
2 6 2 2
3 1 5 3 
4 5 7 3 
3 1 4 60
5 9 2 4 6 
6 4 9 6 27 
7 27 6 4 8 
8 4 9 27 18 
9 8 4 27 18 
10 8 27 36 18 
11 10 9 18 27 
12 18 10 27 28 
13 28 27 18 36 
14 18 10 28 19 
15 36 18 28 19 
16 36 28 37 19 
17 11 10 19 28 
18 19 11 28 29 
19 29 28 19 37 
20 19 11 29 20 
21 37 19 29 20 
22 37 29 38 20 
23 12 11 20 29 
24 20 12 29 30 
25 30 29 20 38 
26 20 12 30 21 
27 38 20 30 21 
28 38 30 39 21 
29 13 12 21 30 
30 21 13 30 31 
31 31 30 21 39 
32 21 13 31 22 
33 39 21 31 22 
34 39 31 40 22 
35 14 13 22 31 
36 22 14 31 32 
37 32 31 22 40 
38 22 14 32 23 
39 40 22 32 23 
40 40 32 41 23 
41 15 14 23 32 
42 23 15 32 33 
43 33 32 23 41 
44 23 15 33 24 
45 41 23 33 24 
46 41 33 42 24 
47 16 15 24 33 
48 24 16 33 34 
49 34 33 24 42 
50 24 16 34 25 
51 42 24 34 25 
52 42 34 43 25 
53 17 16 25 34 
54 25 17 34 35 
55 35 34 25 43 
56 25 17 35 26 
57 43 25 35 26 
58 43 35 44 26 
59 1 17 26 35 
60 26 1 35 5 
61 5 35 26 44 
62 26 1 5 3 
63 44 26 5 3 
64 44 5 7 3 
$EndElements

Then my FEniCSx code:

from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("long_wire.msh", MPI.COMM_WORLD, gdim=3)

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("long_wire.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)

# Now 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")

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)

rho = 1e-9
sigma = 1.0 / rho

kappa = 144

# Reminder of the mesh:

#Physical Volume("the_wire", 13) = {1};
#Physical Surface("left_end", 14) = {5};
#Physical Surface("right_end", 15) = {6};

inj_current_surface = 14
vanish_voltage_surface = 15 # 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.
reading_voltage_surface_0 = 14
reading_voltage_surface_1 = 15

# Define the boundary conditions
from petsc4py.PETSc import ScalarType
left_facets = ft.find(inj_current_surface)
right_facets = ft.find(vanish_voltage_surface)
left_dofs = locate_dofs_topological(ME.sub(1), mesh.topology.dim-1, left_facets)
bc_volt = dirichletbc(ScalarType(V_current_contact_out), left_dofs, ME.sub(1))

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)
bc_temp_left = dirichletbc(ScalarType(300.0), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(320.0), 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)
print(x)
#print(x.shape)
dx = Measure("dx", domain=mesh,subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)
the_current = 20.0 # 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))))

from ufl import dot
from petsc4py import PETSc

def temp_init(x):
    values = np.full(x.shape[1], 300, dtype=PETSc.ScalarType) 
    return values

def volt_init(x):
    values = np.full(x.shape[1], 1.0e-3, 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)
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 + v * J * ds(vanish_voltage_surface)
weak_form = F_T + F_V


# Solve the PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-4
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}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}")

lu_solver = ksp
viewer = PETSc.Viewer().createASCII("lu_output.txt")
lu_solver.view(viewer)
solver_output = open("lu_output.txt", "r")

with io.XDMFFile(mesh.comm, "solution/voltage.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(TempVolt.sub(1))

with io.XDMFFile(mesh.comm, "solution/temperature.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(TempVolt.sub(0))

I get a linear temperature distribution, from 300 K to 320 K, which indicates that the boundary conditions are correctly implemented, but not the influence of J on the temperature. I do expect a parabolic dependence, not a linear one. I do not see what is wrong, it is as if “temp” was not influenced by “J_vector” in my code. Any idea what could be wrong?

2

Its quite hard to give any guidance when there is so much different stuff going on in your code, and you have not provided the reference code you have in legacy dolfin.

Have you checked that

This integral is the same in both FEniCSx and legacy fenics?

So could you try to simplify the problem? what happens if you remove the J term in both legacy fenics and FEniCSx?

1 Like
3

Thanks for the support @dokken.
My original code is even more complex (rho and kappa are tensors, my mesh is much more complex, and I have other terms in the weak form, but I remember than I tested every single term effects in the weak form, and I am 100% sure that the Joule effect had the effect of a warming of the wire, and that it lead to a parabolic-like shape for T(z), which makes physical sense).

Here is the evaluation of the integral in my new code:

print(assemble_scalar(form(1 * v * J * ds(vanish_voltage_surface))))

which gives 0 (I guess something is wrong?). Also:
print(v * J * ds(vanish_voltage_surface)) returns { 200000.0 * v_0[1] } * ds(<Mesh #5>[15], {}), is this normal/expected?

What is strange to me, is that if I comment out the part + v * J * ds(vanish_voltage_surface) in the weak form, I get a totally different “volt” solution (I get almost 0 along the wire), whereas with this term I get a linear voltage increment along the wire, which is the expected behavior. This is strange because this surface integral is worth 0 according to the print() statement, and so shouldn’t have any impact, yet it does… I am getting lost at comprehending what is going on.

It’s been maybe 2 years I haven’t ran the legacy fenics code, I am not even sure I can run it without installing things on my computer. But here it is, just in case…:

from fenics import *
import numpy as np
import meshio
    
the_file_name = 'file'
msh = meshio.read(the_file_name + ".msh")

meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))
meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]},
                                cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))
meshio.write("cf.xdmf", meshio.Mesh(
    points=msh.points, cells={"tetra": msh.cells["tetra"]},
    cell_data={"tetra": {"name_to_read": msh.cell_data["tetra"]["gmsh:physical"]}}))
mesh = Mesh()

with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2) 
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("cf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

dx_sd = Measure("dx",domain=mesh,subdomain_data=cf)
ds = Measure("ds", domain=mesh, subdomain_data=mf)

subdomains = cf
boundary_parts = mf
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
# Mixed function space
ME = FunctionSpace(mesh, P1*P1)
 
# Surfaces where we apply current and read voltages. Valid for configuration 1.
inj_current_surface = 5
vanish_voltage_surface = 6
V_current_contact_out = 0.0
reading_voltage_surface_0 = 3
reading_voltage_surface_1 = 4

# Define parameters from the physics of the problem.
epsi_emiss = 0.49
sigma_SB = 5.670373e-8
A_current_wire = assemble(1 * ds(inj_current_surface, domain=mesh)) # Area.
# Area of contact where current leaves
A_current_out_wire = assemble(1 * ds(vanish_voltage_surface, domain=mesh))

kappa_wire = 140.0
A_voltage_wire = assemble(1 * ds(reading_voltage_surface_0, domain=mesh))
A_voltage_wire_opp_side = assemble(1 * ds(reading_voltage_surface_1, domain=mesh))    
length_current_wire = 0.005
length_voltage_wire =  0.005
K_current_wire = kappa_wire * A_current_wire / length_current_wire
K_current_out_wire = kappa_wire * A_current_out_wire / length_voltage_wire
K_voltage_a_wire = kappa_wire * A_voltage_wire / length_voltage_wire
K_voltage_b_wire = kappa_wire * A_voltage_wire_opp_side / length_voltage_wire
K_leaks = K_current_wire + K_current_out_wire + K_voltage_a_wire + K_voltage_b_wire
rho_tensor = as_matrix([[1e-8, 0, 0], [0, 1e-8, 0], [0, 0, 1e-5]])
sigma_tensor = as_matrix([[1e8, 0, 0], [0, 1e8, 0], [0, 0, 1e5]])     
kappa_tensor = as_matrix([[140, 0, 0], [0, 140, 0], [0, 0, 28]])  
J = the_current / A_current_wire
T_ref = 300.0

bc0 = DirichletBC(ME.sub(1), V_current_contact_out, boundary_parts, vanish_voltage_surface)
bcs = [bc0]    

dTsteadyVsteady = TrialFunction(ME)
usteady, vsteady = TestFunctions(ME)            
TsteadyVsteady = Function(ME)
Tsteady, Vsteady = split(TsteadyVsteady)
     
TsteadyVsteady_init=Expression((str(T_ref),'1.0e-3'),degree=1)
TsteadyVsteady.interpolate(TsteadyVsteady_init)
J_vector = -sigma_tensor * grad(Vsteady)
Fourier_term = dot(kappa_tensor * grad(Tsteady), grad(usteady)) * dx_sd
Joule_term = dot(rho_tensor * J_vector, J_vector) * usteady * dx_sd
radiation_term = epsi_emiss * sigma_SB * (Tsteady ** 4 - T_ref**4) * usteady * (ds(1) + ds(2) + ds(3) + ds(4) + ds(5)) 
heat_leak_terms = K_current_wire * (Tsteady - T_ref) * usteady * ds(inj_current_surface) + K_current_out_wire * (Tsteady - T_ref) * usteady * ds(vanish_voltage_surface) + K_voltage_a_wire * (Tsteady - T_ref) * usteady * ds(reading_voltage_surface_0) + K_voltage_b_wire * (Tsteady - T_ref) * usteady * ds(reading_voltage_surface_1)
FsteadyT =  Fourier_term + Joule_term + heat_leak_terms + radiation_term    
FsteadyV = dot(sigma_tensor * grad(Vsteady), grad(vsteady))*dx_sd + J * vsteady * ds(inj_current_surface)
weak_form = FsteadyT + FsteadyV
Jacsteady = derivative(weak_form,TsteadyVsteady,dTsteadyVsteady)
problem = NonlinearVariationalProblem(weak_form, TsteadyVsteady, bcs, J = Jacsteady)
solver = NonlinearVariationalSolver(problem)
solver.solve()
4

My point is that if you are saying that the code is not giving you the expected answer, but you are not providing a reference as to what you expect as a result (except in writing), it is really hard for someone not familiar with the PDE/setup you are using to debug it.

Your code has many different things in it, and without a reference in legacy dolfin that produces the same result (and has the same kind of implementation), it is hard to give pointers.

This integral does not make sense, as v is a test function, and one should call dolfinx.fem.petsc.assemble_vector to assemble it. This is why I asked if you get the same result when calling
print(assemble_scalar(J*ds(vanish_voltage_surface))) in the two implementations.

Here is a slightly cleaner implementation of what you had (without using meshio):

import numpy as np
from dolfinx import log
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
                         dirichletbc, form, locate_dofs_geometrical,
                         locate_dofs_topological)
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)

mesh, ct, ft = gmshio.read_from_msh(
    "mesh.msh", MPI.COMM_WORLD, gdim=3)


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

rho = 1e-9
sigma = 1.0 / rho

kappa = 144


inj_current_surface = 14
# This corresponds to the curve the current leaves the material.
vanish_voltage_surface = 15
# Voltage value of the curve where the current leaves the material.
V_current_contact_out = 0.0
reading_voltage_surface_0 = 14
reading_voltage_surface_1 = 15

# Define the boundary conditions
left_facets = ft.find(inj_current_surface)
right_facets = ft.find(vanish_voltage_surface)
left_dofs = locate_dofs_topological(
    ME.sub(1), mesh.topology.dim-1, left_facets)
bc_volt = dirichletbc(ScalarType(V_current_contact_out), left_dofs, ME.sub(1))

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)
bc_temp_left = dirichletbc(ScalarType(300.0), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(320.0), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right, bc_volt]

x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)
the_current = 20.0  # 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))))


def temp_init(x):
    values = np.full(x.shape[1], 300, dtype=PETSc.ScalarType)
    return values


def volt_init(x):
    values = np.full(x.shape[1], 1.0e-3, 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)
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 + v * J * ds(vanish_voltage_surface)
weak_form = F_T + F_V


# Solve the PDE.

problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5
solver.report = True


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"
opts[f"{option_prefix}ksp_max_it"] = 10000
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(TempVolt)
log.set_log_level(log.LogLevel.WARNING)

assert(converged)
print(f"Number of interations: {n:d}")

lu_solver = ksp
viewer = PETSc.Viewer().createASCII("lu_output.txt")
lu_solver.view(viewer)
solver_output = open("lu_output.txt", "r")

with XDMFFile(mesh.comm, "solution/voltage.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(TempVolt.sub(1))

with XDMFFile(mesh.comm, "solution/temperature.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(TempVolt.sub(0))

I’ve also changed the solver to "lu" as you are debugging a problem, and tightened the tolerances, as the solution is of magnitude 10^-4.

1 Like
5

Thanks a lot, I’ll take the time to see your code.
Meanwhile, here are more details. The heat equation I solve is \nabla \cdot(\kappa \nabla T) + \rho |\vec J|^2=0, in a 3D geometry looking like a wire, however I can approximate it as a 1D domain (the wire is much longer than the lateral dimensions), I also apply boundary conditions at fixed z, the solution to my simple problem should not depend on x and y. So the eq. can be reduced to \kappa \frac{d^2T}{dz^2}+\rho |\vec J|^2=0 where none of the parameters depend on the temperature in my simple model. The solution is therefore a parabola. If I call T_c (“temperature cold”) the temperature at z=0 and T_h (temperature hot), the temperature at z=L, then the solution can be written as T(z)=-\frac{\rho|\vec J|^2z^2}{2\kappa}+\left( \frac{T_h-T_c}{L} + \frac{\rho|\vec J|^2L}{2\kappa}\right) z + T_c. As you can see, the solution depends on \vec J, the current density, a higher current should produce a more parabolic-like shape of the temperature distribution along z, whereas a smaller \vec J should make the solution depart little from a straight line, the temperature distribution expected if there is no current, and hence, no Joule heating.

Aside this eq. I also solve the PDE (from Maxwell equation) \nabla \cdot \vec J=0, and I make use of Ohm’s law in the wire: \vec J = (1/\rho)\vec E = - (1/\rho)\nabla V. I choose a boundary condition of V=0 at z=0, and I choose V(L) such that the current through the wire is given by the choice I put (through Neumann boundary conditions. I could have chosen another Dirichlet boundary condition instead, it would be equivalent but less convenient, as I prefer to specify the current I pass through the wire rather than apply a given voltage). For any choice of the current I input, the solution is a straight line in z. I do get a straight line with my current dolfinx code (hopefully not due to luck). So in my particular case, \vec J is homogeneous throughout the whole volume/region/wire. Thus Joule heat is also homogeneous everywhere. But my current solution with fenicsx shows no Joule heating at all.

6 
print(assemble_scalar(form(J*ds(vanish_voltage_surface))))

yields “20” with my fenicsx code, which is correct, that’s the current (in ampere) I have set to pass through the wire. So I believe the wrong part of my code is elsewhere. It’s as if the voltage equation and the temperature one were decoupled.

I am checking more closely my old (fenics) code, and I see that I had defined dTsteadyVsteady = TrialFunction(ME) which I then used to define the Jacobian Jacsteady = derivative(weak_form,TsteadyVsteady,dTsteadyVsteady), which was used to solve the PDEs: problem = NonlinearVariationalProblem(weak_form, TsteadyVsteady, bcs, J = Jacsteady). That is something I totally avoided with my new code. Could this be the culprit of having decoupled PDE’s?

Edit: I just implemented it exactly as I had done with fenicsx. Sadly this didn’t change anything, the result is still wrong and insensitive to the Joule heat term.

7

Maybe this helps for debugging:

print(assemble_scalar(form(Joule_term)))
print(assemble_scalar(form(Fourier_term)))

after solving the problem returns:

0.0009999999958870721
6.180024547663976e-11

The total Joule heat should be worth I^2R, with I= 20A and R = \rho L/A = 1e-9\Omega m \times \frac{1m}{1e-4m^2}=1e-5\Omega. This is worth 0.004 W. The code returns 0.001 W, assuming assemble_scalar(form(Joule_term)) computes Joule heat over the whole volume. So, something could be wrong. Nonetheless Joule heat is not vanishing in the code, but its effect on “temp” is null, which should not be the case.

8

I think I mentally figured out what’s wrong, but I’ll have to check tonight with my pc.

I used a very thick wire (1 cm squared cross section), and 20 amperes through it isn’t going to heat up noticeably (estimated 4 mW as above). Maybe the temperature distribution I get is fine, but depart so little from a straight line that I,didn’t notice it. I’ll change my mesh to a 1 mm squared cross section wire and see how.things go.

9

That was it. I also changed signs in the weak form, not sure if it’s equivalent to the one here. I increased the current injected, and I finally got a parabola. So the code was working fine, probably, all along, at least mostly. I just picked bad values for the parameters and didn’t notice anything.

Edit: Everything is not perfect. The computed Joule heat through the whole volume, after the solution is found, is 0.027562500000000215 W, whereas it should be 0.11025 W with my new geometry and injected current. As before (as above), the value is EXACTLY 1/4th of what I think it should be. This cannot be a coincidence… hmm, very strange. Why a factor of 4?

This is how I compute it with fenicsx: assemble_scalar(form(Joule_term)). I just checked the volume with assemble_scalar(form(1 * dx)) and it is correct. Therefore it looks like assemble_scalar(form(Joule_term)) does not return the volume integral of the Joule heat (only a fourth of it, for some reason). Not sure what’s going on.

Edit: Indeed, I miscomputed the Joule heat. It now returns the value I estimated by hand. I also get the correct voltage across the wire. So the problem is fully solved. Here’s the correct way to do the integration for the Joule heat: assemble_scalar(form(dot(rho * J_vector, J_vector) * dx))