Discretize with FEM in space and then do ODE in time : mass matrix changing during time stepping?

Ryan_Vogt2 · July 5, 2022, 5:08pm

Howdy,

So as a warm up I tried to solve the 1-d heat equation where I say that u(x,t) = \sum c_j(t) \phi_j(x). I do the FEM in space, and generate the matices associated for each term in the final system of ODE: M_1 \dot{c} = M_2c + d(t).

When I do the time stepping just doing a simple euler:

	for i in range(1,len(t_list)):
		t = t_list[i-1]

		print(np.linalg.cond(M1))
		next_sol = np.linalg.solve(M1,dt*(M2@sol[i-1,:] +d(t)) + (M1@sol[i-1,:]))
		sol[i,:]=next_sol
		# print(M2@sol[i-1,:])

3.9938387318170525

3.9938387318170525

3.9938387318170525

5960.4321111585605

70872.20443761055

70872.20443761055

70872.20443761055

70872.20443761055

70872.20443761055

70872.20443761055

70872.20443761055

70872.20443761055

70872.20443761055

70872.20443761055

141053.0689730364

141053.0689730364

4108.854366444064

3570269581.8684115

3570269581.8684115

3570269581.8684115

3570269581.8684115

3570269581.8684115

20314608610026.54

558294360160.5804

558294360160.5804

71401822674.21254

71401822674.21254

23072194296.435165

23072194296.435165

5648302300952801.0

5648302300952801.0

5648302300952801.0

5648302300952801.0

5648302300952801.0

5648302300952801.0

3.241721984318394e+18

3.241721984318394e+18

2071418921090024.0

2071418921090024.0

1.8281517530721127e+22

1.8281517530721127e+22

1.8281517530721127e+22

1.8281517530721127e+22

1.8281517530721127e+22

7.455156198284555e+28

7.455156198284555e+28

7.455156198284555e+28

1.8370838355398053e+28

1.8370838355398053e+28

1.8370838355398053e+28

1.8370838355398053e+28

1.8370838355398053e+28

1.8370838355398053e+28

1.8643064584435324e+27

1.8643064584435324e+27

1.8643064584435324e+27

1.8643064584435324e+27

1.8643064584435324e+27

2.9459028094408635e+29

6.347724126735073e+30

5.915766739658558e+36

5.915766739658558e+36

5.915766739658558e+36

5.658535187497739e+40

5.658535187497739e+40

5.658535187497739e+40

5.658535187497739e+40

5.658535187497739e+40

5.658535187497739e+40

5.658535187497739e+40

5.658535187497739e+40

4.949629828329596e+44

4.949629828329596e+44

4.949629828329596e+44

4.949629828329596e+44

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

1.4172700791877978e+24

2.5282685613029403e+28

2.5282685613029403e+28

2.5282685613029403e+28

8.315960439107994e+34

8.315960439107994e+34

8.315960439107994e+34

8.315960439107994e+34

8.315960439107994e+34

8.315960439107994e+34

8.315960439107994e+34

8.315960439107994e+34

2.3362140819840576e+43

2.3362140819840576e+43

9.083946732855834e+60

9.083946732855834e+60

1.4437317885697556e+65

5.635664062124686e+59

5.635664062124686e+59

5.635664062124686e+59

5.635664062124686e+59

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

4.548278265376804e+69

2.33032614244715e+72

3.356795993965915e+83

1.7513137109850495e+97

4.2422509650004495e+104

4.2422509650004495e+104

4.2422509650004495e+104

4.2422509650004495e+104

4.2422509650004495e+104

4.2422509650004495e+104

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

5.065573093182438e+125

inf

5.368240868941827e+177

5.368240868941827e+177

5.368240868941827e+177

1.5074951988067158e+196

1.816253544677133e+214

inf

inf

inf

inf

inf

inf

inf

inf

I am extremely confused because each iteration the condition number of the matrix M1 is changing, even though I am not editing it anywhere during the time stepping. (which means also M1) is changing itself each iteration somehow. This is extremely strange behavior to me?

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 import *
import scipy

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
		# 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):
	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])
	return Term

def IC_b(ic_exp,phi_i):

	term = ic_exp*phi_i*dx
	Term = assemble(term)
	Term = Term.get_local()



	return Term


def ode_prob(c,t,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))



