Howdy,
So as a warm up I solved the 1d heat equation. First I represent the solution as u= \sum c_j(t)\phi_j(x). Next I do a FEM discretization in space (using first-order CG), resulting in an ODE system for the time dependent coifficents. I solve this problem, for increassing finer meshes (each one have 2x the number of vertices as before). I then calculate the convergence or i.e \frac{ln(\frac{E_{n_x}}{E_{2n_x}})}{ln(2)}. I do this up to 256 grid points in space ( I choose \Delta t=\Delta x for the time integration), and produce the following convergence table:
+-----+---------------------+-------------------------+
| N | Max Error | Convergence Order |
+-----+---------------------+-------------------------+
| 32 | 0.26974518677568843 | - |
| 64 | 0.26985983054165796 | -0.0006130261573324018 |
| 128 | 0.2698884906012833 | -0.00015321113308644212 |
| 256 | 0.26989568274959913 | -3.844527950794879e-05 |
+-----+---------------------+-------------------------+
Average Convergence Order:-0.0002682275233089309
It looks like the error is not improving (I would expect for a smooth problem that the convergence order would begin to go towards 2, since I used linear basis functions, but it isn’t changing. When I look at a plot of the numerical solution versus the true solution, then look good in the eyeball norm: here is a comparison of the true and exact solution for 256 points
This makes me believe that somehow I’m calulating the error incorrectly. Here I am constructing my own expression at each time step that represents the solution at time step n i.e u(x,t_n) = \sum c_j(t_n) \phi_j(x). I am then comparing it with the true solution at t_n in the L2 norm. I then take the maximum of all these errors to be used in the congerence order calculation for the error at a spatial discretization of n_x. Is this how one would do it correctly?
max_abs_error = 0.0
for i in range(0,len(t_list)):
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]*(t+1))",t=t_max ,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==len(t_list)-1):
if(i==0):
plot(cur_num_sol_u)
#plt.show()
plt.savefig(str(n_x)+"firstprob.pdf")
plt.clf()
if(i==len(t_list)-1):
plot(cur_num_sol_u, )
u_exaxt = interpolate(true_fun,V)
#plt.show()
plot(u_exaxt,linestyle="dotted")
plt.legend(["numerical","exact"])
plt.savefig(str(n_x)+"lastprob.png")
plt.clf()
Here is my full code if you want to look at it - thank you:
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 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)
#print("value:" + str(value))
#print("c" +str((self.time_step,self.c_list.max())))
#print(self.c_list)
values[0] = value
# print(value)
# quit()
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,first_derv_list_a,first_derv_list_b,left_exp,right_exp,alpha,left_der,right_der):
source_exp.t =t
term = source_exp*phi_i*dx
Term = assemble(term)
Term = Term.get_local()
# Term = PETScMatrix()
# Term = assemble(term,tensor=Term)
# Term =Term.mat()
# Term = Term.convert("dense")
# Term=Term.getDenseArray()
# for i in range(0,len(Term)):
# Term[i] = Term[i] + alpha*(left_exp(t)*first_derv_list_a[i] - right_exp(t)*first_derv_list_b[i])
Term[0] = left_der(t)
Term[-1] = right_der(t)
return Term
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):
# print(M2@c)
# print(d_fun(t))
# print(M1_inv@(M2@c +d_fun(t)))
# print(type(M1_inv@(M2@c +d_fun(t))))
return M1_inv@(M2@c +d_fun(t))
#return d_fun(t)
def heat_1d(a,b,t_max,n_x,alpha,source_exp,true_fun,order_pol,ic_exp,left_exp,right_exp):
np.set_printoptions(precision=17)
order_of_polynomial = order_pol
spatial_mesh= IntervalMesh(n_x,a,b)
tree = BoundingBoxTree()
tree.build(spatial_mesh)
V = FunctionSpace(spatial_mesh,"CG",order_of_polynomial)
left_exp_der = lambda t: a*math.cos(a*(t+1))
right_exp_der = lambda t: b*math.cos(b*(t+1))
phi_i=TestFunction(V)
phi_j = TrialFunction(V)
phi_jx = phi_j.dx(0)
phi_ix = phi_i.dx(0)
phi_ixx=phi_ix.dx(0)
# t=1.0
# print(len(d_1(t,source_exp,phi_i)))
# print(d_1(t,source_exp,phi_i)[-1])
first_derv = lambda i,x:basis_func_derv_i(i,x,1,tree,V,spatial_mesh)
basis_func_temp = lambda i,x: basis_func(i,x,tree,V,spatial_mesh)
first_derv_list_a=np.zeros(V.dim())
first_derv_list_b=np.zeros(V.dim())
basis_list_a=np.zeros(V.dim())
basis_list_b=np.zeros(V.dim())
list_of_basis_functions_der=[]
list_of_basis_functions=[]
point_a = np.array([a])
point_b = np.array([b])
for i in range(0,V.dim()):
list_of_basis_functions_der += [(lambda i: lambda x : first_derv(i,x))(i),]
list_of_basis_functions +=[(lambda i: lambda x : basis_func_temp(i,x))(i),]
first_derv_list_a[i] = list_of_basis_functions_der[i](point_a)
first_derv_list_b[i] = list_of_basis_functions_der[i](point_b)
basis_list_a[i] = list_of_basis_functions[i](point_a)
basis_list_b[i] = list_of_basis_functions[i](point_b)
just_phi = phi_i*dx
phi_assemb = assemble(just_phi)
phi_vec = phi_assemb.get_local()
first_row = np.zeros(len(phi_vec))
last_row = np.zeros(len(phi_vec))
for i in range(0,len(first_row)):
first_row[i] = basis_list_a[i]*phi_vec[i]
last_row[i] = basis_list_b[i]*phi_vec[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*phi_ix*phi_jx*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
for i in range(0,M2_dim[0]):
for j in range(0, M2_dim[1]):
M2[i,j] = M2[i,j] + alpha*(basis_list_b[i]*first_derv_list_b[j] - basis_list_a[i]*first_derv_list_a[j])
zero_row = np.zeros(M2_dim[1])
M2[0,:] = zero_row
M2[-1,:]=zero_row
# (w,v)= np.linalg.eig(M1_inv*M2)
# print(M2)
# print(w)
# quit()
d = lambda t: d_1(t,source_exp,phi_i,first_derv_list_a,first_derv_list_b,left_exp,right_exp,alpha,left_exp_der,right_exp_der)
#ode_func = lambda c,t: ode_prob(c,t,M1_inv,M2,d)
ode_func = lambda t,c: ode_prob(c,t,M1_inv,M2,d)
b_vec = IC_b(ic_exp,phi_i,dx)
#print(b_vec)
#IC = M1_inv@b_vec
IC= np.linalg.solve(ic_mat, b_vec)
#(IC,info)= scipy.sparse.linalg.gmres(M1,b_vec,maxiter=10000)
#print(info)
delta_x = (b-a)/n_x
n_t = math.ceil((t_max*2*alpha)/delta_x**2)
n_t=n_x
t_list = np.linspace(0,t_max,n_t)
sol = odeint(ode_func, IC,t_list)
#sol_ode = solve_ivp(ode_func, (0, t_max), IC, method='Radau')
# for i in range(0,len(t_list)):
# print("--------solution t=" +str(t_list[i])+"-----------")
# print(sol[i,:])
# print("-------------------------------------------------")
# t_list = sol_ode.t
# sol= np.zeros((len(t_list),len(IC)))
# print(sol)
# for i in range(0,len(IC)):
# sol[:,i] = sol_ode.y[i]
# sol = np.zeros((len(t_list),len(IC)))
# sol[0,:]=IC
# dt = t_list[1]-t_list[0]
#print(np.linalg.det(M1))
#quit()
# comp = M1
# for val_i in range(1,len(t_list)):
# t_cur = t_list[val_i-1]
# #print(np.linalg.cond(M1))
# #print(M1.max())
# # print(np.all(M1==comp))
# # comp = M1
# next_sol = np.linalg.solve(M1,dt*(M2@sol[val_i-1,:] +d(t_cur)) + M1@sol[val_i-1,:])
# # next_sol= scipy.sparse.linalg.gmres(M1,dt*(M2@sol[val_i-1,:] +d(t_cur)) + (M1@sol[val_i-1,:]))
# #next_sol= M1_inv*(dt*(M2@sol[val_i-1,:] +d(t_cur)) + (M1@sol[val_i-1,:]))
# # M1=comp1
# # M2=comp2
# print("--------solution ----------")
# print("--------t=" +str(t_list[val_i]) +"----------")
# print("--------order_pol=" +str(order_of_polynomial) +"----------")
# print(next_sol)
# print("------------------")
# time.sleep(1)
# sol[val_i,:]=next_sol
# # print(M2@sol[i-1,:])
max_abs_error = 0.0
for i in range(0,len(t_list)):
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]*(t+1))",t=t_max ,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==len(t_list)-1):
if(i==0):
plot(cur_num_sol_u)
#plt.show()
plt.savefig(str(n_x)+"firstprob.pdf")
plt.clf()
if(i==len(t_list)-1):
plot(cur_num_sol_u, )
u_exaxt = interpolate(true_fun,V)
#plt.show()
plot(u_exaxt,linestyle="dotted")
plt.legend(["numerical","exact"])
plt.savefig(str(n_x)+"lastprob.png")
plt.clf()
#t_list = sol.t
# print(IC)
# print(sol.shape)
# print(M1.shape)
# print(n_x)
# print(sol[:,0])
# quit()
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_1d_main(a,b,t_max,alpha):
#true solution sin((t+1)x)
#source this then
true_fun = Expression("sin(x[0]*(t+1))",t=0 ,degree =2)
source_exp = Expression("cos(x[0]*(t+1))*x[0] +a*sin(x[0]*(t+1))*pow(t+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])",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(5,9):
n_x = 2**i
print("solving 1d heat equation with n=" +str(n_x))
(n_x,max_error)=heat_1d(a,b,t_max,n_x,alpha,source_exp,true_fun,order_pol,ic_exp,left_exp,right_exp)
n_list.append(n_x)
error_list.append(max_error)
generate_convergence_history(error_list,n_list,"heat_out.txt")
heat_1d_main(0.0,1.0,1.0,2.0)
