### Non-analytical source term in time dependend differential equation

173

views

1

I would like to solve the heat equation on a 2d mesh where the heat source (absorbed intensity) is not analytically known (except for t=0), but rather depends on the solution of the heat equation (because the heat change leads to a change in the refractive index and that will alter the intensity). Starting with the heat equation in the form

$\partial_tT\left(x,y,t\right)-K\left(\partial^2_x+\partial^2_y\right)T\left(x,y,t\right)=Q\left(x,y,t\right)$∂

where T is temperature, K is the thermal diffusivity and $Q=\alpha I$

$T^{n+1}-\Delta tK\nabla^2T^{n+1}-\Delta tQ^{n+1}-T^n=0$

In the variational form, this equation translates to

$\int_{\Omega}dx\left[uv+\Delta tK\left(\nabla\cdot u\right)\left(\nabla\cdot v\right)-\Delta tf^{n+1}v-u^nv\right]+\int_{\partial\Omega}ds\frac{Kh\Delta t}{\kappa}uv=0$∫

where I used Robin boundary condition on the edges of my mesh (see my earlier question although this is not of concern for this question and just mentioned for completeness). I also adopted the naming conventions from the tutorial.

My problem is, that $f^{n+1}=Q^{n+1}=Q\left(x,y,t_{n+1}\right)$

$\partial_tT\left(x,y,t\right)-K\left(\partial^2_x+\partial^2_y\right)T\left(x,y,t\right)=Q\left(x,y,t\right)$∂

_{t}`T`(`x`,`y`,`t`)−`K`(∂^{2}_{x}+∂^{2}_{y})`T`(`x`,`y`,`t`)=`Q`(`x`,`y`,`t`)where T is temperature, K is the thermal diffusivity and $Q=\alpha I$

`Q`=`α``I`is the source term (absorption coefficient times intensity distribution). Discretizing in time and applying an implicit Euler discretization yields$T^{n+1}-\Delta tK\nabla^2T^{n+1}-\Delta tQ^{n+1}-T^n=0$

`T`^{n+1}−Δ`t``K`∇^{2}`T`^{n+1}−Δ`t``Q`^{n+1}−`T`^{n}=0In the variational form, this equation translates to

$\int_{\Omega}dx\left[uv+\Delta tK\left(\nabla\cdot u\right)\left(\nabla\cdot v\right)-\Delta tf^{n+1}v-u^nv\right]+\int_{\partial\Omega}ds\frac{Kh\Delta t}{\kappa}uv=0$∫

_{Ω}`d``x`[`u``v`+Δ`t``K`(∇·`u`)(∇·`v`)−Δ`t``ƒ`^{n+1}`v`−`u`^{n}`v`]+∫_{∂Ω}`d``s``K``h`Δ`t``κ``u``v`=0where I used Robin boundary condition on the edges of my mesh (see my earlier question although this is not of concern for this question and just mentioned for completeness). I also adopted the naming conventions from the tutorial.

My problem is, that $f^{n+1}=Q^{n+1}=Q\left(x,y,t_{n+1}\right)$

`ƒ`^{n+1}=`Q`^{n+1}=`Q`(`x`,`y`,`t`_{n+1}) is not an analytical function but depends on the solution u. I am not sure how to deal with this in FEniCS. I came up with this solution strategy:- Construct an initial heat source and temperature $f^0$
`ƒ`^{0}and $u^0$`u`^{0} - Solve the PDE but with $f^n$
`ƒ`^{n}instead of $f^{n+1}$`ƒ`^{n+1}. The result will be $u=u^{n+1}$`u`=`u`^{n+1} - From u construct $f^{n+1}$
`ƒ`^{n+1} - Repeat 2-3 until end of time loop is reached

```
from fenics import *
import matplotlib.pyplot as plt
# Define system parameter
rho0 = 1.176685 # mass density (kg.m^-3)
c_p = 1.0049e3 # specific heat capacity at constant pressure (J.kg^-1.K^-1)
kappa = 1e-10 # thermal conductivity (W.m^-1.K^-1)
w0 = 10e-2 # beam radius
alpha = 1e-6 # absorption coefficient due to aerosols
P0 = 10e3 # beam power (W)
h = 1e-3 # heat transfer coefficient (W.m^-2.K^-1)
# Define spatial grid
w_extend = 10 # weight of the spatial extend
nx = ny = 50 # Number of grid points
# Define time steps
T = 10.0 # final time (s)
num_steps = 50 # number of time steps
dt = T/num_steps # time step size
# Derived parameter
K = kappa/rho0/c_p # thermal diffusivity
beta = alpha/rho0/c_p # prefactor in source term
# Create mesh and define function space
mesh = RectangleMesh(Point(-w_extend*w0/2, -w_extend*w0/2), Point(w_extend*w0/2, w_extend*w0/2), nx, ny, 'crossed')
V = FunctionSpace(mesh, 'P', 1)
# ============================
# Replace the expression I by something that allows me to use a numpy.array
# Define source term
I0 = 2*P0/np.pi/w0**2 # peak intensity
I = Expression('I0*exp(-(pow(x[0], 2) + pow(x[1], 2))/pow(w0, 2))',
degree=2, I0=I0, w0=w0)
# ============================
# Define initial value
u_0 = Constant(0) # ambient temperature with respect to T_0
u_n = interpolate(u_0, V) # save u_0 to u_n since this is kept during time stepping
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = beta*I
F = u*v*dx + dt*K*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx + dt*K*h/kappa*u*v*ds
a, L = lhs(F), rhs(F)
# Create VTK file for saving solution
vtkfile = File('solution/solution2d.pvd')
# Time-stepping
u = Function(V)
t = 0.0
vtkfile << (u, t)
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(a == L, u)
#=====================================
# Here, the source term f (e.g. the expression I) must be recomputed based on u
#=====================================
# Save to file and plot solution
vtkfile << (u, t)
print('t = {0:0.2f}'.format(t))
# Update previous solution
u_n.assign(u)
```

Community: FEniCS Project

### 1 Answer

3

Note that FEniCS is not restricted to solving linear problems. For instance, the demo here

https://fenics.readthedocs.io/projects/dolfin/en/stable/demos/nonlinear-poisson/python/demo_nonlinear-poisson.py.html

illustrates some of the very convenient functionality for solving nonlinear variational problems. (The ability to automatically derive Jacobians for nonlinear solves is one of my favorite features.) It should be relatively easy to adapt that example to an implicit Euler discretization where the nonlinearity is in the source term instead of the diffusivity.

https://fenics.readthedocs.io/projects/dolfin/en/stable/demos/nonlinear-poisson/python/demo_nonlinear-poisson.py.html

illustrates some of the very convenient functionality for solving nonlinear variational problems. (The ability to automatically derive Jacobians for nonlinear solves is one of my favorite features.) It should be relatively easy to adapt that example to an implicit Euler discretization where the nonlinearity is in the source term instead of the diffusivity.

Thanks for yout reply. However, I am already thinking about turning this into a 3d problem and then, your comment is not going to help. This is due to how the heat source is constructed. I basically need to propagate the intensity (electric field) in a given propagation direction (for example along the x-axis) from the first plane (perpendicular to x-axis) of the geometry to the last.

written
5 months ago by
Phillip Springer

Okay, that sounds a bit more complicated than what I was imagining!

Another thought: If you're lucky, you might be able to formulate the source term as the solution to a hyperbolic PDE with temperature-dependent coefficients, then use a mixed element to solve a coupled PDE system for the temperature and source term. However, the possibility of that would depend on the physics of the problem, and would likely need some form of stabilization for the hyperbolic subproblem. (Example: If the problem for the source term on some ray in the propagation direction can be written as an ODE with temperature-dependent coefficients (w.r.t. distance along the ray), you might be able to re-formulate the problem for the source term as a hyperbolic PDE whose characteristics are those rays.)

Another thought: If you're lucky, you might be able to formulate the source term as the solution to a hyperbolic PDE with temperature-dependent coefficients, then use a mixed element to solve a coupled PDE system for the temperature and source term. However, the possibility of that would depend on the physics of the problem, and would likely need some form of stabilization for the hyperbolic subproblem. (Example: If the problem for the source term on some ray in the propagation direction can be written as an ODE with temperature-dependent coefficients (w.r.t. distance along the ray), you might be able to re-formulate the problem for the source term as a hyperbolic PDE whose characteristics are those rays.)

written
5 months ago by
David Kamensky

I am not very familiar with this ansatz. To be honest, I would prefer, at least to get it working, to be able set the source term with something I computed not with FEniCS. That means, taking the solution as an numpy.array (via u.vector().array()) with the respective mesh.coordinates() and use it in a propagation calculation to construct Q.

written
5 months ago by
Phillip Springer

If the numpy array resulting from the non-FEniCS computation has the interpretation of coefficients of some finite element function, then you may be able to use something like

```
Q.vector().set_local(outputArray)
as_backend_type(Q.vector()).vec().ghostUpdate()
```

where Q is a Function. (The second line assumes PETSc linear algebra backend and ensures that the locally-updated values of Q's coefficients are available to other processes in parallel computations.) For more general computations involving a numpy array, you could also resort to a user-defined subclass of Expression, e.g., in my answer here:

https://www.allanswered.com/post/klzvn/#jrzer

where you can do all sorts of things in the overridden eval() method.

written
5 months ago by
David Kamensky

Please login to add an answer/comment or follow this question.