Heat equation with Delta Dirac function as boundary condition

[janduigamer13]

Hello there! Sorry for bothering you guys, but I would like some advices in what I’m doing now. I’m working with the heat equation, the following problem:

And, I’m just trying to apply a delta dirac function in this circle with other boundary conditions. I’ve done the following in fenics:

nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

class UnitaryCircle(SubDomain):
    def inside(self, x, on_boundary):
        r = 0.25  # Radius of the circle
        return on_boundary and near((x[0]-0.5)**2 + (x[1]-0.5)**2, r**2)

# First boundary equation
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)

unitary_circle = UnitaryCircle()

# Second boundary condition (Delta dirac expression)
b = assemble(Constant(1.0))
point_source = PointSource(V, unitary_circle, Constant(1.0))
point_source.apply(b)

# Mark boundary subdomians
circle = CompiledSubDomain("unitary_circle")
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

bcl1 = DirichletBC(V, u_D, left)
bcl2 = DirichletBC(V, u_D, right)
bcr = DirichletBC(V, b, circle)
bcs = [bcl1, bcl2, bcr]
...
 solve(a == L, u, bcs)

It is not compiling, I’ve got the following error:

Traceback (most recent call last):
  File "/home/jandui/Documents/MATH code/testes_3/heat_test.py", line 42, in <module>
    b = assemble(Constant(1.0))
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 64, in _create_dolfin_form
    raise TypeError("Invalid form type %s" % (type(form),))
TypeError: Invalid form type <class 'dolfin.function.constant.Constant'>

I would like to know, please if:

1. Is it correct the way I stated the subdomains and the functions? ;

2. Is it correct the way I stated the Delta Dirac Function? I Basically would like to have an “impulse of heat” in this subdomain, and I tired doing it using the Delta Dirac.

3. To state subdomains problem using different spaces, we need to interpolate all functions at the same space ?

[dokken]

@janduigamer13

b = assemble(Constant(1.0))

Unclear why this is needed and what you are trying to achieve.
If you want a rhs with only values 1, i would suggest that you create a dolfin.Function, and access the underlying vector and set all values to 1.

[janduigamer13]

Sorry for being unclear. My Right hand side in this problem, is not associated with Delta Dirac Function. I’ll try to explain in a better way. I have the heat equation and what I’m trying to do is use different Dirichlet boundary conditions in this equation. Specifically, I’m trying to set the corners of the square with a heat source stated by me as this u_d expression in this domain, and in a circle of radius 0.25 centered at this square domain, an impulse represented by the Delta Dirac function. Basically, what I tried with this b, is based on the topic:

Point Sources (Redux)

Following my post above, with a minor adaptations, a point source can be implemented as follows: # Create a point source for Poisson problem # Author: Jørgen S. Dokken # SPDX-License-Identifier: MIT from mpi4py import MPI from petsc4py import PETSc import dolfinx import dolfinx.fem.petsc import numpy as np import ufl def compute_cell_contributions(V, points): # Determine what process owns a point and what cells it lies within mesh = V.mesh _, _, owning_points, cells = dolfinx.cp…

But I’m trying to change this of acting in a point to act on this circular subdomain. And, sorry, I have noticed now that my formulation is not correct, it is supposed to be a function f at the right hand side of Heat Equation:

The function f is actually:

f = Constant(beta - 2 - 2*alpha)

[dokken]

Please modify your code appropriately, and make sure that the code can reproduce the error. Just be reading error messages and snippets of your code it is really hard to give you any guidance.

[janduigamer13]

Ok, To reproduce the error, you have to run the following:

from __future__ import print_function
from dolfin import File
from fenics import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

T = 2.0            # final time
num_steps = 20     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

class UnitaryCircle(SubDomain):
    def inside(self, x, on_boundary):
        r = 0.25  # Radius of the circle
        return on_boundary and near((x[0]-0.5)**2 + (x[1]-0.5)**2, r**2)

# First boundary equation
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)

unitary_circle = UnitaryCircle()

# Second boundary condition (Delta dirac expression)
b = assemble(Constant(1.0))
point_source = PointSource(V, unitary_circle, Constant(1.0))
point_source.apply(b)

