Hello,
I am currently working on an extension of the FEniCSx tutorial on heat equation for a 2D mesh with two subdomains.
I’m trying to verify that my code works by comparing the energy that is supplied to the domain versus the energy that is stored inside. The domain has only one Dirichlet boundary condition where energy is supplied by conduction, and all other boundaries are natural (0 flux).
I therefore expect that the cumulative energy that is supplied to the domain is always equal to the energy stored inside. For one material, this is verified but when I add a second material (2 subdomains), energy seems to be “lost” and therefore the energy supplied to the domain is higher than the energy stored.
I guess the problem comes from my variational formulation or from the way I define the different materials but after a long search for a possible bug I didn’t find anything.
Do you have any idea where the problem comes from?
Thank you very much for your suggestions.
Here is the code I use, which is very similar to the heat equation tutorial:
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh, plot
from dolfinx.io import XDMFFile
from dolfinx.fem import FunctionSpace, Function
import ufl
from ufl import Measure, FacetNormal
from petsc4py.PETSc import ScalarType
# Read the mesh
with XDMFFile(MPI.COMM_WORLD, "2Dmesh/2D_diffusion_advection_mesh.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
cell_tags = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "2Dmesh/2D_diffusion_advection_facet_mesh.xdmf", "r") as xdmf:
facet_tags = xdmf.read_meshtags(mesh, name="Grid")
# Create output file
xdmf = XDMFFile(mesh.comm, "solution2D.xdmf", "w")
xdmf.write_mesh(mesh)
# Define temporal parameters
t = 0 # Start time
T = 100000000. # Final time
num_steps = 200
dt = T / num_steps # Time step
# Define initial and imposed temperatures
Text = 20 # initial temperature in the whole domain
Tinj = 70 # temperature of Dirichlet boundary condition
# Define subdomains (2 different materials)
# for metal (assume steel)
rho_m = 8000 # density of metal in kg/m³
c_m = 420 # specific heat of metal in J/(kg*K)
cond_m = 45 # thermal conductivity of metal in W/(m·K)
alpha_m = cond_m/(rho_m*c_m)
print("Thermal diffusivity of metal:" + str(cond_m/(rho_m*c_m)))
# for water
rho_w = 997 # density of water in kg/m³
c_w = 4182 # specific heat of water in J/(kg*K)
cond_w = 0.598 # thermal conductivity of water in W/(m·K)
alpha_w = cond_w/(rho_w*c_w)
print("Thermal diffusivity of water:" + str(cond_w/(rho_w*c_w)))
# give those different properties to the domain
M = FunctionSpace(mesh, ("DG", 0))
fluid_mask = (cell_tags.values == 101)
metal_mask = (cell_tags.values == 102)
rho = Function(M)
rho.x.array[metal_mask] = np.full(metal_mask.sum(), rho_m)
rho.x.array[fluid_mask] = np.full(fluid_mask.sum(), rho_w)
c = Function(M)
c.x.array[metal_mask] = np.full(metal_mask.sum(), c_m)
c.x.array[fluid_mask] = np.full(fluid_mask.sum(), c_w)
cond = Function(M)
cond.x.array[metal_mask] = np.full(metal_mask.sum(), cond_m)
cond.x.array[fluid_mask] = np.full(fluid_mask.sum(), cond_w)
# Define integration measure
ds = Measure("ds", domain=mesh, subdomain_data=facet_tags)
n = FacetNormal(mesh)
# Finite element function space
V = FunctionSpace(mesh, ("CG", 3))
# Set initial conditions
T_n = Function(V)
T_n.name = "T_n"
T_n.x.array[:] = np.full(len(T_n.x.array), Text)
# Define time dependant output
T_h = T_n.copy()
T_h.name = "T_h"
xdmf.write_function(T_h, t)
# Trial and test functions
T, r = ufl.TrialFunction(V), ufl.TestFunction(V)
# Boundary conditions
boundary_dofs_inj = fem.locate_dofs_topological(V, mesh.topology.dim-1, facet_tags.indices[facet_tags.values == 1])
bc_inj = fem.dirichletbc(ScalarType(Tinj), boundary_dofs_inj, V)
bc_tot = [bc_inj]
# Variational problem
F = T * r * ufl.dx + dt * ufl.inner(cond/(rho*c) * ufl.grad(T), ufl.grad(r)) * ufl.dx \
- (T_n * r * ufl.dx)
a = ufl.lhs(F)
L = ufl.rhs(F)
# Preparation of linear algebra structures for time dependent problems
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = fem.petsc.assemble_matrix(bilinear_form, bcs=bc_tot)
A.assemble()
b = fem.petsc.create_vector(linear_form)
# Linear solver
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
E_stored = [0]
E_input_cond = [0]
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)
fem.petsc.assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
fem.petsc.apply_lifting(b, [bilinear_form], [bc_tot])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, bc_tot)
# 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
# Write solution to file
xdmf.write_function(T_h, t)
# Compute energy stored in the system
I = fem.form(rho * c * (T_n-Text) * ufl.dx)
E_stored.append(fem.assemble_scalar(I))
# Compute input energy from conduction (Fourier)
I2 = fem.form(dt * (ufl.dot(cond*ufl.grad(T_n), n)) * ds(1))
E_input_cond.append(fem.assemble_scalar(I2))
if i % 10 == 0:
print(i)
xdmf.close()
# Plot input energy to the system versus stored energy
plt.plot(np.array(E_stored), label = "total stored energy")
plt.plot(np.array(E_input_cond).cumsum(), label = "cumulative input energy conduction")
plt.legend(loc="upper left")
plt.xlabel('Time step')
plt.ylabel('Energy [J/m]')
plt.show()
The mesh I use is accessible from the following link: