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

**URL:** <https://fenicsproject.discourse.group/t/newbie-trying-to-solve-a-3d-moving-heat-source-problem-with-crank-nicholson-for-a-month/14990>\
**Category:** General\
**Tags:** dolfinx\
**Created:** [June 14, 2024, 2:17pm UTC](https://fenicsproject.discourse.group/t/newbie-trying-to-solve-a-3d-moving-heat-source-problem-with-crank-nicholson-for-a-month/14990 "2024-06-14T14:17:02Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![suparnob](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/suparnob/32/7130_2.png) [@suparnob](https://fenicsproject.discourse.group/u/suparnob)\
**Post date:** [June 14, 2024, 2:17pm UTC](https://fenicsproject.discourse.group/t/newbie-trying-to-solve-a-3d-moving-heat-source-problem-with-crank-nicholson-for-a-month/14990/1 "2024-06-14T14:17:02Z")

</div>

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!!

```python
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

```

```python
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType # type: ignore

```

```python
# !pip install pyvista

```

#### Define Problem domain

```python
# 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)

```

```python
point2

```

```
array([0.01 , 0.001 , 0.00025])

```

#### Visualize Domain

```python
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()

```

```python
V = fem.functionspace(domain, ("CG", 1))
x = ufl.SpatialCoordinate(domain)

```

#### Define Material Property

```python
# 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>

```

```python
CFL = fr*dt/cell_dim
print(CFL)

```

```
0.8

```

#### Define initial condition

```python
# Create initial condition
def initial_condition(x):
    return x[0]*0+295.0

```

```python
# Define the initial condition
# ic = InitialCondition()
T_n = fem.Function(V)
T_n.name = "T_n"
T_n.interpolate(initial_condition)

```

```python
T_n.vector.array.shape

```

```
(25326,)

```

### Boundary Condition

#### Neumann BC

```python
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

```python
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]))]

```

```python
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])

```

```python
import os
notebook_dir = os.getcwd()

```

Debugging boundary condition

```python
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)

```

```python
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

```

We can now create a general boundary condition class

```python
T = TrialFunction(V)
v = TestFunction(V)

```

```python
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

```python
# Define solution variable, and interpolate initial solution for visualization in Paraview
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)

```

```python
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

```

```python
# 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

```python
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)

```

```python
# List to store solutions and time steps
solutions = []
solutions.append(np.copy(uh.x.array))

time_steps = []
time_steps.append(0.0)

```

```python
num_steps

```

```
20

```

#### Solve

```python
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

```

```python
# 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.")

```

```python
sol=np.array(solutions)
sol.shape

```

```python
minsol=np.zeros(21)
for i in range(10):
    minsol[i]=np.min(sol[i])

```

```python
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. ])

```

```python
pv.set_jupyter_backend('static')

```

```python
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')

```

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [June 14, 2024, 7:00pm UTC](https://fenicsproject.discourse.group/t/newbie-trying-to-solve-a-3d-moving-heat-source-problem-with-crank-nicholson-for-a-month/14990/2 "2024-06-14T19:00:20Z")

</div>

Unfortunately your question and code example is extremely long. Consider: [Read before posting: How do I get my question answered?](https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21)

This forum is populated by industry workers, academics and students who volunteer their help as a contribution to the community. By posting a much more terse question with only the salient points, you’re much more likely to receive a response.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [June 16, 2024, 5:04am UTC](https://fenicsproject.discourse.group/t/newbie-trying-to-solve-a-3d-moving-heat-source-problem-with-crank-nicholson-for-a-month/14990/3 "2024-06-16T05:04:56Z")

</div>

Revised in [3D moving heat source problem with crank-Nicholson ODE solver algorithm](https://fenicsproject.discourse.group/t/3d-moving-heat-source-problem-with-crank-nicholson-ode-solver-algorithm/15001).

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [June 16, 2024, 5:05am UTC](https://fenicsproject.discourse.group/t/newbie-trying-to-solve-a-3d-moving-heat-source-problem-with-crank-nicholson-for-a-month/14990/4 "2024-06-16T05:05:02Z")

</div>