Discarding Point Source by applying Negative Point Source

1

Dear all,

I have a time-dependent thermal problem. The boundary conditions are not time-dependent but there is a moving heat source, that I apply each iteration. I’m trying to find a way to make the code more efficient by trying to avoid reassembly of the system, namely b. I have tried two different approaches which yield the same wrong result.

Firstly I tried calculating the mass matrix and the stiffness matrix and recalculating “b” via simple matrix operations rather than assembly, using the reference here; Time-Dependent Problems — A FEniCS Tutorial 1.0 documentation. Unfortunately, this didn’t work; the thermal effect of previous point sources was missing.

Secondly, I applied a more straightforward method, applying the negative of the point source after the system was solved. Here is a simplified working example;

from fenics import *
import os
import numpy as np


num_steps = 100
#define simulation parameters
rho = 4420 #kg/m^3 
c_p = 750 #J/kgK
k = 27 #W/mK
alpha = k/(rho*c_p) #m^2/s
power = 175 #W
tol = 1e-14
dt = 1.6/100


mesh = BoxMesh(Point(-1,-1,-1), Point(1,1,1), 20,20,20)

#determine the limits of the domain in m
x_max = np.max(mesh.coordinates()[:,0])
x_min = np.min(mesh.coordinates()[:,0])
y_max = np.max(mesh.coordinates()[:,1])
y_min = np.min(mesh.coordinates()[:,1])
z_max = np.max(mesh.coordinates()[:,2])
z_min = np.min(mesh.coordinates()[:,2])


#define markers for the boundary faces
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 4)

#define boundary subclasses for bottom, top and side surfaces, mark the boundaries
class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], z_min, tol)
b_bottom = BoundaryBottom()
b_bottom.mark(boundary_markers, 0)

class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], z_max, tol)
b_top = BoundaryTop()
b_top.mark(boundary_markers, 1)

class BoundarySides(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0], x_min, tol) or near(x[0], x_max, tol) or near(x[1], y_min, tol) or near(x[1], y_max, tol))
b_sides = BoundarySides()
b_sides.mark(boundary_markers, 2)

#define function space for v
V = FunctionSpace(mesh, "CG", 1)

#redefine d in therm of the boundary markers
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

#define Dirichlet initial condition
T_d = Constant(0)
#define Neumann boundary conditions
g_top = Constant(0.0)
g_sides = Constant(0)

#list of boundary conditions for convenience
boundary_conditions = {0: {'Dirichlet': T_d},
                       1: {'Neumann': g_top},
                       2: {'Neumann': g_sides}}

#interpolate the initial condition
T_0 = interpolate(Constant(0.0), V)

#define test and trial function and the source term
T = TrialFunction(V)
v = TestFunction(V)
#f = Constant(0)

#gather all Dirichlet boundary conditions in one list
bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
        bcs.append(bc)
         
#gather all Neumann boundary conditions in weak form
integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        g = boundary_conditions[i]['Neumann']
        integrals_N.append(g*v*ds(i))
            

#write the weak form and seperate into right- and left-hand sides
F = T*v*dx + alpha*dt*dot(grad(T), grad(v))*dx - T_0*v*dx - alpha*dt*sum(integrals_N)
a, L = lhs(F), rhs(F)

#define the function and time
T = Function(V)
t = 0

xdmf_file = XDMFFile(MPI.comm_world, "results/temperature.xdmf")
xdmf_file.parameters["flush_output"] = True
xdmf_file.parameters["functions_share_mesh"] = True
xdmf_file.parameters["rewrite_function_mesh"] = False

A, b = assemble_system(a, L, bcs)

x = np.linspace(-0.8, 0.8, 100)
y = np.linspace(-0.8, 0.8, 100)


#solve the weak form for each timestep
for i in range(num_steps):
    
    #assemble the equations
    #b = assemble(L, tensor=b)

    try:
        delta = PointSource(V, Point(x[i], y[i], z_max-tol), dt*power/(rho*c_p)) 
        delta.apply(b)
    except RuntimeError:
        pass

    #solve the system   
    solve(A, T.vector(), b, 'cg', 'sor')


    #assign the temperature as the initial condition for the next step
    T_0.assign(T)
    
    #store the temperature at points of interests
    try:
        delta = PointSource(V, Point(x[i], y[i], z_max-tol), -dt*power/(rho*c_p)) 
        delta.apply(b)
    except RuntimeError:
        pass
    #write temperature data
    T.rename("Temperature", "FEM Temperature")
    xdmf_file.write(T, t)
    t += dt

This approach only works if the line

b = assemble(L, tensor=b)

is uncommented. But then, I reassemble L and recalculate b at every iteration, which defies my initial purpose of efficiency. I think when the line. I checked whether b is the same when a point source is added and “subtracted” and it is the same.

There is definitely something funny going on that I cannot understand and I don’t know much about what’s happening under the hood in dolfin. If anyone encountered a similar problem or has experience with PointSource applications, I would appreciate every bit of help.

Cheers,

rek

2

The reason this MWE doesn’t work without re-assembling b is that the update to the previous time step solution T_0 (i.e., T_0.assign(T)) is not reflected in the vector b, which was originally assembled with T_0 equal to the initial condition. I believe a procedure like the linked tutorial should be applicable to this problem. The “memory” of previous point sources should enter into the current time step through T_0, so you might look for a bug related to the update of T_0.