Source term ineffective
My goal is to model heat diffusion at the gold / water interface during and after perturbation by an ultrafast laser pulse. I am starting with Fenics and I have been stuck on how to properly use a source term for the incoming pulse. All I get setting time 't' to whatever value is a flat green square:
I would expect to see a pulse appearing near the edge... If I manage to get this cleared up, I would move on to define a time-stepping loop and create a subdomain for water.
I look forward to your comments and suggestions!
Here's a minimum (non)working example:
from dolfin import * from matplotlib import pyplot import numpy as np parameters["plotting_backend"] = "matplotlib" t_i = 0. # initial time / s t_f = 2e-12 # final time / s t_max = 1e-12 # max of the pulse / s num_steps = 1000 # number of time steps dt = (t_f - t_i) / num_steps # time step size initial_temperature = 0. # K C_metal = .129 # heat capacity, gold / J g-1 K-1 K_metal = 3.18 # thermal conductivity, gold / W cm-1 K-1 rho_metal = 19.3 # density, gold / g cm-3 c = K_metal / C_metal / rho_metal tau = 0.1e-12 # pulse width / s w = 0.1 # pulse diameter / cm alpha = 8.4197e5 # absorption coefficient / cm-1 # Create mesh and define function space nx = 32 ny = 32 width = 1e-6 height = 1e-6 mesh = RectangleMesh(Point(-width/2., -height/2.), Point(width/2., height/2.), nx, ny) V = FunctionSpace(mesh, 'P', 1) # Define boundary condition u_D = Constant(initial_temperature) class Boundary(SubDomain): def inside(self, x, on_boundary): return on_boundary boundary = Boundary() bc = DirichletBC(V, Constant(initial_temperature), boundary) # Initial condition u_n = interpolate(u_D, V) # Define variational problem u = TrialFunction(V) v = TestFunction(V) # Source term f = Expression('exp(-pow((x-bw)/w, 2)) * exp(-pow(((t-t0)/tau), 2)) * exp(-alpha*x)', degree=2, w=w, tau=tau/2., alpha=alpha, t=0., t0=t_max, bw=width/2. ) # Write down the equations F = u*v*dx + dt*dot(c*grad(u), grad(v))*dx - (u_n + dt*f)*v*dx a, L = lhs(F), rhs(F) # Solve u = Function(V) t = 1e-12 file = File('one_domain_pulsed.pvd') f.t = t solve(a == L, u, bc) u_e = interpolate(u_D, V) error = np.abs(u_e.vector().array() - u.vector().array()).max() print('t = %.3g: error = %.3g' % (t, error)) file << (u, t) plot(u)
I have updated my Docker image using the command
curl -s https://get.fenicsproject.org | bashfound in the Quickstart section of the FEniCSContainers documentation (that I also used previously). In that fresh image, I got the same null result...
However, using the command
docker pull quay.io/fenicsproject/stable:latestthen
docker run -ti quay.io/fenicsproject/stable:latestI finally got it working!
What is the difference between the two commands / images? Maybe the documentation at https://fenics.readthedocs.io/projects/containers/en/latest/quickstart.html#quickstart should be updated to reflect that you do not get the same Docker images.