Source term ineffective


306
views
0
8 months ago by

Hello!

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[0]-bw)/w, 2)) * exp(-pow(((t-t0)/tau), 2)) * exp(-alpha*x[1])',
               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)
Community: FEniCS Project
Check the visualization scale (or print u.vector().array()), if I run the code everything is as expected.
written 8 months ago by Michal Habera  
Thanks for your comment. However, no matter how closely I rescale, I always get the same thing. Also, running:
for i in u.vector().array():
    print('%.3e' % i)​


gives zero everywhere:

0.000e+00
0.000e+00
0.000e+00
0.000e+00
0.000e+00
0.000e+00
(...)
0.000e+00
0.000e+00
0.000e+00
0.000e+00
0.000e+00
written 8 months ago by François Lapointe  
I also think it's a scaling issue. Have you looked at your pvd file? If I run the code it looks as expected.

Check out the output of
np.any(u.vector().array())​

For me it returns true.

written 8 months ago by klunkean  
A scaling issue would make sense, but I get False as the output of that command.
written 8 months ago by François Lapointe  
This is a strange issue. You obviously get all-zeros function when you just initialized it, i.e. u = Function(V).
Really make sure, that you are executing the exact copy of what you posted here. Then try to ged rid of state-dependent issues (clear cache with command instant-clean). If any of that works try running the exact script in docker (stable) environment - http://fenics.readthedocs.io/projects/containers/en/latest/.
written 8 months ago by Michal Habera  
I copy-pasted the code from this post in a clean Jupyter Python2 notebook, and I still get the same result.

I tried the instant-clean command:
fenics@31d870d536fe:~$ instant-clean
Removing 48 modules from Instant cache...
Removing 3 error logs from Instant cache...
Removing 0 lock files from Instant cache...
with the same output...

I will try later on the latest Docker (stable) image and keep you posted.
written 8 months ago by François Lapointe  

1 Answer


2
8 months ago by
OK, I got something.

I have updated my Docker image using the command curl -s https://get.fenicsproject.org | bash found 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:latest then docker run -ti quay.io/fenicsproject/stable:latest I 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.

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

Similar posts:
Search »