FiPy模拟一维平流-扩散系统（冰体流动场景）时的非物理间断问题求助

我正在尝试用FiPy模拟冰体向下流动的一维平流-扩散问题，当平流速度在网格中发生突变时，求解结果在突变位置出现了非物理间断，但平流速度恒定的情况下代码运行完全正常，希望能得到大家的技术帮助。

结果图说明 #

• 顶部图：平流系数随深度的变化曲线
• 中间图：冰板温度剖面随深度的变化曲线，边界条件为两端273K和223K，非物理间断出现在平流系数突变处
• 底部图：平流速度恒定为-0.05 m/y时的求解结果，其中绿色虚线为初始条件，蓝色曲线为不同时间步的剖面演化过程，红色曲线为t=1e5年的最终解

模拟代码 #

import numpy as np
from scipy import interpolate, signal, ndimage, integrate
import openpyxl
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from matplotlib.ticker import MultipleLocator, FormatStrFormatter, FixedLocator, FixedFormatter
from matplotlib.backends.backend_pdf import PdfPages
import time
import pprint
import fipy

plt.rc("legend", fontsize = 16)
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['axes.labelsize'] = 16
plt.rcParams['axes.titlesize'] = 16
plt.rcParams['xtick.labelsize'] = 16
plt.rcParams['xtick.direction'] = "in"
plt.rcParams['ytick.labelsize'] = 16
plt.rcParams['ytick.direction'] = "in"
plt.rc('text.latex', preamble=r'\usepackage{cmbright}')
plt.rcParams["figure.autolayout"] = False
plt.ion()
plt.close("all")

H = 2764
dx = 10
dt = 100
t_total = 20000
plotit = True

llib_dict = {"H_llib": 2579, "a_llib": 0.016, "p_llib": 5.5}
lliboutry = False
sim_dict= {"H": H, "dx": dx, "dt": dt, "t_total": t_total, "llib_dict": llib_dict}

t1 = time.time()

if plotit:
    plt.close("all")
    f1, axes1 = plt.subplots(nrows = 1, ncols = 1, num = 5522, figsize = (6,6), tight_layout = True)
    axes1.set_ylabel(r'Temperature [K]')
    axes1.set_xlabel(r'z from bedrock(m)')
    f2, axes2 = plt.subplots(nrows = 1, ncols = 1, num = 5523, figsize = (6,6), tight_layout = True)
    axes2.set_ylabel(r'Advection [ma-1]')
    axes2.set_xlabel(r'z from bedrock (m)')

nx = H/dx
sim_dict["nx"] = nx
mesh = fipy.Grid1D(dx=dx, nx=nx)
sim_dict["mesh"] = mesh
X = mesh.faceCenters[0]
Xc = mesh.cellCenters[0]

one_yr = 365.25*24*3600
t_total_s = t_total*one_yr
sim_dict["t_total_s"] = t_total_s
dt_s = dt*one_yr
sim_dict["dt_s"] = dt_s
n_steps = t_total_s/dt_s
sim_dict["n_steps"] = n_steps

D_ice = 1.13e-6 #m2s-1
sim_dict["D_ice"] = D_ice
k_ice = 2.1 #Jm-1K-1s-1

if lliboutry == True:
    p = llib_dict["p_llib"]
    d = np.arange(0, llib_dict["H_llib"] + 1, 1)
    lamda_lliboutry = llib_dict["a_llib"]*(1-(p+2)/(p+1)*d/llib_dict["H_llib"] + 1/(p+1)*(d/llib_dict["H_llib"])**(p+2))
    conv_coeff_arr = fipy.FaceVariable(mesh = mesh, name = "conv_coeff", value = [-np.interp(mesh.faceCenters.numericValue[0], d, lamda_lliboutry, right = 1e-7)[::-1]])
else:
    conv_coeff_arr = np.zeros_like(mesh) - 0.05 #convection coeff
    conv_coeff = fipy.FaceVariable(mesh = mesh, value = [conv_coeff_arr/one_yr])
    conv_coeff.setValue(-1e-1/one_yr, where = (X<=H/2))

sim_dict["conv_coeff"] = conv_coeff

F = (D_ice*dt)/(dx**2)
sim_dict["F"] = F
sim_dict["Pf"] = D_ice/np.mean(conv_coeff)

print("\n")
print("Dice: %0.3e" %D_ice)
print("conv_coeff_mean: %0.3e" %np.mean(conv_coeff))
print(("F number: %0.3e" %F))

temp_left = 270.15
temp_right = 223.15
flux_left = 0

phi = fipy.CellVariable(mesh = mesh, name = "Temperature", value = 218.15)
phi.setValue(250.15)

eqX = fipy.TransientTerm() == fipy.DiffusionTerm(coeff=D_ice) -fipy.PowerLawConvectionTerm(coeff = conv_coeff)

phi.constrain(temp_left, where=mesh.facesLeft)
phi.constrain(temp_right, where=mesh.facesRight)

if plotit:
    axes1.plot(Xc.value, phi.value, linewidth = 0.9, color = "g", linestyle = ":")

for i in np.arange(n_steps+1):
    eqX.solve(var = phi, dt = dt_s, solver = fipy.LinearLUSolver(tolerance = 1.e-15))
    print("\t%i/%i steps - %i y" %(i, n_steps, i*dt_s/one_yr), end = "\r")
    if i%(n_steps/4)==0:
        if plotit:
            axes1.plot(Xc.value, phi.value, linewidth = 0.7, color = "b")

if plotit:
    axes1.plot(Xc.value, phi.value, linewidth = 0.7, color = "r")
    axes2.plot(X.value, conv_coeff[0].value*one_yr, linewidth = 0.8, color = "k")

print("\n")
exec_time = time.time() - t1
sim_dict["exec_time"] = exec_time
sim_dict["phi"] = phi
pprint.pprint(sim_dict)

内容的提问来源于stack exchange，提问作者VG_nbi