Hi everyone,
I have written a code which solves for a transient heat conduction equation with neumann and robin boundary conditions. I am trying to save the simulation to a gif file but the kernel dies automatically, every time I run the code. Can someone tell me what’s wrong? I guess the problem is in the last 5 lines of the code when I am updating the gif. The code runs without an issue when I comment out those lines
I am running the code on Windows installation via docker!
Also, there are no dirichlet boundary conditions for my problem. When I apply lifting on the b vector, I have written it as apply_lifting(b, [a], [[]]). Is that correct?
from dolfinx import default_scalar_type
import dolfinx
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, form, locate_dofs_topological)
from dolfinx.fem.petsc import (LinearProblem, assemble_vector, assemble_matrix, create_vector,
apply_lifting, set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, locate_entities, locate_entities_boundary, meshtags, CellType
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, div, dot, dx, grad,
inner, lhs, rhs)
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pyvista
nx, ny = 50, 50
mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]),np.array([5,5])], [nx, ny], CellType.triangle)
V = FunctionSpace(mesh, ("Lagrange",1))
t = 0 # Start time
t_end = 20 # End time
num_steps = 50 # Number of time steps
dt = (t_end-t)/num_steps # Time step size
x = SpatialCoordinate(mesh)
T_n = Function(V)
T_n.interpolate(lambda x: np.full((x.shape[1],), 323))
T_h = Function(V)
T_h.interpolate(lambda x: np.full((x.shape[1],), 323))
T_inf = Constant(mesh,default_scalar_type(298)) # Ambient temperature
Q_gen = Constant(mesh,default_scalar_type(100)) # Heat generation inside CV
q_tabs = Constant(mesh,default_scalar_type(10)) # Heat flux to tabs
rho = Constant(mesh,default_scalar_type(2500)) # Density
C_p = Constant(mesh,default_scalar_type(900)) # Specific heat capacity of the cell
k = Constant(mesh,default_scalar_type(25)) # Thermal conductivity
h = Constant(mesh,default_scalar_type(8)) # Convection coefficient
boundaries = [(1, lambda x: np.isclose(x[0],0)),
(2, lambda x: np.isclose(x[1],5)),
(3, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],4),np.less_equal(x[1],5))),
(4, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],3),np.less(x[1],4))),
(5, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],2),np.less(x[1],3))),
(6, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],1),np.less(x[1],2))),
(7, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],0),np.less(x[1],1))),
(8, lambda x: np.isclose(x[1],0))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
u, v = TrialFunction(V), TestFunction(V)
a = form(rho * C_p * u * v * dx + dt * inner(k * grad(u), grad(v)) * dx +
dt * h * u * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))
L = form((rho * C_p * T_n + dt * Q_gen) * v * dx - dt * q_tabs * v * (ds(4) + ds(6)) +
dt * h * T_inf * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))
A = assemble_matrix(a, bcs=[])
A.assemble()
b = create_vector(L)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
plotter = pyvista.Plotter()
plotter.open_gif("heat_equation.gif")
grid.point_data["T_h"] = T_h.x.array
# warped = grid.warp_by_scalar("T_h", factor=1)
grid.set_active_scalars("T_h")
viridis = plt.cm.get_cmap("viridis", 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(grid, show_edges=True, lighting=False,
cmap=viridis, scalar_bar_args=sargs,
clim=[0, max(T_h.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, L)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [a], [[]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [])
# Solve linear problem
solver.solve(b, T_h.vector)
T_h.x.scatter_forward()
# Update solution at previous time step (u_n)
T_n.x.array[:] = T_h.x.array
# Update plot
# warped = grid.warp_by_scalar("T_h", factor=1)
grid.set_active_scalars("T_h")
plotter.update_coordinates(grid.points.copy(), render=False)
plotter.update_scalars(T_h.x.array, render=False)
plotter.write_frame()
plotter.close()