Solve the 3D unsteady convection-diffusion equation

General
1

Hi, everyone!

This is what I wrote with fenicsx to solve the three-dimensional unsteady convection-diffusion equation, but the convergence order I obtained has always been incorrect.

But after I annotated the diffusion term, the model became a thermal equation, and the convergence order obtained was reasonable. Could you please help me figure out where the problem lies?

from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import mesh, fem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
import time 
from ufl import ds, dx, grad, inner, sin, cos, pi, exp, dot
# 精确解


def ExactRho(x,t):
    factor = np.exp(-t) + 1.0
    return factor * (x[0]**2 - x[0]) * (x[1]**2 - x[1]) * (x[2]**2 - x[2])

    

def  ExactU(x,t):
    factor = np.exp(-t) + 1.0
    temp1=factor * np.sin(np.pi * x[0]) * np.sin(np.pi * x[1]) * np.sin(np.pi * x[2])
    temp2=np.zeros_like(x[0])
    temp3=factor * (x[2]**2 - x[2]) * np.sin(np.pi * x[1]) * np.sin(np.pi * x[2])
    return [temp1,temp2,temp3]
    

def rhs_expr(x, t, D1):
    exp_term = np.exp(-t)
    factor = exp_term + 1.0
    return -exp_term*(x[0]**2 - x[0]) * (x[1]**2 - x[1]) * (x[2]**2 - x[2])  \
           -D1*factor*(2*((x[1]**2 - x[1]) * (x[2]**2 - x[2])+(x[0]**2 - x[0])* (x[2]**2 - x[2])+(x[0]**2 - x[0]) * (x[1]**2 - x[1])))\
           +factor**2*((x[1]**2 - x[1]) * (x[2]**2 - x[2])*((2*x[0]-1.0)*np.sin(np.pi*x[0])+np.pi*(x[0]**2 - x[0])*np.cos(np.pi*x[0]))*np.sin(np.pi*x[1])*np.sin(np.pi*x[2]))\
           +factor**2*((x[0]**2 - x[0]) * (x[1]**2 - x[1])*2*(x[2]**2 - x[2])*(2*x[2]-1)*np.sin(np.pi*x[0])*np.sin(np.pi*x[1]))
        


# 固定终止时间
T = 1.0
# 结果记录
hs = []
errors_L2 = []
errors_max = []  
start_time = time.time()
for i in range(3):
    n = 2 ** (i + 2)               # 空间划分数 n = 2^{i+2}
    h = 1.0 / n
    num_steps = 2 ** (2 * i + 4)   # 时间划分数 tau = 2^{2i+4}
    dt = T / num_steps
    D1 = dt**2
    t = 0.0
# 创建 3D 网格（单位立方体）
    domain = mesh.create_unit_cube(MPI.COMM_WORLD, n, n, n, mesh.CellType.tetrahedron)
    V = fem.functionspace(domain, ("Lagrange", 1))  # 线性 Lagrange 元
    Vv=fem.functionspace(domain, ("Lagrange", 1, (3, )))

    rho_n = fem.Function(V)
    rho_n.interpolate(lambda x: ExactRho(x, t))
    rho_D = fem.Function(V)
    rho_D.interpolate(lambda x: ExactRho(x, t))
    
    
    u_old = fem.Function(Vv)
    u_old.interpolate(lambda x: ExactU(x, t))

    f_fun = fem.Function(V)

# 施加 Dirichlet 边界条件（3D 边界）
    tdim = domain.topology.dim  
    fdim = tdim - 1  
    domain.topology.create_connectivity(fdim, tdim)
    boundary_facets = mesh.exterior_facet_indices(domain.topology)
    bc = fem.dirichletbc(rho_D, fem.locate_dofs_topological(V, fdim, boundary_facets))

# 变分形式（3D 形式相同，但 dx 现在是 3D 积分）
    rho, w = ufl.TrialFunction(V), ufl.TestFunction(V)
  #  a_form = (rho * w + dt * D1 * dot(grad(rho), grad(w))) * dx- dt* rho * dot(u_old, grad(w))* dx
    F1 = (rho * w + dt * D1 * dot(grad(rho), grad(w))) * dx- dt * rho * dot(u_old, grad(w))* dx
    F2 = rho_n * w* dx + dt *f_fun* w * dx
    a = fem.form(F1)
    L = fem.form(F2)

# 组装矩阵和向量
    A = assemble_matrix(a, bcs=[bc])
    A.assemble()
    b = create_vector(L)
    rho_h = fem.Function(V)

# 定义求解器（直接法：LU 分解）
    solver = PETSc.KSP().create(domain.comm)
    solver.setOperators(A)
    solver.setType(PETSc.KSP.Type.PREONLY)
    solver.getPC().setType(PETSc.PC.Type.LU)

# 时间步进
    for n_step in range(num_steps):
      t += dt
      rho_D.interpolate(lambda x: ExactRho(x, t))
      bc = fem.dirichletbc(rho_D, fem.locate_dofs_topological(V, fdim, boundary_facets))
      f_fun.interpolate(lambda x: rhs_expr(x, t, D1))
      # 更新矩阵
      A.zeroEntries()
      assemble_matrix(A, a, bcs=[bc])
      A.assemble()

      # 重新组装右端向量
      with b.localForm() as loc_b:
        loc_b.set(0)
      assemble_vector(b, L)

      # 施加边界条件
      apply_lifting(b, [a], [[bc]])
      b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
      set_bc(b, [bc])

      # 求解线性系统
      solver.solve(b, rho_h.x.petsc_vec)
      rho_h.x.scatter_forward()

      # 更新 u_n
      rho_n.x.array[:] = rho_h.x.array
      u_old.interpolate(lambda x: ExactU(x, t))
     

# 计算误差（在更高阶空间）
    V_ex = fem.functionspace(domain, ("Lagrange", 2))
    rho_ex = fem.Function(V_ex)
    rho_ex.interpolate(lambda x: ExactRho(x,T))

# L2 误差
    error_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((rho_h - rho_ex)**2 * ufl.dx)), op=MPI.SUM))

    if MPI.COMM_WORLD.rank == 0:
        hs.append(h)
        errors_L2.append(error_L2)
        print(f"[i={i}] n={n}, num_steps={num_steps}, dt={dt:.3e}, h={h:.3e}")
        print(f"        → L2 error = {error_L2:.3e}")

end_time = time.time()  # 计时结束
elapsed = end_time - start_time
# 输出收敛阶
if MPI.COMM_WORLD.rank == 0:
    print("\nEstimated convergence rates (L2):")
    for j in range(1, len(hs)):
        r = np.log(errors_L2[j-1]/errors_L2[j]) / np.log(hs[j-1]/hs[j])
        print(f"h={hs[j-1]:.3e} → h={hs[j]:.3e}, rate ≈ {r:.2f}")

    print(f"\nTotal runtime: {elapsed:.2f} seconds")
2

A few things I would consider:

  1. Write out the strong form (and the corresponding analytical solution) to ensure that it matches your code).
  2. Instead of deriving the solutions by hand, use the differentiation in UFL. This is for instance highlighted in: The standard way of compiling code with DOLFINx — FEniCS Workshop
  3. Visualize the exact solution and approximate solution side by side. What eyeball-norm differences do you spot?