Hi all,
I would like to effectively stop my heat load at a given time ‘t_off’ during my time-dependent simulation.
I am able to see the source term, ‘gsn_heat_load.dT’ go to zero at my given time step but when I plot the temperature profile at the end of the simulation, it appears as though the system stops receiving heat at the first time step. I will provide a simplified version of the code below. Any help would be greatly appreciated. Thank you.
import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
from scipy.stats import norm
import time as tm
import copy
import warnings
from IPython.display import Latex
from fenics import*
from mshr import*
%matplotlib inline
import matplotlib.pyplot as plt
import plotly.graph_objects
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from IPython.display import Markdown, display
def printmd(string):
display(Markdown(string))
hfontLarge = {'fontname':'Latin Modern Roman', 'size' : 24, 'weight' : 'bold'}
hfontMed = {'fontname':'Latin Modern Roman', 'size' : 19, 'weight' : 'bold'}
hfontMed2 = {'fontname':'Cambria', 'size' : 19, 'weight' : 'bold'}
# geometry dimensions
leng = 2.5 # cylinder length [cm]
diam = 1.0 # cylinder diameter [cm]
# derived parameters
rad = diam / 2 # radius [cm]
lh = leng / 2 # half-length [cm]
rad2 = rad * rad # radius squared [cm^2]
# create a cylindrical geometry
# cylinder = Cylinder('coordinate of center of the top circle', 'coordinate of center of the bottom circle', 'radius of the circle at the top', 'radius of the circle at the bottom')
cylinder = Cylinder(Point(0, 0, lh), Point(0, 0, -lh), rad, rad)
geometry = cylinder
# define mesh
mden = 25 # mesh density
mesh = generate_mesh(geometry, mden)
# save the mesh
File("cylinder.pvd") << mesh
print('rad: %s cm' %(rad))
# plot mesh
print('Geometry')
print('length: %s cm, radius: %s cm' %(leng, rad))
print('Mesh')
print('mesh density: %s, no. of tetrahedral elements: %s, no. of vertices: %s' %(mden, mesh.num_cells(), mesh.num_vertices()))
mesh
# constants
T0 = 0.0 # deg C : T_0
# dT = 20.0 # deg C : T_1 - used in earlier heat load expressions
# wdT = 0.25 # laser radius [cm] (1/2 crystal radius)
wdT = 0.162 # laser radius [cm] - calculated via w_p = sqrt(2) * c_1
# ldT = 0.40 # cm
# sg_px = 8 # 1 : x-component of super gaussian exponent
# sg_py = 1 # 1 : y-component of super gaussian exponent
# rad_ratio_sq = (rad**2) / (wdT**2)
# log_rrsq = np.log(rad_ratio_sq)
# print('log of square of radius ratio: %s' %(log_rrsq))
# define function space on mesh
V = FunctionSpace(mesh, 'P', 1)
# define Dirichlet boundary condition for sides at r_max
tol = 1.e-13
def boundary_D(x, on_boundary):
return on_boundary and near(x[0]*x[0] + x[1]*x[1], rad2, tol)
bc = DirichletBC(V, Constant(T0), boundary_D)
# Gaussian (BELLA pump laser)
pix_to_deg = 6.9 # pix to micron ratio 1 pixel/6.9 micron
# a1 = 3.84e4 # Gaussian amplitude
b1 = 0.0 # 361.5
c1 = 166.0 * pix_to_deg / 1e4 # variance [cm]
print('variance, c1: %s [cm]' %(c1))
# absorption length inferred from Chenais paper: alpha = ~0.4 1/cm
# absorption length for BELLA pi-polarized pump: 1.2 1/cm (95.5% absorption)
alpha_h = 1.2 # absorption length [1/cm]
# calculate incremental temperature deposition
PumpPower = 20 # pump laser power [W]
P_abs = PumpPower * (800-532)/800 # absorbed power [W] (previously used 35 * 2 in top hat heat load)
print('Absorbed power: %f [W]' %(P_abs))
# V_abs = np.pi * (wdT/100)**2 * (leng/100) # absorption volume [m^3] (1/100 to convert cm to m)
V_eff = np.pi * c1**2 * (1 - np.exp(-alpha_h * leng)) / alpha_h # effective volume [cm^3]
K_c_tisaph = 33/100 # thermal conductivity [W/cm/K] https://www.rp-photonics.com/titanium_sapphire_lasers.html
K_c_ndyag = 0.1 # thermal conductivity of Nd:YAG [W/cm/K]
rho_tisaph = 3.98 # density of Ti:Sapphire [gm/cm^3]
rho_ndyag = 4.56 # density of Nd:YAG [gm/cm^3]
rho = rho_ndyag
c_p_tisaph = 0.756 # specific heat of Ti:Sapphire [J/gm/K]
c_p_ndyag = 0.59 # specific heat of Nd:YAG [J/gm/K]
c_p = c_p_ndyag
# dQ_incr = P_abs / V_abs # incremental heat deposition [W/m^3]
dQ_incr = P_abs / V_eff # incremental heat deposition [W/m^3]
dT_incr = dQ_incr / rho / c_p # incremental temperature deposition [K/s]
print('dQ_incr: %f [W/cm^3], dT_incr: %f [K/s]' %(dQ_incr, dT_incr))
# define heat load expressions
gsn_bella_heat_load = Expression('dT * exp( -(pow((x[0] - b1) / c1, 2) + pow((x[1] - b1) / c1, 2)) ) * exp(-alpha_h * x[2])',
degree = 1, dT = dT_incr, b1 = b1, c1 = c1, alpha_h = alpha_h)
# define initial value
u_n = interpolate(Constant(T0), V)
tol = 2e-2
u_min = u_n(rad-tol, 0., 0.)
u_max = u_n(0., 0., lh-tol)
u_n(rad-tol, 0., 0.), u_n(0., 0., lh-tol)
# define crystal properties
a_Al2O3 = 7.98e-2 # cm^2/s diffusion constant of sapphire (Al2O3): (thermal conductivity [W/cm/K] / rho [gm/cm^3] / specific heat capacity(J/gm/K))
a_perp = 7.64e-2 # perpendicular to C
a_para = 8.57e-2 # parallel to C
a_tisaph = K_c_tisaph / c_p_tisaph / rho_tisaph
a_ndyag = K_c_ndyag / c_p_ndyag / rho_ndyag
a_crystal = a_ndyag
# a_crystal = a_tisaph
# define time stepping
T = 1.5e1 # total simulation time [s]
n_steps = 1000 # number of time steps
dt = T / n_steps # size of time step [s]
nip = 20 # number of intervals between records
print('dt: %s s' %dt)
# define variational problem
u = TrialFunction(V)
v = TestFunction(V)
# source term
f = gsn_bella_heat_load # Gaussian heat load for BELLA pump laser
# f = Constant(0) # no source term
# f = uniform_heat_load # uniform cylindrical heat load
# f = long_abs_heat_load # Chenais heat load formula (2.1.6) without log term
# f = chenais_heat_load_logterm # chenais heat load eqn. 2.1.6 with log term
# f = bolus # source term bolus
# f = expsg_heat_load # source term hse
# F = u*v*dx + dt*a_Al2O3*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx + gi*v*ds # w/ Neumann BC
# F = u*v*dx + dt*a_crystal*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx # w/ Dirichlet + initial condition
F = u*v*dx + dt*a_crystal*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx # w/ Dirichlet + initial condition
a, L = lhs(F), rhs(F)
# define time-evolution function
def evolve():
# g_src = gsn_bella_heat_load
# report initial state
yield u_n
# time-stepping
u = Function(V)
t = 0
# source_on = True
# eom_set = False
for n in range(1, n_steps + 1):
# update current time
t += dt
# uncomment below to switch off pump laser at t = t_off [s]
if t_off < t:
# print('turning off source')
gsn_bella_heat_load.dT = 0.
# print('t: %s [s], gsn_bella_heat_load: %s' %(t, gsn_bella_heat_load))
# compute solution
set_log_level(30)
solve(a == L, u, bc)
# report current state
if n % nip == 0:
yield u
# update previous solution
u_n.assign(u)
# for 3D plots, get the facets, and build an array
# containing the indices of their coordinates
inds = []
for item in dolfin.cpp.mesh.facets(mesh):
inds.append(item.entities(0).tolist())
# we will provide these indices to plotly so it can draw proper surfaces
inds = np.array(inds)
ii = inds[:, 0]
jj = inds[:, 1]
kk = inds[:, 2]
# get node coördinate values and ranges
xvals = mesh.coordinates()[:,0]
yvals = mesh.coordinates()[:,1]
zvals = mesh.coordinates()[:,2]
xmin, xmax = xvals.min(), xvals.max()
ymin, ymax = yvals.min(), yvals.max()
zmin, zmax = zvals.min(), zvals.max()
n_rows = 1
n_cols = 2
# fig_wd = 15
# default sizing here yields unit aspect ratio
# plt.figure(figsize = (fig_wd, fig_wd * n_rows // n_cols))
# plt.subplot(n_rows, n_cols, idx)
fig = plt.figure(figsize=(13,8))
fig1 = fig.add_subplot(n_rows, n_cols, 1)
fig2 = fig.add_subplot(n_rows, n_cols, 2)
tol = 2.e-2 # avoid points outside domain (for a coarse mesh, increase this value) - already defined above
xv = np.linspace(xmin * (1 - tol), xmax * (1 - tol), 201)
zv = np.linspace(zmin * (1 - tol), zmax * (1 - tol), 201)
radpts = [(x_, 0, 0) for x_ in xv]
axipts = [(0, 0, z_) for z_ in zv]
# define laser shut-off time [s]
t_off = 0.25 #0.7
idx = 0 # number of time steps to complete simulation, idx*dt will yield final simulation time
T_center = []
T_entrance = []
t_vals = []
t0 = tm.time()
for u in evolve():
idx += 1
ux = np.array([u(pt) for pt in radpts])
uz = np.array([u(pt) for pt in axipts])
fig1.plot(xv, ux, lw=1, label='{0:5.3f} s'.format(idx*dt))
fig2.plot(zv, uz, lw=1)
T_center.append(u(0, 0, 0))
T_entrance.append(u(0, 0, zv[0]))
t_vals.append((idx - 1) * dt)
t1 = tm.time()
fig1.plot(xv, ux, 'k', lw=1)
fig1.legend()
fig2.plot(zv, uz, 'k', lw=1)
#fig2.legend()
fig1.set_xlabel(r'transverse position $x\,/\,\mathrm{cm}$')
fig1.set_ylabel(r'$T\,/\,{}^\circ\mathrm{C}$')
fig2.set_xlabel(r'longitudinal position $z\,/\,\mathrm{cm}$')
fig2.set_ylabel(r'$T\,/\,{}^\circ\mathrm{C}$')
fig.savefig('heat_decay_profiles.pdf')
fig.show()
print('mesh density: %s, no. of tetrahedral elements: %s, no. of vertices: %s' %(mden, mesh.num_cells(), mesh.num_vertices()))
print("simulation time: %4.3f minutes" % ((t1 - t0)/60.0))
#print("simulation time: %4.3f seconds" % (t1 - t0))
# plot center point temperature profiles at entrance and half-length
fig = plt.figure(figsize=(8,6.5))
ax = fig.gca()
ax.plot(t_vals, T_center,'r-', lw = 1, label = 'longitudinal center point')
ax.plot(t_vals, T_entrance, 'b-', lw = 1, label = 'center entrance point')
# ax.vlines(t_vals, 0, T_center, linestyle="dashed")
ax.axvline(x=t_off, color='k', ls='--')
# ax.axvline(x=0.015, color='k', ls='--')
ax.set_xlabel(r'time [s]')
ax.set_ylabel(r'$T\,/\,{}^\circ\mathrm{C}$')
print('t_off: %s seconds' %t_off)
ax.legend()
fig.show()
You can see that despite having set t_off equal to 0.25 seconds and updating gsn_heat_load to zero at this time, it appears as though the heat load ceases at the first time step of 0.015 seconds.