# Doubt about heat problem

**URL:** <https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586>\
**Category:** General\
**Created:** [February 11, 2024, 7:33pm UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586 "2024-02-11T19:33:44Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![janduigamer13](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/janduigamer13/32/6547_2.png) [@janduigamer13](https://fenicsproject.discourse.group/u/janduigamer13)\
**Post date:** [February 11, 2024, 7:33pm UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586/1 "2024-02-11T19:33:44Z")

</div>

Hello there! I’m studying the heat equation from this website: [The heat equation — FEniCSx tutorial](https://jsdokken.com/dolfinx-tutorial/chapter2/heat_equation.html).

```auto
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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 11, 2024, 7:41pm UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586/2 "2024-02-11T19:41:49Z")

</div>

> [@janduigamer13](#):
>
> ```auto
> 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.
> 
> ```

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.

> [@janduigamer13](#):
>
> 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)

Heat sources are covered in

> [@Point Sources (Redux)](https://fenicsproject.discourse.group/t/point-sources-redux/13496/6):
>
> Following my post above, with a minor adaptations, a point source can be implemented as follows: # Create a point source for Poisson problem # Author: Jørgen S. Dokken # SPDX-License-Identifier: MIT from mpi4py import MPI from petsc4py import PETSc import dolfinx import dolfinx.fem.petsc import numpy as np import ufl def compute\_cell\_contributions(V, points): # Determine what process owns a point and what cells it lies within mesh = V.mesh \_, \_, owning\_points, cells = dolfinx.cp…

---

<div class="post-metadata">

**Author:** ![janduigamer13](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/janduigamer13/32/6547_2.png) [@janduigamer13](https://fenicsproject.discourse.group/u/janduigamer13)\
**Post date:** [February 11, 2024, 7:51pm UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586/3 "2024-02-11T19:51:20Z")

</div>

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

```auto
  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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 11, 2024, 8:28pm UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586/4 "2024-02-11T20:28:27Z")

</div>

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](https://github.com/dexplo/dataframe_image/issues/74)  
Which can be found by googling

> [@janduigamer13](#):
>
> `AttributeError: 'ColormapRegistry' object has no attribute 'get_cmap'`

> <https://stackoverflow.com/questions/52251582/matplotlib-has-no-attribute-cm-when-deploying-an-app>

---

<div class="post-metadata">

**Author:** ![janduigamer13](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/janduigamer13/32/6547_2.png) [@janduigamer13](https://fenicsproject.discourse.group/u/janduigamer13)\
**Post date:** [February 13, 2024, 10:23am UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586/5 "2024-02-13T10:23:59Z")

</div>

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

```auto
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)

```auto
  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?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 13, 2024, 3:04pm UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586/6 "2024-02-13T15:04:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Yanjun](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/yanjun/32/6462_2.png) [@Yanjun](https://fenicsproject.discourse.group/u/Yanjun)\
**Post date:** [March 18, 2024, 10:21pm UTC](https://fenicsproject.discourse.group/t/doubt-about-heat-problem/13586/7 "2024-03-18T22:21:06Z")

</div>

I solved this issue by

```auto
pip install scanpy

```