Got unexpected negative values when solving the convection-diffusion PDE

variational formulation
1

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:

image
image928×118 8.38 KB

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

I don’t have time to read your lengthy example, but is your initial condition smooth? (I.e. does your initial condition or boundary conditions have a discontinuity?)

1 Like
3

Thank you very much for your reply.

My initial condition is a 32×32 grid, with the value of 100 in the central 5×5 region, while the rest of the grid has a value of 0.
Mathematically, this kind of initial condition is discontinuous.

I would appreciate the opportunity to discuss this in more depth.
For the finite element method, is it necessary to ensure that the initial conditions are continuous? Or, when the initial conditions are discontinuous, should we consider other discretization methods instead of the Galerkin method? The Galerkin method discretization variational form is as follow:

    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

If there are any related demos or resources that I could refer to, I would be very grateful. Thank you very much.

4

It looks like you’re using a continuous finite element basis

    # define grid
    mesh = UnitSquareMesh(row - 1, row - 1) #这里区间是0到1的
    element = FiniteElement('P', triangle, 1)
    V = FunctionSpace(mesh, element)

which cannot represent a discontinuous initial state. So I imagine the negative values you’re seeing are a result of the Gibbs phenomenon.

1 Like
5

Thank you very much for your reply. If I want to use discontinuous initial values(such as a 32×32 grid, with a value of 100 in the central 5×5 region, while the rest of the grid has a value of 0.), how should I set up the grid and the weak form of the convection-diffusion equation?

6

You can still use a continuous basis. You just have to be aware of the implications of the Gibbs phenomenon in your model. In your case that would be the negative values you see.

You could consider some kind of regularisation of your initial condition such that the discontinuous step is smoothed out. See for example Heaviside step function - Wikipedia.

If you’ve got plenty of time and willpower, you could consider a discontinuous Galerkin method. These methods are more resilient to the Gibbs phenomenon; however, not immune. This would require a more complicated variational form, and greater care with system setup, which may be out of the scope of your project.

7

Thank you for your reply. I will try the regularisation of the initial condition and let you know the result.

Additionally, if I apply a smoothing function to the initial condition, should I consider “unsmoothing” the solution afterward?

8

Hi Nate, the updated information: I have used some smooth methods, such as a blurring filter, which helps with the negative value magnitude, under the same very discontinous initial condition.

1 Like