Incorrect temperature in laser heating PDE

Hi,
I am trying to solve a PDE to simulate laser-heating of a material. The heat-source (Gaussian beam of diameter d) moves with speed v, falling on top of a thin cuboidal strip, and some heat flows out because of convection.

I am able to formulate and solve the PDE, but the values of Temperature are coming ( ~12000 K) far away from my expected experimental value (~ 1800K). Can you please look into my implementation and point any errors?

Some of my anticipations:
a) Is my solver not converging to the correct solution?
b) Am I messing up the Robin/Neumann BCs somewhere?
c) Perhaps an incorrect Variational Formulation.
d) Incorrectly setting up u_old. Do I need to use project() or interpolate() as well (apologies - I don’t udnerstand when to use these)
e) On checking min(u.vector()), I obtain -64K: which is infeasible. This could point to the possible bug.

MWE:

# Define Geometry:
length_s = 0.050
width_s = 0.020
height_s = 0.010
object_dim = np.array([length_s, width_s, height_s])

# Set physical parameters
thermal_conductivity = 7.7
T_amb = 300    # temperature of ambience
convec_coeff = 50
T0 = Constant(300.)
k = Constant(7.7)  # thermal conductivity of Ti6Al4v
rho = Constant(4380)     # density
cV = Constant(570)*rho  # specific heat per unit volume 
speed = 0.0066
laser_power = 1700
laser_diameter = 0.005
init_x = -0.40*object_dim[0]/2
init_y = 0
dt = 0.01
T = (-init_x-init_x)/speed         # final time
num_steps = math.ceil(T/dt)   # number of time steps
n_probes = 5



# Create Geometry
box = Box(Point(-object_dim[0]/2, -object_dim[1]/2, -object_dim[2]/2), Point(object_dim[0]/2, object_dim[1]/2, object_dim[2]/2))
mesh = generate_mesh(box, 32)

# Define expressions

laser_flux = Expression('pow(x[0]-curr_x, 2) + pow(x[1]-curr_y, 2) < pow(diameter/2, 2) ? (-8*P/(pi*pow(diameter,2)*k))*exp(-8*(pow(x[0]-curr_x, 2) + pow(x[1]-curr_y, 2))/pow(diameter,2)) : 0',
                       degree=2,
                       diameter= laser_diameter,
                       curr_x = init_x,
                       curr_y = init_y,
                       P = laser_power,
                       k = k)

r = Expression('h*k',
               degree = 2,
               h = convec_coeff,
               k = k)

s = Expression('s',
               degree = 2,
               s = T_amb)

T_amb_exp = Expression('T_amb',
                       degree = 2,
                       T_amb = T_amb)

f = Constant(0)   #internal heat source = 0

V = FunctionSpace(mesh, 'P', 1)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
u_old = Function(V)

u_old = interpolate(T_amb_exp, V)



boundary_conditions = {0: {'Neumann': laser_flux},
                       1: {'Robin':     (r, s)},
                       2: {'Robin':     (r, s)}}

class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], object_dim[2]/2, tol)

class BoundaryBot(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], -object_dim[2]/2, tol)

class BoundarySide(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0], -object_dim[0]/2, tol) or
                                near(x[0], object_dim[0]/2, tol) or
                                near(x[1], -object_dim[1]/2, tol) or
                                near(x[1], object_dim[1]/2, tol))

# Mark boundaries
boundary_markers = MeshFunction("size_t", mesh, 2)
boundary_markers.set_all(9999)
bz1 = BoundaryTop()
bz2 = BoundaryBot()
bz3 = BoundarySide()

bz1.mark(boundary_markers, 0)
bz2.mark(boundary_markers, 1)
bz3.mark(boundary_markers, 2)

# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)


# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                          boundary_markers, i)
        bcs.append(bc)

# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']
            integrals_N.append(k*g*v*ds(i))

# Simpler Robin integrals
integrals_R = []
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        r, s = boundary_conditions[i]['Robin']
        integrals_R.append(k*r*(u - s)*v*ds(i))

F = (cV*(u-u_old)/dt*v + \
    dot(k*grad(u), grad(v)))*dx + \
    sum(integrals_R) + sum(integrals_N)

a, L = lhs(F), rhs(F)

# Compute solution
u = Function(V)

T_result = np.zeros((n_probes, num_steps))
x_probe = np.linspace(init_x, -init_x, n_probes)
curr_x = init_x
curr_y = 0


t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    curr_x += speed*dt
    laser_flux.curr_x = curr_x


    # Compute solution
    solve(a == L, u, bcs)
    T_result[:, n] = [u(xi, 0, object_dim[2]/2) for xi in x_probe]


    # Update previous solution
    u_old.assign(u)  
    print('%.2f'%(100*n/num_steps), '% completed')

# Completion of simulation

# Plot the results
t = np.linspace(0, T, num_steps)
plt.figure()
plt.plot(t,T_result[0,:])
plt.plot(t,T_result[1,:])
plt.plot(t,T_result[2,:])
plt.plot(t,T_result[3,:])
plt.plot(t,T_result[4,:])
plt.xlabel("t")
plt.ylabel("Temperature in K")
plt.legend(["$x={:.4f}$".format(xi) for xi in x_probe], ncol=2)
plt.show()

min(u.vector())
# gives value = - 64

Hi Fenics Community,

Can someone please help me around this? I have been stuck at this point for quite a long time. It would be great if someone can assist in here.

Thanks in Advance,