Shut off source term at given time step

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

image

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.

Please format your code with markdown syntax using 3x backtick (`) encapsulation, such that it renders as

```python
def test(x):
    return x
```

Please also include a code that can reproduce the issue. The code you have supplied is missing all core definitions of meshes, function spaces etc, which makes it really challenging for people to help you.

Hi dokken,

I have updated the code in my post. Please let me know if that works - I tried to trim it down as much as possible.

Thanks,
Nick

You can verify that the expression is properly updated by adding:

        print(t, t_off, assemble(gsn_bella_heat_load*dx(domain=mesh)))

inside your for loop (after updating gsn).
Note that by not using yield, but plain for loops, the code runs as expected:

# define laser shut-off time [s]
t_off = 0.7     #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()


# 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))
    print(t, t_off, assemble(dt*gsn_bella_heat_load*dx(domain=mesh)), assemble(L).norm("l2"))
    # compute solution
    set_log_level(30)
    solve(a == L, u, bc)



    # update previous solution
    u_n.assign(u)



    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)

I suggest you either skip the usage of yield or read up on how it works.

Your original code is working correctly, and the heat load is shutting off at the correct time. However, in your for loop, you are not calculating the current time correctly:

The evolve function only yields a value once every 20 timesteps (because nip = 20), so instead of idx += 1, you should use idx += nip. Alternatively, you could rewrite your evolve function as follows:

def evolve():

    # g_src = gsn_bella_heat_load
    
    # time-stepping
    u = Function(V)
    t = 0
    
    # report initial state
    yield u_n, t

    # 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, assemble(gsn_bella_heat_load*dx(domain=mesh))) )
        
        # compute solution
        set_log_level(30)
        solve(a == L, u, bc)

        # report current state
        if n % nip == 0:
            yield u, t

        # update previous solution
        u_n.assign(u)

and change your for loop to match:

for u, curr_t 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(curr_t)
1 Like

Hi @dokken and @conpierce8,

Your solutions worked well for me - thank you very much for the help!

Best,
Nick