Localise radiation source on upper face of cylindrical mesh

Hi everyone,

I am currently trying to model the heating of a cylindrical rod via a Gaussian laser on its upper face. I am also accounting for radiative losses during heating. The equation of the laser is:

P(x,y,z) = \frac{2AP_{max}}{m}\exp{-\frac{r^2}{w^2}}\exp{-\frac{z^2}{w_2^2}}

Where A is the absorptivity of the rod, P_{Max} is the peak power of the laser, m is the mass of the rod, w is the width of the Gaussian beam and w_2 is the width of said beam along the z direction (I have set this to be very small such that the beam is incident on the top surface, but I am unsure if this is correct).

Whilst this seemingly works well, my simulation end up not replicating the experiments I have performed, wherein I found a significant temperature difference (of around 300K) between the top and bottom of the rod, but my result indicates little temperature difference. This leads me to suspect that I have done something wrong with localizing the laser to just illuminate the top surface.

Is there a easy way to be able to insert the Gaussian beam to just hit the top of the rod, rather than having to localise it via the w_2 I used in the equation above.
In addition, is there a way in which i can extract the temperatures at the center of the top and bottom surfaces so I plot the temperature difference as a function of time?

My full code is shown below:

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
from ufl import as_tensor
from ufl import Index
import math
import mshr

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True


t_start = 0.0
t_end = 200
nstep = 200
dtime = (t_end - t_start)/ nstep  #time step
#Configure for dimensions of system
#----------THERMAL--PROPERTIES--[Element: W]--------
kappa = 0.0172         #conductivity [W/mm K] 
c = 0.132         #Specific Heat Capacity [J/gK]
rho = 0.01925         #Density [g/mm^3]
const = kappa /(c * rho)
tau_T = 0.0 #Temperature Gradient Lag
tau_q = 0.0 #Heat Flux Lag

pi = 3.141592653589793
T_am = 298 #ambient vacuum temperature (K)
T_a4 = T_am**4
epsilon = 0.1  # material emissivity
sigma = 5.67E-14 # W/(mm**2.K**4)
es = epsilon*sigma
Pmax = 400 #peak power of laser in W
w = 0.5
R = 1.5     #mm
area = pi*R*R  #area of target
depth = 8.0 #thickness of target
volume = area * depth #define a disk through which the laser is absorbed
Pden = Pmax / (volume* rho) #need to deine the power per g to balance units
absorb_depth = 1e-6

Laser = Expression('2*A*Pden*exp(-pow((x[0] - 0), 2)/(w*w)-pow((x[1]-0), 2)/(w*w)-pow((x[2]-0), 2)/(w2*w2) )',degree=3, A=epsilon, Pden=Pden,w= 0.5, w2= absorb_depth,z1=z1) #power (w2 localises the z-coordinates)
#A =! emissivity for system in thermal equillibrium

geometry = mshr.Cylinder(Point(0,0,0),Point(0,0,-8),R,R) #8mm long target 
# Create mesh
mesh = generate_mesh(geometry, 40)          # generate a mesh from the given geometry
Space = FunctionSpace((mesh), 'P', 1)      #define finite element function space, defined via basis functions
VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
cells = MeshFunction('size_t',mesh,mesh.topology().dim()) #codimension of 0
facets = MeshFunction('size_t',mesh,mesh.topology().dim()-1) #codimension of 1
da = Measure('ds', domain=mesh, subdomain_data = facets)  #area element
dv = Measure('dx', domain=mesh, subdomain_data = cells)   #volume element

initial_T = Expression("Tini", Tini=T_am, degree=3) # extrapolate an expression for the temperature before heating
T0 = interpolate(initial_T, Space)
T = Function(Space)         # Temperature
V = TestFunction(Space)     # Test Function used for FEA
dT = TrialFunction(Space)   # Temperature Derivative
q0 = Function(VectorSpace)  
i = Index()
G = as_tensor(T.dx(i), (i))  #gradient of T
G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step 

q = as_tensor(dtime/(dtime + tau_q) * (tau_q/dtime*q0[i] - kappa*(1+tau_T/dtime)*G[i] + kappa*tau_T/dtime*G0[i]),(i)) #heat
F = (rho*c/dtime*(T-T0)*V - q[i]*V.dx(i) - rho*Laser*V ) * dv + es*(T**4 - T_a4)*V*da   #final form to solve
Gain = derivative(F, T, dT)    # Gain will be usedas the Jacobian required to determine the evolution of a dynamic system 

