Hi, FEniCS community,

I have been stuck and confused with what seems to be a very simple problem for nearly two months, and I have debugged from several views, even I have calculated by hands for a 1×1 grid case. But still got negative values under some system parameters. I have looked up some papers but have not found the proper answer yet, so I came here to help.

I have a dispersion process PDE, which is the fundamental one:

and I discretize this PDE by the standard Galerkin method.

I code this in fenics, but get the negative solution under some system parameters, which is not expected. The code is here:

```
# -*- coding: utf-8 -*-
# @Time : 05-08-2024 19:06
# @Author : snow
# @FileName: example_compare.py
# @Software: PyCharm
from fenics import *
import matplotlib.pyplot as plt
from common_util.common_parameter_function import theta_transform
from fem_helper import assign_initial_conditions
import random
import math
import numpy as np
random.seed(0)
def initial_condition_constant(row, col, center_initial_size, value):
initial_grids = [[0 for _ in range(col)] for _ in range(row)]
center_x, center_y = row // 2, col // 2 # 中心位置
half_size = center_initial_size // 2 # size的一半，用于确定中心区域的范围
if center_initial_size % 2 == 0: # 如果边长是偶数
for i in range(center_x - half_size, center_x + half_size):
for j in range(center_y - half_size, center_y + half_size):
initial_grids[i][j] = value
else: # 如果边长是奇数
for i in range(center_x - half_size, center_x + half_size + 1):
for j in range(center_y - half_size, center_y + half_size + 1):
initial_grids[i][j] = value
return initial_grids
class FowVelocityExpression2(UserExpression):
"""
https://fenicsproject.org/olddocs/dolfin/2017.2.0/python/programmers-reference/functions/expression/Expression.html
"""
def __init__(self, v_x_list, v_y_list, **kwargs):
self.v_x_list = v_x_list
self.v_y_list = v_y_list
self.current_step = 0
super().__init__(**kwargs)
def eval(self, values, x):
values[0] = self.v_x_list[self.current_step] # x方向的风速分量
values[1] = self.v_y_list[self.current_step] # y方向的风速分量
def value_shape(self):
return (2,)
def plot_node_fem(V, c_n, row, col, length, t, plt_same_scale, plt_vmin, plt_vmax):
"""
plot node value
"""
show_grid = [[0 for j in range(row)]for i in range(col)]
c_n_values = c_n.vector().get_local()
dof_coordinates = V.tabulate_dof_coordinates().reshape((-1, 2))
for i in range(len(dof_coordinates)):
x, y = dof_coordinates[i]
ix, iy = min(int(x * length), length), min(int(y * length), length)
ix_matrix = length - iy
iy_matrix = ix
show_grid[ix_matrix][iy_matrix] = c_n_values[i]
plt.figure()
if plt_same_scale:
plt.imshow(show_grid, cmap='viridis_r', vmin=plt_vmin, vmax=plt_vmax)
else:
plt.imshow(show_grid, cmap='viridis_r')
plt.colorbar(label='Concentration values')
plt.title(f'Solution $c(p,t)$ at time step={t}')
plt.show()
if __name__ == '__main__':
plt_same_scale = True
plt_vmin = -10
plt_vmax = 110
# system parameter
row = 32
col = 32
d = 0.002 #Diffusion coefficient
f = 0 # external_source
dt = 0.01
cell_length = 1 / row
current_grids = initial_condition_constant(row, col, 5, 100)
## Set flow velocity and directions
theta_v_dict = [(30, 0.5), (30, 0.5), (30, 0.5), (30, 0.5), (30, 0.5), (30, 0.5), (30, 0.5), (30, 0.5)]
theta_in_degrees_origin_list = []
theta_in_radians_list = []
v_list = []
v_x_list = []
v_y_list = []
for theta, v in theta_v_dict:
theta_in_degrees_origin_list.append(theta)
theta_in_radians_list.append(math.radians(theta_transform(theta)))
v_list.append(v)
v_x_list.append(v * np.cos(math.radians(theta)))
v_y_list.append(v * np.sin(math.radians(theta)))
# define velocity
velocity_expression = FowVelocityExpression2(v_x_list, v_y_list, degree=1)
# define grid
mesh = UnitSquareMesh(row - 1, row - 1) #这里区间是0到1的
element = FiniteElement('P', triangle, 1)
V = FunctionSpace(mesh, element)
# define varibal and test function
c = TrialFunction(V)
w = TestFunction(V)
# assign initial condition, which is a 32×32 grid, w
# initial value is 5×5 center region with value 100, others 0,
c_n = assign_initial_conditions(V, current_grids, row - 1)
# define the variational form
diffusion_term = d * 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
# F = c*w*dx + dt * d *dot(grad(c), grad(w))* dx - c_n * w *dx + dt * dot(velocity_expression, grad(c)) * w * dx
a, L = lhs(F), rhs(F)
c = Function(V)
# solve the equation
for n in range(len(v_list)):
velocity_expression.current_step = n
solve(a == L, c)
c_n.assign(c)
plot_node_fem(V, c_n, row, col, row - 1, n + 1, plt_same_scale, plt_vmin, plt_vmax)
```