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
Torch_rad = 0.5e-4
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])
point2 = np.array([200 * Torch_rad, 20 * Torch_rad, 5 * Torch_rad]) # Adjust these values to define the size of your parallelepiped
# 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
grid = pv.read("workspace/domain.xdmf")
# Plot the mesh using PyVista
plotter = pv.Plotter()
plotter.add_mesh(grid, show_edges=True)
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:
def __init__(self,t, Q=100.0, Torch_rad=Torch_rad, fr=fr):
self.Torch_rad = Torch_rad # Radius of the heat source
self.fr = fr # Feed rate of the moving heat source
self.q = Q/(np.pi*Torch_rad**2)
# 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)
# being an implicit we start with t+dt
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])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
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()
plotter.add_mesh(grid, show_edges=False, scalars="u", cmap="viridis")
plotter.show()
plotter.export_html('solution.html')