file_T = File('target3D/solution.pvd')
for t in np.arange(t_start,t_end,dtime):
	print( "Time", t)
	solve(F==0, T, [], J = Gain, solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, form_compiler_parameters={"cpp_optimize": True, "representation": "quadrature","quadrature_degree": 2} )
	file_T << (T,t)
	q0 = project(q, VectorSpace)
	T0.assign(T)      #change so that T0 is equal to T for the next time step

You cannot do this, as you are now creating a new variable q0, not assigning values to the previously defined function as you want to do. This should be done in the same manner as for T0.
To plot such values, I recommend using Paraview, as it has features as “Plot over line”.

I have struggled with this before, but I have found that when I replace

q0 = project(q, VectorSpace)



I end up with an AssertionError

raceback (most recent call last):
  File "target3D.py", line 84, in <module>
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 414, in assign
    linear_comb = _check_and_contract_linear_comb(rhs, self, multi_index)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 154, in _check_and_contract_linear_comb
    linear_comb = _check_and_extract_functions(expr, multi_index=multi_index)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 127, in _check_and_extract_functions
    linear_comb = _check_mul_and_division(e, linear_comb, scalar_weight, multi_index)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py", line 81, in _check_mul_and_division
    assert isinstance(expr, Indexed)

Then do it in two steps,

q_tmp = project(q, VectorSpace)

That fixes the problem. Thank you very much for your help

Hi everyone:

I am currently trying to replicate the use of the radiative boundary condition \frac{\partial T}{\partial n} = \varepsilon \sigma_{r} ( T^{4} - T_{ext}^{4} , with dolfinx, but as there is no solve function I am struggling a bit …

if this boundary condition is applied on a part \Gamma of the boundary \partial \Omega, I believe a term \int_{\Gamma} \varepsilon \sigma_{r} T^{4} \varphi dS should be added to the weak formulation.

( I know that there is a constant term but it is less problematic)

I tried

V = FunctionSpace(mesh, ("Lagrange", 1))
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(mesh, PETSc.ScalarType(0))
# Robin Radiation
Radiation_LHS = inner(h_rad*u**4,v)*ds(1)

# Weak form
# Add various terms to the weak formulation
# step by step 
a = RhoCp*u*v*dx 
a += dt*dot(k*grad(u), grad(v))*dx
a += dt*Radiation_LHS

L = (RhoCp*Tinit + dt * f) * v * dx

F = a - L
problem = NonlinearProblem(F, u, bcs=[bc])

with obvious notations, but it yield the error

Traceback (most recent call last):
  File "/mnt/c/Users/quent/Documents/Codes/FENICSX/CaseSetup/HeatTransfer/BarConduction_RobinTypeBC/ExtMesh/Run.py", line 153, in <module>
    problem = NonlinearProblem(F, u, bcs=[bc])
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 712, in __init__
    self._L = _create_form(F, form_compiler_options=form_compiler_options,
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 188, in form
    return _create_form(form)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 183, in _create_form
    return _form(form)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 141, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 88, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, options) for form in forms)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 88, in <genexpr>
    form_data = tuple(_analyze_form(form, options) for form in forms)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ffcx/analysis.py", line 163, in _analyze_form
    form_data = ufl.algorithms.compute_form_data(
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py", line 415, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py", line 198, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py", line 179, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 34, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 100, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/home/quentin_digimind/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py", line 53, in sum
    raise ArityMismatch(f"Adding expressions with non-matching form arguments {_afmt(a)} vs {_afmt(b)}.")
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ('v_0', 'v_1').

I did not added the bits related to the mesh but it is generated with the classical gmsh workflow to gte tags and works fine when used on a "classical " heat equation, meaning with this boundary condition

Thanks by advance,


See: Invalid coefficient type for derivative - #2 by dokken

1 Like

Works by defining

u = Function(V)

instead of

u = TrialFunction(V)

may I ask what is the difference between these two definitions ?



is a DOLFINx object, that can the values of degrees of freedom, accessed through for instance u.x.array.
It has known values, and are supported over the whole domain.

Indiciates that this is an unknown value, which has no coefficient.
When paired with a TestFunction(V) in a variational form, it instructs DOLFINx to assemble a bi-linear form into a matrix, see for instance:

Thanks for the clarification.