fipy problem with discontinuity in an 1D advection-diffusion system

Question

I am trying to model a 1-D advection-diffusion problem involving a variable advection velocity. It concerns a block of ice flowing downwards. I experience problems with an unphysical discontinuity occurring at the point in the grid where the velocity changes (plot attached). The code works fine for the case of constant velocity. Any help will be very much appreciated

Top plot the advection coeff as a function of depth Middle plot: the temperature profile of the cie slab as a function of depth. The boundary conditions are 273 and 223 for the two edges and the discontinuity appears at the point where the advection coefficient changes Bottom plot: The same problem solved for a constant advection velocity equal to -0.05 m/y = green dashed curve is the initial condition, blues curves show the evolution of the profile at various time steps and the red curve is the final solution for t = 1e5 years.

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)

Answer 1

I appreciate that the result seems counterintuitive. It might be instructive to look at the solution as it approaches the purely convective case. When D_ice = 1.13e-8 #m2s-1 (and increasing the elapsed time by a factor of ten),

we observe that the solution approaches temp_right on the right, but some much lower constant temperature on the left. This is because more heat is fluxing out of the domain on the left than is fluxing in on the right, precisely because of the Dirichlet boundary conditions. The flux in at the right is temp_right * conv_coeff[..., mesh.facesRight], which is about 11 m K / yr, whereas the flux out at the left is temp_left * conv_coeff[..., mesh.facesLeft], or 27 m K / yr. This discrepancy in heat flux has to come from somewhere, and it is satisfied by cooling down the domain.

With higher diffusivity, the system seeks a compromise between this purely convective solution and the purely diffusive solution of a diagonal line between temp_left and temp_right, leading to the solution you asked about.

Part of the issue in seeing this may be that what you've written appears to be a conservative equation for temperature

but temperature is not a conserved quantity. The more typical form of the heat equation

makes clear that heat is being conserved and that resulting temperature excursions are not unphysical.