# Newbie trying to solve a 3D moving heat source problem with crank-Nicholson for a month

I do not know what I am doing wrong. Perhaps something basic. But the Temperature is initially very high and then suddenly drops to 0.0!!

import dolfinx
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.fem import (Constant,  Function, functionspace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.mesh import create_unit_square, locate_entities, meshtags

from ufl import *
from dolfinx import default_scalar_type
from dolfinx.io import XDMFFile

import pyvista
import matplotlib as mpl
import matplotlib.pyplot as plt
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType  # type: ignore
# !pip install pyvista

#### Define Problem domain

# Parameters
cell_dim = Torch_rad # Cell dimension for the grid

# Define the corner points of the parallelepiped and the number of cells
point1 = np.array([0.0, 0.0, 0.0])

# Number of cells in each direction
num_cells = [int(point2[0] / cell_dim), int(point2[1] / cell_dim), int(point2[2] / cell_dim)]

# Create the parallelepiped mesh
domain = mesh.create_box(MPI.COMM_WORLD, [point1, point2], num_cells, cell_type=mesh.CellType.tetrahedron)
point2
array([0.01   , 0.001  , 0.00025])

#### Visualize Domain

import pyvista as pv
pv.set_jupyter_backend('static')
from dolfinx.io import XDMFFile

with XDMFFile(MPI.COMM_WORLD, "workspace/domain.xdmf", "w") as xdmf_file:
xdmf_file.write_mesh(domain)

# Load the mesh with pyvista

# Plot the mesh using PyVista
plotter = pv.Plotter()
plotter.show()
V = fem.functionspace(domain, ("CG", 1))
x = ufl.SpatialCoordinate(domain)

#### Define Material Property

# Define material properties
rho = fem.Constant(domain, PETSc.ScalarType(7870.0))  # Density
cp = fem.Constant(domain, PETSc.ScalarType(500.0))    # Specific heat capacity
k = fem.Constant(domain, PETSc.ScalarType(16.0))     # Thermal conductivity

# Define temporal parameters
t = 0.0  # Initialize time as a UFL constant
t_end = 0.008 # Final time
# dt = 0.000025
theta = 1.0

# num_steps = 100
dt = 0.0004
num_steps = int(t_end/dt)
print(f"<num_steps>")
# t_end / num_steps  # time step size

fr = 0.1
<num_steps>
CFL = fr*dt/cell_dim
print(CFL)
0.8

#### Define initial condition

# Create initial condition
def initial_condition(x):
return x[0]*0+295.0
# Define the initial condition
# ic = InitialCondition()
T_n = fem.Function(V)
T_n.name = "T_n"
T_n.interpolate(initial_condition)
T_n.vector.array.shape
(25326,)

### Boundary Condition

#### Neumann BC

class MovingHeatSource:

self.fr = fr  # Feed rate of the moving heat source
# print(self.q)
self.t = t

def __call__(self, x):
# Calculate the current position of the moving heat source
xc = point2[0]*0.1 + self.fr * t
yc = point2[1]/2 # Assume a fixed y-coordinate for simplicity
zc = point2[2]
# print('I am here')
# Calculate the distance from the heat source center
r = ((x[0] - xc)**2 + (x[1] - yc)**2)**0.5

# Set sigma based on the torch radius for rapid decay
in_z_plane = np.isclose(x[2], zc)

sigma = self.Torch_rad / 3  # Adjust sigma as needed for faster decay
heat_expr = self.q * np.exp(-r**2 / (2 * sigma**2)) * in_z_plane

return heat_expr

#### Define domain boundaries

boundaries = [(1, lambda x: np.isclose(x[2], point2[2])),# upper z surface where heat source is moving
(2, lambda x: np.isclose(x[2], point1[2])),
(3, lambda x: np.isclose(x[1], point2[1])),
(4, lambda x: np.isclose(x[1], point1[1])),
(5, lambda x: np.isclose(x[0], point2[0])),
(6, lambda x: np.isclose(x[0], point1[0]))]
facet_indices,facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(domain, 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(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
import os
notebook_dir = os.getcwd()

Debugging boundary condition

domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
with XDMFFile(domain.comm, notebook_dir +"/facet_tags.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(facet_tag, domain.geometry)
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

We can now create a general boundary condition class

T = TrialFunction(V)
v = TestFunction(V)
class BoundaryCondition():
def __init__(self, type, marker, values,dt=None):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
# self._bc = dt * inner(values, v) * ds(marker) #notice dt
self._bc = -inner(values, v) * ds(marker) #notice dt
elif type == "Robin":
# self._bc = values[0] * inner(T-values[1], v)* ds(marker)
self._bc = values[0] * inner(T*theta +(1-theta)*T_n -values[1], v)* ds(marker)

else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc

@property
def type(self):
return self._type

bcs = []

#### Setup the linear system

# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)
F = (T * v / dt + theta * k * dot(grad(T), grad(v))) * dx - (T_n * v / dt - (1 - theta) * k * inner(grad(T_n), grad(v))) * dx
# Define the Dirichlet condition
T_inf = lambda x: x[0]*0+298.0
T_infty = T_inf(x)

#### Using petsc4py to create a linear solver

solver = PETSc.KSP().create(domain.comm)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.setTolerances(rtol=1e-6, atol=1e-8)
# List to store solutions and time steps
solutions = []
solutions.append(np.copy(uh.x.array))

time_steps = []
time_steps.append(0.0)
num_steps
20

#### Solve

for i in range(num_steps):

# num_steps
print(i)

heat_source_tn = MovingHeatSource(t)

t += dt

heat_source_tnp1 = MovingHeatSource(t)

g = Function(V)
g.interpolate(lambda x: theta*heat_source_tnp1(x) + (1-theta)*heat_source_tn(x))

boundary_conditions = [BoundaryCondition("Robin", 2, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 3, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 4, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 5, (25.0, T_infty),dt=dt),
BoundaryCondition("Robin", 6, (25.0, T_infty),dt=dt),
BoundaryCondition("Neumann", 1, g, dt=dt)]

for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F+= condition.bc

# # Define the left-hand side (LHS) and right-hand side (RHS) of the equation
L = rhs(F)
a = lhs(F)

bilinear_form = fem.form(a)
A = assemble_matrix(bilinear_form, bcs=bcs) # we write this because we do not have a dirichlet BC.
A.assemble()
solver.setOperators(A)

linear_form = form(L)
b = create_vector(linear_form)

with b.localForm() as loc_b:
loc_b.set(0)

assemble_vector(b, fem.form(L))

# Apply Dirichlet boundary condition to the vector

apply_lifting(b, [bilinear_form], bcs=[bcs])

set_bc(b, bcs=bcs)

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

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

# Store the solution and the current time step
solutions.append(np.copy(uh.x.array))
time_steps.append(t)

0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# import os
# import pyvista as pv

# # Initialize PyVista
# pv.start_xvfb()

# # Assuming V and plot.vtk_mesh(V) are predefined
# grid = pv.UnstructuredGrid(*plot.vtk_mesh(V))

# # Directory to save the VTK files
# output_dir = "workspace/vtk_files_static_laser_point_refined"
# os.makedirs(output_dir, exist_ok=True)
# j=0

# # Iterate over the stored solutions and time steps for writing to VTK
# for i, (solution, t) in enumerate(zip(solutions, time_steps)):
#     if i%10==0:
#         # Update the grid with the current solution
#         print(f"Writing time step {j} to VTK file")
#         grid.point_data["uh"] = solution

#         # Define the filename
#         filename = f"{output_dir}/solution_{j:04d}.vtk"

#         # Write the grid to a VTK file
#         grid.save(filename)
#         j+=1

# print("All files have been written to the VTK format.")

sol=np.array(solutions)
sol.shape
minsol=np.zeros(21)
for i in range(10):
minsol[i]=np.min(sol[i])
minsol
array([   295.        , 100068.37581255, 107831.83703142, 116929.01594176,
135269.08819064, 152687.62496717, 164027.43308872, 168148.15419566,
171988.70673892, 180126.98577475,      0.        ,      0.        ,
0.        ,      0.        ,      0.        ,      0.        ,
0.        ,      0.        ,      0.        ,      0.        ,
0.        ])
pv.set_jupyter_backend('static')
import pyvista as pv
from dolfinx.plot import vtk_mesh

# pv.set_jupyter_backend('pythreejs')

# Convert the FEniCSx solution to a PyVista mesh
topology, cell_types, geometry = vtk_mesh(V)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)

# Add the solution to the mesh
grid.point_data["u"] = uh.x.array

# Plot the solution
plotter = pv.Plotter()