Doubt about heat problem

General
1

Hello there! I’m studying the heat equation from this website: The heat equation — FEniCSx tutorial.

import matplotlib.pyplot as plt
import pyvista
import ufl
import numpy as np

from petsc4py import PETSc
from mpi4py import MPI

from dolfinx import fem, mesh, io, plot
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc

# Define temporal parameters
t = 0  # Start time
T = 1.0  # Final time
num_steps = 50
dt = T / num_steps  # time step size

# Define mesh
nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-2, -2]), np.array([2, 2])],
                               [nx, ny], mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("Lagrange", 1))

# Create initial condition
def initial_condition(x, a=5):
    return np.exp(-a * (x[0]**2 + x[1]**2))


u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

# Create boundary condition
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)

xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(domain)

# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
xdmf.write_function(uh, t)

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))
a = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx

bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

import matplotlib as mpl
pyvista.start_xvfb()

grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(V))

plotter = pyvista.Plotter()
plotter.open_gif("u_time.gif", fps=10)

grid.point_data["uh"] = uh.x.array
warped = grid.warp_by_scalar("uh", factor=1)

viridis = mpl.colormaps.get_cmap("viridis").resampled(25)
sargs = dict(title_font_size=25, label_font_size=20, fmt="%.2e", color="black",
             position_x=0.1, position_y=0.8, width=0.8, height=0.1)

renderer = plotter.add_mesh(warped, show_edges=True, lighting=False,
                            cmap=viridis, scalar_bar_args=sargs,
                            clim=[0, max(uh.x.array)])

for i in range(num_steps):
    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, linear_form)

    # Apply Dirichlet boundary condition to the vector
    apply_lifting(b, [bilinear_form], [[bc]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, [bc])

    # Solve linear problem
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array

    # Write solution to file
    xdmf.write_function(uh, t)
    # Update plot
    warped = grid.warp_by_scalar("uh", factor=1)
    plotter.update_coordinates(warped.points.copy(), render=False)
    plotter.update_scalars(uh.x.array, render=False)
    plotter.write_frame()
plotter.close()
xdmf.close()

But this error:

HDF5-DIAG: Error detected in HDF5 (1.10.7) MPI-process 0:
  #000: ../../../src/H5F.c line 366 in H5Fcreate(): unable to create file
    major: File accessibility
    minor: Unable to open file
  #001: ../../../src/H5Fint.c line 1713 in H5F_open(): unable to lock the file
    major: File accessibility
    minor: Unable to lock file
  #002: ../../../src/H5FD.c line 1675 in H5FD_lock(): driver lock request failed
    major: Virtual File Layer
    minor: Unable to lock file
  #003: ../../../src/H5FDsec2.c line 990 in H5FD__sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
    major: Virtual File Layer
    minor: Unable to lock file
Traceback (most recent call last):
  File "/home/jandui/Documents/MATH code/AAAAAAAAAAAAa/heat_updated.py", line 40, in <module>
    xdmf = io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
RuntimeError: Failed to create HDF5 file.

My questions are:

Why is this error about xmdf appearing? I don’t know what I should install now; most of the libraries come with fenicsx.

Also, how can I create multiple heat points with specific shapes at the square domain with different functions and constants? (I mean, maybe put two heat sourcers (one at center in a circular domain and another one in the left corner with sin or exponential function, just an example)

Thanks for your attention.

2

This usually means that you already have a file named diffusion.h5, that has been opened by another process and has not been closed. Please check if this is the case, and remove the file.

Heat sources are covered in

3

Ok, thanks, I’ll read it now. Just one little doubt please:

  File "/home/jandui/Documents/MATH code/testes_2/Heat.py", line 76, in <module>
    viridis = mpl.colormaps.get_cmap("viridis").resampled(25)
AttributeError: 'ColormapRegistry' object has no attribute 'get_cmap'

Should I replace ‘get_cmap’ with something else? I’m not very familiar with how to plot the solutions. I’m using matplotlib 3.5 to avoid conflicts with the Paraview version.

4

This is a pure matplotlib error. Please see related questions such as: AttributeError: 'ColormapRegistry' object has no attribute 'get_cmap' · Issue #74 · dexplo/dataframe_image · GitHub
Which can be found by googling

5

Ok, I changed following what the sites said, but is not working yet:

import matplotlib.pyplot as plt

viridis = plt.cm.get_cmap('Reds')

sargs = dict(title_font_size=25, label_font_size=20, fmt="%.2e", color="black",
             position_x=0.1, position_y=0.8, width=0.8, height=0.1)

renderer = plotter.add_mesh(warped, show_edges=True, lighting=False,
                            cmap=viridis, scalar_bar_args=sargs,
                            clim=[0, max(uh.x.array)])

And this message appears (not an error message)

  warnings.warn(

WARNING:py.warnings:/home/jandui/.local/lib/python3.10/site-packages/pyvista/plotting/plotter.py:4633: PyVistaDeprecationWarning: This method is deprecated and will be removed in a future version of PyVista. Directly modify the scalars of a mesh in-place instead.
  warnings.warn(

And I thought that it could be running the program, and I waited for an hour, but I didn’t get the results (my PC is definitely not that weak). What do you think could be the issue?

6

I would suggest removing the pyvista plotting and write directly to file (“xdmf” or “bp”) and visualize it that way. It is hard for me to debug error related directly to pyvista.

7

I solved this issue by

pip install scanpy
1 Like