Howdy,
So I’m trying to write a simple example code of solving the 2d heat equation with source i.e
u_t = \alpha \Delta u +s. I’m doing this with the method of manufactured solution where I assume u_{T}(x,y,t) = \sin(xy(t+1)), and enforce neurmann boundary conditions on a box such that the neumann boundary condition is consistent with the true manufactured solution.
I’m doing a semi-discrete approach, where first I discretized in space with fenics, then solve an ODE in time. If I assume that the solution takes the form u(x,y,t)=\sum_{j} c_j(t) \phi_j(x,y), then the ODE system looks like M_1 \dot{c}(t) = M_2c(t) + d(t), where \bigg[M_1\bigg]_{i,j} = \int_{\Omega} \phi_i \phi_j d\Omega, \bigg[M_2\bigg]_{i,j} = -\alpha \int_{\Omega} \nabla\phi_i \cdot \nabla \phi_j d\Omega, and d_i(t) = \int_{\Omega} s \phi_i d \Omega + \alpha \int_{\partial \Omega_1} -(u_T)_x \phi_i d \partial \Omega_1 ++\alpha \int_{\partial \Omega_2} (u_T)_y \phi_i d \partial \Omega_2+ \alpha \int_{\partial \Omega_3} (u_T)_x \phi_i d \partial \Omega_3 + \alpha \int_{\partial \Omega_4} -(u_T)_y \phi_i d \partial \Omega_4; where \partial \Omega_1, \partial \Omega_2, \partial \Omega_3, \partial \Omega_4, are the left, top, right, and bottom of the box respectively.
In my code when I do the convergence analysis if I set \alpha=0, then I obtain the correct order of convergence I would expect, but when \alpha >0 it is not correct, making me think that I am not doing the boundary integral correctly in my fenics?
When doing the ODE time stepping, I invert the M_1 matrix to solve the ODE system: \dot{c}(t) = M_1^{-1}\bigg(M_2c(t) + d(t)\bigg)
I defined the portions of the box:
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
and then mark them
spatial_mesh= RectangleMesh(p1,p2,size,size)
boundaries = MeshFunction("size_t", spatial_mesh,spatial_mesh.topology().dim()-1)
boundaries.set_all(0)
left = Left()
right = Right()
bottom = Bottom()
top = Top()
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
I then in the source term calculate the boundary integral:
def d_1(t,source_exp,phi_i,alpha,true_x,true_y):
source_exp.t =t
true_x.t=t
true_y.t=t
ac = Constant(alpha)
boundary = ac*(-true_x*phi_i*ds(1) + true_y*phi_i*ds(2) +true_x*phi_i*ds(3) -true_y*phi_i*ds(4))
boundary_vec = assemble(boundary)
boundary_vec = boundary_vec.get_local()
term = source_exp*phi_i*dx
Term = assemble(term)
Term = Term.get_local()
rhs = np.zeros(len(Term))
for i in range(0,len(rhs)):
rhs[i]= Term[i] + boundary_vec[i]
return rhs
With picking the component of the spatial gradient \nabla u to be consistent with what \frac{\partial u}{\partial n} would be for the outward normal vector of the boundary I am on.
Am I not thinking about this correctly somehow?
Here is my complete code:
from fenics import *
import numpy as np
import math
from numpy.linalg import inv
from scipy.integrate import odeint
from scipy.integrate import solve_ivp
from scipy import *
import scipy
import time
from prettytable import PrettyTable
import pylab as plt
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
class NumericalSolution(UserExpression):
def __init__(self,time_step,c_list, basis_funcs,V, **kwargs):
super().__init__(**kwargs)
self.c_list = c_list
self.V=V
self.time_step=time_step
basis_vals = np.zeros(V.dim())
self.basis_vals=basis_vals
self.basis_funcs=basis_funcs
def eval(self, values, x):
V=self.V
value=0.0
for i in range(0,V.dim()):
value = value+self.c_list[self.time_step,i]*self.basis_funcs[i](x)
values[0] = value
def value_shape(self):
return (1,)
def basis_func_derv_i(i,x,n,tree,V,spatial_mesh):
cell_index = tree.compute_first_entity_collision(Point(*x))
cell_global_dofs = V.dofmap().cell_dofs(cell_index)
for local_dof in range(0,len(cell_global_dofs)):
if(i==cell_global_dofs[local_dof]):
cell = Cell(spatial_mesh, cell_index)
values = np.array([0,])
return V.element().evaluate_basis_derivatives(local_dof,n,x,cell.get_vertex_coordinates(),cell.orientation())
# If none of this cell's shape functions map to the i-th basis function,
# then the i-th basis function is zero at x.
return 0.0
def basis_func(i,x,tree,V,spatial_mesh):
cell_index = tree.compute_first_entity_collision(Point(*x))
cell_global_dofs = V.dofmap().cell_dofs(cell_index)
for local_dof in range(0,len(cell_global_dofs)):
if(i==cell_global_dofs[local_dof]):
cell = Cell(spatial_mesh, cell_index)
values = np.array([0,])
return V.element().evaluate_basis(local_dof,x,
cell.get_vertex_coordinates(),
cell.orientation())
# If none of this cell's shape functions map to the i-th basis function,
# then the i-th basis function is zero at x.
return 0.0
def d_1(t,source_exp,phi_i,alpha,true_x,true_y):
source_exp.t =t
true_x.t=t
true_y.t=t
ac = Constant(alpha)
boundary = ac*(-true_x*phi_i*ds(1) + true_y*phi_i*ds(2) +true_x*phi_i*ds(3) -true_y*phi_i*ds(4))
boundary_vec = assemble(boundary)
boundary_vec = boundary_vec.get_local()
term = source_exp*phi_i*dx
Term = assemble(term)
Term = Term.get_local()
rhs = np.zeros(len(Term))
for i in range(0,len(rhs)):
rhs[i]= Term[i] + boundary_vec[i]
return rhs
def IC_b(ic_exp,phi_i,dx):
term = ic_exp*phi_i*dx
Term = assemble(term)
Term = Term.get_local()
return Term
def ode_prob(t,c,M1_inv,M2,d_fun):
return M1_inv@(M2@c +d_fun(t))
def heat_2d(a,b,c,d,t_max,n_x,alpha,source_exp,true_fun,order_pol,ic_exp,left_exp,right_exp,true_x,true_y):
p1 = Point(a,c)
p2 = Point(b,d)
np.set_printoptions(precision=17)
order_of_polynomial = order_pol
size=n_x
spatial_mesh= RectangleMesh(p1,p2,size,size)
boundaries = MeshFunction("size_t", spatial_mesh,spatial_mesh.topology().dim()-1)
boundaries.set_all(0)
left = Left()
right = Right()
bottom = Bottom()
top = Top()
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
tree = BoundingBoxTree()
tree.build(spatial_mesh)
V = FunctionSpace(spatial_mesh,"CG",order_of_polynomial)
basis_func_temp = lambda i,x: basis_func(i,x,tree,V,spatial_mesh)
list_of_basis_functions=[]
for i in range(0,V.dim()):
list_of_basis_functions +=[(lambda i: lambda x : basis_func_temp(i,x))(i),]
m1=phi_i*phi_j*dx
M = PETScMatrix()
assemble(m1, tensor=M)
Mmat = M.mat()
Mmat_dense = Mmat.convert("dense")
M1=Mmat_dense.getDenseArray()
# M1[0,:] = basis_list_a
# M1[-1,:]=basis_list_b
M1_inv=inv(M1)
ic_int =phi_i*phi_j*dx
ic_mat_pet= PETScMatrix()
assemble(ic_int, tensor=ic_mat_pet)
ICmat = ic_mat_pet.mat()
ICmat_dense = ICmat.convert("dense")
ic_mat = ICmat_dense.getDenseArray()
nalpha_f = Constant(-alpha)
m2 = nalpha_f*np.inner(grad(phi_i),grad(phi_j))*dx
# M2 = PETScMatrix()
# M2 = assemble(m2,tensor=M2)
# M2 = M2.mat()
# M2 = M2.convert("dense")
# M2=M2.getDenseArray()
B = PETScMatrix()
assemble(m2, tensor=B)
MBmat = B.mat()
MBmat_dense = MBmat.convert("dense")
M2=MBmat_dense.getDenseArray()
M2_dim = M2.shape
d = lambda t: d_1(t,source_exp,phi_i,alpha,true_x,true_y)
ode_func = lambda t,c: ode_prob(c,t,M1_inv,M2,d)
b_vec = IC_b(ic_exp,phi_i,dx)
IC= np.linalg.solve(ic_mat, b_vec)
delta_x = (b-a)/n_x
#n_t = math.ceil((t_max*2*alpha)/delta_x**2)
n_t=n_x
#n_t = math.ceil(4*alpha*n_x**2)
t_list = np.linspace(0,t_max,n_t)
sol = odeint(ode_func, IC,t_list)
max_abs_error = 0.0
for i in range(0,len(t_list)):
if(i==0 or i==len(t_list)-1):
cur_num_sol = NumericalSolution(i,sol,list_of_basis_functions,V,element=V.ufl_element())
cur_num_sol_u =interpolate(cur_num_sol, V)
true_fun = Expression("sin(x[0]*x[1]*(t+1))",t=t_list[i] ,degree =2)
error_cur = math.sqrt(assemble(np.inner(cur_num_sol_u - true_fun, cur_num_sol_u- true_fun)*dx))
print("error" + str((i,error_cur)))
if(error_cur > max_abs_error):
max_abs_error=error_cur
if(i==0):
plt.subplot(1, 2, 1)
plot(cur_num_sol_u)
plt.subplot(1, 2, 2)
u_exaxt = interpolate(true_fun,V)
plot(u_exaxt)
#plt.show()
plt.savefig(str(n_x)+"heat2dfirstprob.png")
file = File (str(n_x)+"numheat2dfirstprob.pvd")
file << cur_num_sol_u
file = File (str(n_x)+"exaheat2dfirstprob.pvd")
file << u_exaxt
plt.clf()
if(i==len(t_list)-1):
plt.subplot(1, 2, 1)
plot(cur_num_sol_u)
plt.subplot(1, 2, 2)
u_exaxt = interpolate(true_fun,V)
plot(u_exaxt)
plt.savefig(str(n_x)+"head2dlastprob.png")
file = File (str(n_x)+"numhead2dlastprob.pvd")
file << cur_num_sol_u
file = File (str(n_x)+"exthead2dlastprob.pvd")
file << u_exaxt
plt.clf()
return (n_x,max_abs_error)
def generate_convergence_history(error_list,n_list,file_name):
x = PrettyTable()
x.field_names = ["N", "Max Error","Convergence Order"]
x.add_row([n_list[0], error_list[0],"-"])
con_list = np.zeros(len(error_list)-1)
for i in range(0,len(error_list)-1):
con_order = math.log(error_list[i]/error_list[i+1])/math.log(2.0)
x.add_row([n_list[i+1], error_list[i+1],con_order])
con_list[i] = con_order
print(x)
print("Average Convergence Order:" + str(np.mean(con_list)))
f = open(file_name, "a")
f.write(x.get_string() + str("\n"))
f.write("Average Convergence Order:" + str(np.mean(con_list)))
f.close()
def heat_2d_main(a,b,c,d,t_max,alpha):
#true solution sin((t+1)x)
#source this then
true_fun = Expression("sin(x[0]*x[1]*(t+1))",t=0 ,degree =2)
source_exp = Expression("cos(x[0]*x[1]*(t+1))*x[0]*x[1] +a*sin(x[0]*x[1]*(t+1))*pow(t+1,2)*(pow(x[0],2)+pow(x[1],2))",t=0,a=alpha ,degree =2)
left_exp = lambda t: math.sin(a*(t+1))
right_exp = lambda t: math.sin(b*(t+1))
ic_exp = Expression("sin(x[0]*x[1])",t=0 ,degree =2)
true_x =Expression("cos(x[0]*x[1]*(t+1))*x[1]*(t+1)",t=0 ,degree =2)
true_y =Expression("cos(x[0]*x[1]*(t+1))*x[0]*(t+1)",t=0 ,degree =2)
#grid size 2^n for n=5....10 so 64 to 1024 grid points
order_pol=1
n_list=[]
error_list=[]
for i in range(2,5):
n_x = 2**i
print("solving 1d heat equation with n=" +str(n_x))
(n_x,max_error)=heat_2d(a,b,c,d,t_max,n_x,alpha,source_exp,true_fun,order_pol,ic_exp,left_exp,right_exp,true_x,true_y)
n_list.append(n_x)
error_list.append(max_error)
generate_convergence_history(error_list,n_list,"heat2d_out.txt")
heat_2d_main(0.0,1.0,0.0,1.0,1.0,2.0)