Mixed elements, setting the problem to get the weak form

General
1

I would like to convert a working legacy FEniCS code to FEniCSx. It’s a mixed elements problem, where I solve div(J) = 0 and the heat equation (Fourier + Joule heat terms = 0) at the same time. It is a time-independent problem. Therefore, the solutions I seek are the voltage and the temperature.

I am almost there, but I am stuck at initializing values for the voltage and temperature (otherwise Newton’s method fail to converge). Here is my full code:

Sorry, I cannot share the file.msh, too many characters for discourse website.
The mesh is generated with Gmsh, the file.geo is:

SetFactory("OpenCASCADE");
Rectangle(1) = {0, 0, 0, 1, 2, 0};
Rectangle(2) = {-1, 0, 0, 1, 2, 0};
Point(11) = {-1, 1.4, 0, 1.0};
Point(12) = {-1, 1.1, 0, 1.0};
Point(13) = {1, 0.1, 0, 1.0};
Point(14) = {1, 0.2, 0, 1.0};
Line(20) = {11, 12};
Line(21) = {14, 13};

BooleanFragments{ Surface{1}; Delete; }{ Surface{2}; Delete; }
BooleanFragments{ Curve{20}; Delete; }{Surface{2}; Delete;}
BooleanFragments{ Curve{21}; Delete; }{Surface{1}; Delete;}
Physical Surface("material A", 11) = {2};
Physical Surface("material B", 12) = {1};
Physical Curve("left_mat_A", 13) = {20};
Physical Curve("right_mat_B", 14) = {21};

My FEniCSx code shown below contains some doubts I have, as comments.

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

proc = MPI.COMM_WORLD.rank 
import meshio
def create_mesh(mesh, cell_type, prune_z=True):
    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/rectangles_smaller_contacts.msh")
   
    # Create and save one file for the mesh, and one file for the facets 
    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write("meshes/mesh_smaller_contacts.xdmf", triangle_mesh)
    meshio.write("meshes/mt_smaller_contacts.xdmf", line_mesh)

# Now read the mesh files produced above.
from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "meshes/mesh_smaller_contacts.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, "meshes/mt_smaller_contacts.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 / rho

kappa = 144

# Reminder of the mesh:
# Physical Surface("material A", 11) = {2};
# Physical Surface("material B", 12) = {1};
# Physical Curve("left_mat_A", 13) = {8};
# Physical Curve("right_mat_B", 14) = {2};
inj_current_curve = 13
vanish_voltage_curve = 14 # 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 = 13
reading_voltage_surface_1 = 14

# We define the boundary conditions
from petsc4py.PETSc import ScalarType
left_facets = ft.find(inj_current_curve)
left_dofs = locate_dofs_topological(ME.sub(1), mesh.topology.dim-1, left_facets)
bcs = [dirichletbc(ScalarType(V_current_contact_out), left_dofs, ME.sub(1))]

from dolfinx.fem import assemble_scalar, form
from ufl import Measure
x = SpatialCoordinate(mesh)  # Doubt here. Do I need this? And is it really "mesh"? proc = 0.
print(x)
#print(x.shape)
dx = Measure("dx", domain=mesh,subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)
the_current = 2.0 # Current, in amperes.
J = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
print(J)
print('Length of curve where current is injected', assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh))))
print('Length of curve where current leaves the wire', assemble_scalar(form(1 * ds(vanish_voltage_curve, domain=mesh))))

from ufl import dot
def temp_init(x):
    values = np.zeros((300.0, x.shape))
    return values

def volt_init(x):
    values = np.zeros((1.0e-3, x.shape))
    return values

# Weak form.

TempVolt.sub(0).interpolate(temp_init)
TempVolt.sub(1).interpolate(volt_init)

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_curve)
weak_form = F_T + F_V

From what I have read at How to define initial conditions in mixed formulation using dolfinx? - #3 by domi and dolfinx/demo_mixed-poisson.py at d376b20b4b7bf5b85bd3a26582ec9683864bd23e · FEniCS/dolfinx · GitHub, I need to interpolate initial values, by using “x” which is a spatial coordinate of the mesh. I am not sure I have implemented this part correctly in my code.

The error returned is:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py:343, in Function.interpolate(self, u, cells)
    341 try:
    342     # u is a Function or Expression (or pointer to one)
--> 343     _interpolate(u, cells)
    344 except TypeError:
    345     # u is callable

File /usr/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
    886     raise TypeError(f'{funcname} requires at least '
    887                     '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py:319, in Function.interpolate.<locals>._interpolate(u, cells)
    318 """Interpolate a cpp.fem.Function"""
--> 319 self._cpp_object.interpolate(u, cells)

TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfinx.cpp.fem.Function_float64, f: numpy.ndarray[numpy.float64], cells: numpy.ndarray[numpy.int32]) -> None
    2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32]) -> None
    3. (self: dolfinx.cpp.fem.Function_float64, expr: dolfinx::fem::Expression<double>, cells: numpy.ndarray[numpy.int32]) -> None

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7f883236f1f0>, <function temp_init at 0x7f88323b36d0>, array([   0,    1,    2, ..., 1773, 1774, 1775], dtype=int32)

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In [35], line 12
      8     return values
     10 # Weak form.
---> 12 TempVolt.sub(0).interpolate(temp_init)
     13 TempVolt.sub(1).interpolate(volt_init)
     15 #TempVolt_init=Expression(('300.0','1.0e-3'),degree=1)
     16 #TempVolt.interpolate(TsteadyVsteady_init)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py:348, in Function.interpolate(self, u, cells)
    346 assert callable(u)
    347 x = _cpp.fem.interpolation_coords(self._V.element, self._V.mesh, cells)
