Hi, Dear FEniCS community,
I have recently encountered a numerical error issue while solving a convection-diffusion equation, and I am seeking some help from the community.
Here’s the situation: In the equation,
there are velocity and diffusion coefficient parameters. However, when I set both the velocity and diffusion coefficient to 0, I noticed that the solution value during the iteration process showed a slight deviation from the initial value.
I searched and found that FEniCS uses float64 precision, which provides 15 significant digits. However, the deviation I observed during the iteration process is around the order of e-14, which does not seem to be caused by float precision limitations.
This phenomenon confuses me, and I’m not sure what could be causing this numerical error. I would greatly appreciate it if anyone with knowledge in this area could help explain this. Any suggestions are welcome!
Below is a snippet of the code for reference:
from fenics import *
from Finite_element_method.fem_helper import assign_initial_conditions
from Finite_element_method.fem_helper import return_solutions
import random
import numpy as np
import copy
random.seed(2)
"""
When both the speed and diffusion coefficient are set to 0,
the comparison of the solution value with the initial value during the iteration process.
"""
def fem_return_solutions():
solution_matrix = []
# parameters read
r = 0
f = 0
dt = 0.0001
row = 6
col = 6
if_keep_decimal_point = False
#set the initial grid
initial_grids = [[1 for j in range(row)] for i in range(col)]
for i in range(2, 4):
for j in range(2, 4):
initial_grids[i][j] = 100
current_grids = initial_grids
inital_for_compare = np.array(copy.deepcopy(current_grids))
# define velocity
velocity_expression = Constant((0, 0))
# Creating the grid and function space
mesh = UnitSquareMesh(row, col) # Here the range is from 0 to 1
element = FiniteElement('P', triangle, 1)
V = FunctionSpace(mesh, element)
# Defining variables and testing functions
c = TrialFunction(V)
w = TestFunction(V)
# Define c_n as a Function of the solution at the previous time step
c_n = assign_initial_conditions(V, current_grids, row)
current_solution = return_solutions(V, c_n, row, col, row, 0, if_keep_decimal_point,16)
solution_matrix.append(current_solution)
rmse_values_limited = np.sqrt(np.mean((np.array(current_solution) - inital_for_compare) ** 2, axis=(0, 1)))
print(rmse_values_limited)
# define variation form
diffusion_term = r * dot(grad(c), grad(w))
convection_term = dot(velocity_expression, grad(c)) * w
source_term = f * w
F = ((c - c_n) / dt * w + convection_term + diffusion_term - source_term) * dx
a, L = lhs(F), rhs(F)
c = Function(V)
# for every time step solve
for n in range(10):
velocity_expression.n = n
solve(a == L, c)
# update previous step solution
c_n.assign(c)
current_solution = return_solutions(V, c_n, row, col, row, n+1, if_keep_decimal_point,16)
solution_matrix.append(current_solution)
rmse_values_limited = np.sqrt(np.mean((np.array(current_solution) - inital_for_compare) ** 2, axis=(0, 1)))
print(f'rmse between initial condition and solutions in time step {n+1} is: ', rmse_values_limited)
return solution_matrix
solutions = fem_return_solutions()
Run the code, we can get a slight deviation from the initial value:
