Hello!

# Background info about simulation

I’ve been trying to setup a simulation which solves the following heat equation:

with the initial condition

and the following Neumann boundary condition:

where

u(\vec{X},t) is the temperature of the cell at location \vec{X} and time t [K].

\alpha(\vec{X},u) is the thermal diffusivity of the material at location \vec{X} and temperature u [m^2/s].

f(\vec{X},t) is the increase in temperature at location \vec{X} and time t due to internal processes [K].

k(\vec{X},u) is the thermal conductivity of the material at location \vec{X} and temperature u [W/mK]

q(\vec{X},t) is the heat flux at location \vec{X} and time t [W/m^2].

for this I’ve derived following variational formulation:

# Description of issue

I am struggling to understand how to define the function q(\vec{X},t) using dolfinx given the input data I have: From a trajectory propagation and using an engineering model (Quinn 2000) I obtain the heat flux on each boundary panel of an arbitrary mesh. In other words: I have a value for q at the centroids of each boundary panel. I have not found a way to pass this data into the simulation. I have looked at creating a function space which is only defined on the centroids of the boundary panels but have not fully understood how to do this (if this is even possible). I would then create a function in this function space to which I would write my values for q.

Here is a MWE, I’ve marked the places where I would need a function space which is only defined on the centroids of the boundary panels with # TODO: (line 50 and 89)

```
import numpy as np
import ufl
from dolfinx import fem, mesh, io, plot
import pyvista
from petsc4py import PETSc
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector
import matplotlib as mpl
from mpi4py import MPI
cubeMesh = mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)
def boundary_definition_func(X):
out = np.sum(np.array([
np.isclose(X[0], 0),
np.isclose(X[0], 1),
np.isclose(X[1], 0),
np.isclose(X[1], 1),
np.isclose(X[2], 0),
np.isclose(X[2], 1),
]), axis=0) > 0
return out
boundary_definition_func.tag = 1
externalPanelRegionTags = [boundary_definition_func.tag]
externalPanelTagsList = [mesh.locate_entities(cubeMesh, cubeMesh.topology.dim-1, boundary_definition_func)]
externalPanelMarkers = []
for regionMarker, panelTags in zip(externalPanelRegionTags, externalPanelTagsList):
externalPanelMarkers.append(np.full_like(panelTags, regionMarker))
allBoundaryPanelTags = np.hstack(externalPanelTagsList).astype(np.int32)
allBoundaryPanelRegionTags = np.hstack(externalPanelMarkers).astype(np.int32)
boundaryPanelMeshTagsObj = mesh.meshtags(cubeMesh, cubeMesh.topology.dim-1, allBoundaryPanelTags, allBoundaryPanelRegionTags)
x = ufl.SpatialCoordinate(cubeMesh)
ds = ufl.Measure("ds", domain=cubeMesh, subdomain_data=boundaryPanelMeshTagsObj)
# defining function space for the temperature field (nodes at the corners of the tetrahedra)
nodeSpace = fem.FunctionSpace(cubeMesh, ("Lagrange", 1))
# defining function space for the material properties (nodes at the centre of the tetrahedra)
tetraCentreSpace = fem.FunctionSpace(cubeMesh, ("DG", 0))
# defining function space for the boundary conditions (boundary condition is the heat flux at the centroids of the boundary panels)
# TODO: define correct function space
boundarySpace = fem.FunctionSpace(cubeMesh, ("DG", 0))
# defining the solution at previous time step
u_lastStep = fem.Function(nodeSpace)
u_lastStep.name = "u_lastStep"
# defining the solution at current time step
u_h = fem.Function(nodeSpace)
u_h.name = "u_h"
# apply initial condition
u_h.x.array[:] = np.ones_like(u_h.x.array) * 273
u_lastStep.x.array[:] = np.ones_like(u_h.x.array) * 273
# defining material property field (alpha)
alpha_field = fem.Function(tetraCentreSpace)
alpha_field.name = "alpha_field"
# for this example, we will set the alpha field to a constant value (in reality a function of the material properties and temperature field)
alpha_field.x.array[:] = np.ones_like(alpha_field.x.array) * 0.00018291
# defining material property field (k)
k_field = fem.Function(tetraCentreSpace)
k_field.name = "k_field"
# for this example, we will set the k field to a constant value (in reality a function of the material properties and temperature field)
k_field.x.array[:] = np.ones_like(k_field.x.array) * 413
# defining the source term field (f)
f_field = fem.Function(tetraCentreSpace)
f_field.name = "f_field"
# for this example, we will set the source term to zero
f_field.x.array[:] = np.zeros_like(f_field.x.array)
# defining the boundary condition field (q)
q_field = fem.Function(boundarySpace)
q_field.name = "q_field"
# specifying a random heat flux at the boundary -> in reality this will be an input from a different code module (trajectory propagator)
# TODO: use correct function space for this boundary condition (only the heat flux at the centroids of the boundary panels is known)
q_field.x.array[:] = np.random.uniform(5000000, 8000000, q_field.x.array.shape[0])
# defining the time step
dt = 0.1
t_eval = 0
no_steps = 10
# setup variational problem
u, v = ufl.TrialFunction(nodeSpace), ufl.TestFunction(nodeSpace)
a = u * v * ufl.dx + dt * alpha_field * ufl.dot( ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_lastStep + dt * f_field) * v * ufl.dx
# applying neuamnn boundary condition
L -= dt * alpha_field * -1/k_field * q_field * v * ds
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form)
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(cubeMesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
# setup file to write solution to
xdmf = io.XDMFFile(cubeMesh.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(cubeMesh)
xdmf.write_function(u_h, t_eval)
pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(nodeSpace))
plotter = pyvista.Plotter()
plotter.open_gif("u_time.gif", fps=10)
grid.point_data["u_h"] = u_h.x.array
warped = grid.warp_by_scalar("u_h", factor=1)
viridis = mpl.colormaps.get_cmap("jet").resampled(50)
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, 2*max(u_h.x.array)])
for i in range(no_steps):
print(f"Time step {i} / {no_steps} completed \t\t", end="\r")
t_eval += dt
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Solve linear problem
solver.solve(b, u_h.vector)
u_h.x.scatter_forward()
# Update solution at previous time step (u_n)
u_lastStep.x.array[:] = u_h.x.array
# Write solution to file
xdmf.write_function(u_h, t_eval)
# Update plot
warped = grid.warp_by_scalar("u_h", factor=1)
plotter.update_scalars(u_h.x.array, render=False)
plotter.write_frame()
plotter.close()
xdmf.close()
```

I hope this is all clear, if not I’m happy to elaborate more