point source in RHS of heat eq

Question

I need to solve heat equation but source term includes delta function and I couldn't add it fenics although I saw code PointSource. I also try it with Gauss (delta(x)=(exp(-x^2/epsilon^2)/(epsilon*sqrt(pi)). the simplest example

u_t = u_rr + u_psi psi + delta(r)delta(psi)sigma
r [0,R] and psi is [-pi,pi)
u(R,psi,t)=0 and u(r,psi,0)=0

thanks

Answer 1

Hi, what do you mean by I couldn't add it to fenics although...? The following uses PointSource to represent a point charge and computes its electric potential. Maybe you find it helpful

from dolfin import *

mesh = RectangleMesh(-1, -1, 1, 1, 100, 100)
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

a = inner(grad(u), grad(v))*dx
L = Constant(0)*v*dx
bc = DirichletBC(V, Constant(0), DomainBoundary())
A, b = assemble_system(a, L, bc)

delta = PointSource(V, Point(0., 0.,), 100)
delta.apply(b)

u = Function(V)
solve(A, u.vector(), b)

plot(u, interactive=True) 

Comments

Thank you MiroK for really helpful reply,

I am new in Fenics and Python. your solving is for poission right such that

-Nabla^2 u = delta(x)*delta(y)

but how can I improve f (f is the function in the RHS of equation)

for example multiply f by x
f(x,y)=delta(x)delta(y)x

or for heat equation

f(x,y,t)=delta(x)delta(y)g(t)
(t is time)

Thank you again

---

You can do that by modifying the third argument of PointSource which determines the
magnitude of the delta-function. Consider this silly example

from dolfin import *
from math import sin, cos

mesh = RectangleMesh(-1, -1, 1, 1, 100, 100)
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

a = inner(grad(u), grad(v))*dx
L = Constant(0)*v*dx
bc = DirichletBC(V, Constant(0), DomainBoundary())
A, b = assemble_system(a, L, bc)

u = Function(V)
n_steps = 20
dt = 1./n_steps
t = 0
for step in range(n_steps):
    cn = cos(2*pi*step*dt)
    sn = sin(2*pi*step*dt)
    den = 1 + sn**2

    x = 0.5*cn/den
    y = 0.5*cn*sn/den
    t += dt

    delta = PointSource(V, Point(x, y), abs(x*y)*t*100)
    delta.apply(b)
    solve(A, u.vector(), b)

    plot(u, interactive=True)

    # reset rhs
    b.zero()
    bc.apply(b)