# Mark boundary subdomians
circle = CompiledSubDomain("unitary_circle")
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

bcl1 = DirichletBC(V, u_D, left)
bcl2 = DirichletBC(V, u_D, right)
bcr = DirichletBC(V, b, circle)
bcs = [bcl1, bcl2, bcr]

# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0

solutions = []
for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t

    # Compute solution
    solve(a == L, u, bcs)

    solutions.append(u.copy())
    # Plot solution
    plot(u)

    # Compute error at vertices
    u_e = interpolate(u_D, V)
    error = np.abs(u_e.vector().get_local() - u.vector().get_local()).max()
    print('t = %.2f: error = %.3g' % (t, error))

    # Update previous solution
    u_n.assign(u)

[dokken]

You are still keeping the

b = assemble(Constant(1.0))
point_source = PointSource(V, unitary_circle, Constant(1.0))
point_source.apply(b)

which I already told you does not make sense.

As this is a continuous point source, you would have to make a discrete choice of the number of point sources you want in your domain.

What I would do is to make your mesh in a way that conforms with the boundary of the circle with radius sqrt(0.25) from 0.5, 0.5.
Then you can set the bcs explicitly (as any other bc) on the degrees of freedom on this boundary.

[janduigamer13]

Ok, thanks. I tried the following, as my problem was defining the Delta Dirac function, I found one discussion from the Forum about the implementation, and I used it:

from __future__ import print_function

from dolfin import File
from dolfin import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

T = 2.0            # final time
num_steps = 20     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 20
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# First boundary equation
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0)

# Mark boundary subdomians

class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs)

    def eval(self, values, x):
        eps = self.eps
        values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)

    def value_shape(self): return ()

delta = Delta(eps=1E-4, x0=np.array([0.5, 0.5]), degree=5)
f = delta

def boundary_point1(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], x[0]+tol,tol)

def boundary_point1(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0]+0.2, x[0]+tol+0.2,tol)

def boundary_L(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 0 ,tol)

def boundary_R(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 1 ,tol)

bcl1 = DirichletBC(V, u_D, boundary_L)
bcl2 = DirichletBC(V, u_D, boundary_R)
bcr1 = DirichletBC(V, f, boundary_point1)
bcr2 = DirichletBC(V, f, boundary_point1)

bcs = [bcl1, bcl2, bcr1,bcr2]
# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0

solutions = []
for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t

    # Compute solution
    solve(a == L, u, bcs)

    solutions.append(u.copy())
    # Plot solution
    plot(u)

    # Compute error at vertices
    u_e = interpolate(u_D, V)
    error = np.abs(u_e.vector().get_local() - u.vector().get_local()).max()
    print('t = %.2f: error = %.3g' % (t, error))

    # Update previous solution
    u_n.assign(u)

But the error is really big, I got it:

Solving linear variational problem.
t = 2.00: error = 7.4

My question is: Since the error is big (it started with 5.12 and was increasing at each time step), is the way I impose the boundary conditions not correct? I tried to choose a single point using this “boundary_points”. Also, to see the solution in Paraview, i need to write:

vtkfile = File("heat_equation_solution.pvd")
vtkfile << u

Outside the time step loop? Because I’m doing this and it is not saving.
Also, really thanks for your time and attention and sorry for sending a lot of messages. Programmation is not my strong side and I’m trying to do my best to learn it.

[dokken]

Are you now trying to have a single dirac delta function? If you have a dirac delta function at a single point, you could have a look at:

Point Sources (Redux)

Following my post above, with a minor adaptations, a point source can be implemented as follows: # Create a point source for Poisson problem # Author: Jørgen S. Dokken # SPDX-License-Identifier: MIT from mpi4py import MPI from petsc4py import PETSc import dolfinx import dolfinx.fem.petsc import numpy as np import ufl def compute_cell_contributions(V, points): # Determine what process owns a point and what cells it lies within mesh = V.mesh _, _, owning_points, cells = dolfinx.cp…

that uses a single point source in either dolfin or dolfinx.