--> 348 self.interpolate(u(x), cells)

Cell In [35], line 3, in temp_init(x)
      2 def temp_init(x):
----> 3     values = np.zeros((300, x.shape))
      4     return values

TypeError: 'tuple' object cannot be interpreted as an integer
2

You need to read the documentation of what np.zeros does.
https://numpy.org/doc/stable/reference/generated/numpy.zeros.html

The input to x will be of shape (3, num_interpolation_points), thus what you want is to use
np.full(x.shape[1], 300, dtype=PETSc.ScalarType) and similarly for volt_init

1 Like
3

Thank you very much dokken, I wasn’t even aware of np.full(), nor about what “x” looks like (you do a .shape[1], but when I print() it, it just shows as “x”).

I do face the problem of Newton’s method not converging though. if I use your suggestion and add the following code:

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

The error is:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In [15], line 3
      1 from dolfinx import log
      2 log.set_log_level(log.LogLevel.WARNING)
----> 3 n, converged = solver.solve(TempVolt)
      4 assert(converged)
      5 print(f"Number of interations: {n:d}")

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py:41, in NewtonSolver.solve(self, u)
     38 def solve(self, u: fem.Function):
     39     """Solve non-linear problem into function u. Returns the number
     40     of iterations and if the solver converged."""
---> 41     n, converged = super().solve(u.vector)
     42     u.x.scatter_forward()
     43     return n, converged

RuntimeError: Newton solver did not converge because maximum number of iterations reached

How would I approach the debug of my program, at this point? I don’t think Newton’s method is supposed to diverge on the problem I think I am solving. Therefore I guess I am not solving the problem I have in mind…

Edit: Nevermind. I think I should keep a region at a fixed temperature, for example, otherwise the thermal problem will diverge because Joule heat is generated but never evacuated anywhere, so the temperature should be infinite. I guess this explains why Newton’s method diverges.

4

Dokken, I think my code works, it seems to compute voltages and temperatures fine. However, I would like to visualize it, but the code that works (for dolfinx) for a single space, does not work for mixed elements as is. Specifically the line:

u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(V)

needs to be changed. I tried to change “V” by “ME”, but then running the code returns the error RuntimeError: Can only create meshes from continuous or discontinuous Lagrange spaces. I tried some other viants, to no avail. Do you have an idea that you think would work, dokken?

Nevermind… I have found the solution, I believe, in the links I myself provided above.

Nevermind the nevermind… I get stuck at the line u_grid.point_data["volt"] = volt.x.array[dofs].real, which returns the error AttributeError: 'Indexed' object has no attribute 'x'.

At this point, my code is:

import pyvista
pyvista.set_jupyter_backend("pythreejs")
from dolfinx import plot
V0, dofs = ME.sub(1).collapse()

u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(V0)

u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["volt"] = volt.x.array[dofs].real
u_grid.set_active_scalars("volt")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()

Yielding:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [31], line 37
     17 '''
     18 if have_pyvista:
     19     # Create a VTK 'mesh' with 'nodes' at the function dofs
   (...)
     33 u0.x.array[:] = u.x.array
     34 '''
     36 u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
---> 37 u_grid.point_data["volt"] = volt.x.array[dofs].real
     38 u_grid.set_active_scalars("volt")
     39 u_plotter = pyvista.Plotter()

AttributeError: 'Indexed' object has no attribute 'x'

I haven’t found a remedy so far.

5

and you are trying to call

While you should call
TempVolt.x.array[dofs].real

For the sake of transparity, it would be better if you post a new post when you resolve any errors, instead of editing your previous posts.

1 Like
6

Thanks again a lot dokken. OK about the edits, I will proceed like this in the future. I can finally visualize both the voltage and the temperature (voltage looks OK but the temperature is a little strange, I will investigate if I didn’t forget something in the weak form).

Probably lastly, I wish to compute and plot the electric field. My code doesn’t return any error, but I don’t see any vector plotted. In other words, I only see the mesh being displayed, without any vector at any point, whereas the same code yields the correct electric field in the same problem without mixed elements (without the temperature part).

# 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(volt), W.element.interpolation_points())
E_field.interpolate(E_field_expr)

from dolfinx.mesh import compute_midpoints
plotter = pyvista.Plotter()
plotter.set_position([0,0,10])

# We include ghosts cells as we access all degrees of freedom (including ghosts) on each process
top_imap = mesh.topology.index_map(mesh.topology.dim)
num_cells = top_imap.size_local + top_imap.num_ghosts
midpoints = compute_midpoints(mesh, mesh.topology.dim, range(num_cells))

num_dofs = W.dofmap.index_map.size_local +  W.dofmap.index_map.num_ghosts
assert(num_cells == num_dofs)
values = np.zeros((num_dofs, 3), dtype=np.float64)
values[:, :mesh.geometry.dim] = E_field.x.array.real.reshape(num_dofs, W.dofmap.index_map_bs)
cloud = pyvista.PolyData(midpoints)
cloud["E_field"] = values
glyphs = cloud.glyph("E_field", factor=0.4)
actor = plotter.add_mesh(grid, style="wireframe", color="k")
actor2 = plotter.add_mesh(glyphs)

if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    B_fig = plotter.screenshot("E_field.png")

I don’t see what could be wrong here.

7

Start by printing the min and Max values of this array.

Then try changing the scaling factor here.

I guess you have taken inspiration from Electromagnetics example — FEniCSx tutorial

1 Like
8

Thank you very much @dokken for all the pointers, tips and so on. This was very helpful. I had to increase the scaling factor by 1e9 (I had blindly tried 1e6 but still saw nothing. I did not expect it to differ from my other code, so I still have things to figure out).