def heat_1d(a,b,t_max,n_x,alpha,source_exp,true_fun,order_pol,ic_exp,left_exp,right_exp):
	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)
	

	phi_i=TestFunction(V)
	phi_j = TrialFunction(V)
	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])

	m1=phi_i*phi_j*dx

	M = PETScMatrix()
	assemble(m1, tensor=M)

	Mmat = M.mat()
	

	Mmat_dense = Mmat.convert("dense")
	print("Dense")
	M1=Mmat_dense.getDenseArray()
	print(M1)
	M1_inv=inv(M1)

	alpha_f = p = Constant(alpha)
	m2 = alpha_f*phi_ixx*phi_j*dx
	# M2 = PETScMatrix()
	# M2 = assemble(m2,tensor=M2)
	# M2 = M2.mat()
	# M2 = M2.convert("dense")
	# M2=M2.getDenseArray()


	M = PETScMatrix()
	assemble(m2, tensor=M)
	Mmat = M.mat()
	Mmat_dense = Mmat.convert("dense")
	M2=Mmat_dense.getDenseArray()
	
	



	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)

	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])
	

	d = lambda t: d_1(t,source_exp,phi_i,first_derv_list_a,first_derv_list_b,left_exp,right_exp,alpha)


	ode_func = lambda c,t: ode_prob(c,t,M1_inv,M2,d)
	b_vec = IC_b(ic_exp,phi_i)
	#print(b_vec)

	#IC = M1_inv@b_vec
	IC= np.linalg.solve(M1, 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*3*alpha)/delta_x**2)
	n_t=2*n_x
	t_list = np.linspace(0,t_max,num=10000)

		

	#sol = odeint(ode_func, IC,t_list)
	sol = np.zeros((len(t_list),len(IC)))
	sol[0,:]=IC
	dt = t_list[1]-t_list[0]
	print("mat")
	print(M1)
	#print(np.linalg.det(M1))
	
	#quit()
	comp = M1
	for i in range(1,len(t_list)):
		t = t_list[i-1]

		print(np.linalg.cond(M1))
		next_sol = np.linalg.solve(M1,dt*(M2@sol[i-1,:] +d(t)) + (M1@sol[i-1,:]))
		sol[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)
	# 	error_cur = math.sqrt(assemble(inner(cur_num_sol_u - true_fun, cur_num_sol_u- true_fun)*dx))
		
	# 	if(error_cur > max_abs_error):
	# 		max_abs_error=error_cur


	# print(max_abs_error)

	#t_list = sol.t
	
	# print(IC)
	# print(sol.shape)
	# print(M1.shape)
	# print(n_x)
	# print(sol[:,0])
	# quit()

	return 



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
	for i in range(5,11):

		n_x = 2**i
		print("solving 1d heal equation with n=" +str(n_x))
		heat_1d(a,b,t_max,n_x,alpha,source_exp,true_fun,order_pol,ic_exp,left_exp,right_exp)



heat_1d_main(0.0,1.0,1.0,1.0)

Ryan_Vogt2 · July 6, 2022, 4:31pm

Found the problem - it seems to be a reference bug in PETSC

What I had before:

	M = PETScMatrix()
	assemble(m1, tensor=M)

	Mmat = M.mat()
	

	Mmat_dense = Mmat.convert("dense")
	print("Dense")
	M1=Mmat_dense.getDenseArray()
	print(M1)
	M1_inv=inv(M1)

	alpha_f = p = Constant(alpha)
	m2 = alpha_f*phi_ixx*phi_j*dx
	# M2 = PETScMatrix()
	# M2 = assemble(m2,tensor=M2)
	# M2 = M2.mat()
	# M2 = M2.convert("dense")
	# M2=M2.getDenseArray()


	M = PETScMatrix()
	assemble(m2, tensor=M)
	Mmat = M.mat()
	Mmat_dense = Mmat.convert("dense")
	M2=Mmat_dense.getDenseArray()

I made a PETSC matrix and called it M, to create M1 and then made another, which again I called M, and then made M2. Looks like some weird referencing thing is occuring where the pointer was pointing to wrong locations.

When I make unique variable names i.e now:

	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_inv=inv(M1)

	alpha_f = p = Constant(alpha)
	m2 = alpha_f*phi_ixx*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()

Things are fine, and the matrix doesn’t change during time stepping.

Thanks.

Topic		Replies	Views
Convergence order of semi-distretization not yielding reasonable convergence order, yet the numerical and true solution look good	4	525	July 19, 2022
FEM for Space-Time variational formulation	1	571	March 23, 2022
Code for DG time stepping scheme for heat equation General	1	207	January 24, 2024
Heat equation - error for longer times	0	408	March 19, 2021
Heat equation with explicit euler General	10	425	July 1, 2024

Discretize with FEM in space and then do ODE in time : mass matrix changing during time stepping?

Related topics