Inaccuracy when applying internal fixed gradient: evident leakage through zero-flux boundaries

[taoyun]

I implemented zero-flux internal boundary condition in 2d space using the official method, i.e., through the use of an implicit source term (https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#internal-fixed-value). In the simple diffusion example below, there exists an internal rectangular "barrier" of which the boundaries are reflective. However, it is clear from the resultant plots (short animation attached) that the solution "leaks" across these supposedly shielded internal boundaries. Is this a limitation of FiPy or might there be a workaround to this issue?

import numpy as np
from fipy import *
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt

plt.rc('text', usetex=True)

fig =plt.figure()
txt = fig.suptitle('place holder', fontsize=14)
ax = plt.subplot()
plt.show()

# generate mesh
l = 6.
nx = 301 
ny = nx
dx = l/nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]]

# initial condition
x_init, y_init = -2., 2.
x_var, y_var = 0.1, 0.1
distr = bivariate_normal(mesh.x, mesh.y, x_var, y_var, x_init, y_init)
phi = CellVariable(name="solution variable", mesh=mesh, value=distr)

# function to reshape solution for plotting
def data_process(arg):
    p1 = arg.value.reshape(mesh.shape)
    return p1

# baseline diffusion coefficient
d0 = 0.1

# positional variables
xy_cell = mesh.cellCenters()
xy_face = FaceVariable(mesh=mesh, value=mesh.faceCenters(), rank=1)

# internal zero-flux boundaries
xx, yy = xy_face
mask = ((xx > -1.7) & (xx < -1.5) & (yy > 1.7))
mask = FaceVariable(mesh, value=mask)
d = FaceVariable(mesh=mesh, value=d0)
d[mask.value] = 0.

largeValue = 1e+10

# main equation
eqX = (TransientTerm(var=phi) == DiffusionTerm(coeff=d, var=phi) + \
                                 DiffusionTerm(coeff=largeValue * mask, var=phi) - \
                                 ImplicitSourceTerm(coeff=(mask * largeValue * 0. * mesh.faceNormals).divergence, var=phi))

""" SIMULATION """
# time steps
dt = 0.004
steps = 100

# plot at initial time
p1 = data_process(phi)
txt.set_text(r'Time Step 0')
img = ax.imshow(p1, origin='lower', extent=[-l/2,l/2,-l/2,l/2])
cb = fig.colorbar(img)
plt.draw()

# solving transient solutions
for step in range(steps):
    print 'step', step
    eqX.solve(var=phi, dt=dt)
    
    p1 = data_process(phi)
    txt.set_text(r'Time Step %s'%(step+1))
    
    ax.cla()
    img = ax.imshow(p1, origin='lower', extent=[-l/2,l/2,-l/2,l/2])
    plt.draw()

internalBC_leak.mp4

[guyer]

There are a few issues here, although I don't ever observe the leakage you show here.

• You appear to be running pretty old versions of things. Updating your packages may resolve many or all of your issues:
  • Your print syntax is from Python 2.x
  • matplotlib.mlab.bivariate_normal() was removed four years ago.

Once those issues are resolved, there's no observable leakage into the barrier, but conservation is poor, with about 27% loss of mass.

The issue there is that the internal fixed gradient technique doesn't make any sense when the gradient is zero (the ImplicitSourceTerm adds zeroes to the matrix, which is the same thing as doing nothing). It is better to just zero the diffusivity in the barrier and leave off the rest, resulting in very good conservation.

The code I used is

import numpy as np
from fipy import *
from scipy.stats import multivariate_normal 
import matplotlib.pyplot as plt

plt.rc('text', usetex=True)

fig =plt.figure()
txt = fig.suptitle('place holder', fontsize=14)
ax = plt.subplot()
plt.show()

# generate mesh
l = 6.
nx = 301 
ny = nx
dx = l/nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]]

# initial condition
x_init, y_init = -2., 2.
x_var, y_var = 0.1, 0.1
# Replacement for matplotlib.mlab.multivariate_normal
# from https://stackoverflow.com/a/67595065/2019542
rv = multivariate_normal([x_init, y_init], [[x_var, 0], [0, y_var]])
distr = rv.pdf(np.dstack((mesh.x, mesh.y)))
phi = CellVariable(name="solution variable", mesh=mesh, value=distr)

# function to reshape solution for plotting
def data_process(arg):
    p1 = arg.value.reshape(mesh.shape)
    return p1

# baseline diffusion coefficient
d0 = 0.1

# positional variables
xy_cell = mesh.cellCenters()
xy_face = FaceVariable(mesh=mesh, value=mesh.faceCenters(), rank=1)

# internal zero-flux boundaries
xx, yy = xy_face
mask = ((xx > -1.7) & (xx < -1.5) & (yy > 1.7))
mask = FaceVariable(mesh, value=mask)
d = FaceVariable(mesh=mesh, value=d0)
d[mask.value] = 0.

phi[..., (mesh.x > -1.7) & (mesh.x < -1.5) & (mesh.y > 1.7)] = 0

largeValue = 1e+10

# main equation
eqX = (TransientTerm(var=phi) == DiffusionTerm(coeff=d, var=phi))

""" SIMULATION """
# time steps
dt = 0.004
steps = 100

viewer = Viewer(vars=phi, datamin=0, datamax=1.5)
viewer.plot()
print('step', 0, phi.cellVolumeAverage)

# solving transient solutions
for step in range(steps):
    eqX.solve(var=phi, dt=dt)
    
    viewer.plot()
    print('step', step, phi.cellVolumeAverage)

I used FiPy's standard viewer because your imshow() weren't displaying for me in Jupyter (matplotlib animation updates are really sensitive to versions and exactly what the display environment is). I also zeroed out the variable inside the barrier, but that's easily removed if undesirable.

I used FiPy 3.4.3, Python 3.10.6, and petsc4py 3.17.3.

[taoyun]

Hi Jon, Thank you for the extremely thorough explanation and carefully annotated sample code - I'll take a few days to mull over and experiment with all the details. I assume that this approach applies to a convection-diffusion model as well: there, I should just zero both the convection and diffusion coefficients at the barrier and not worry about implicit/explicit terms. Much appreciated! Yun

…

On Mon, Apr 24, 2023 at 1:43 PM Jonathan Guyer ***@***.***> wrote: There are a few issues here, although I don't ever observe the leakage you show here. - You appear to be running pretty old versions of things. Updating your packages may resolve many or all of your issues: - Your print syntax is from Python 2.x - matplotlib.mlab.bivariate_normal() was removed four years ago <https://github.com/matplotlib/matplotlib/blob/d77801ea94b06b990fc2b535a42522b45143e630/doc/api/prev_api_changes/api_changes_3.1.0.rst#L578> . Once those issues are resolved, there's no observable leakage into the barrier, but conservation is poor, with about 27% loss of mass. [image: image] <https://user-images.githubusercontent.com/1830725/234071452-bb7c7397-7bfc-4e0f-8921-46df9422116f.png> The issue there is that the technique from doesn't make any sense when the gradient is zero (the ImplicitSourceTerm adds zeroes to the matrix, which is the same thing as doing nothing). It is better to just zero the diffusivity in the barrier and leave off the rest, resulting in very good conservation. [image: image] <https://user-images.githubusercontent.com/1830725/234072455-7992662f-55e3-43e5-bb63-30a5c3dfc5b9.png> The code I used is import numpy as npfrom fipy import *from scipy.stats import multivariate_normal import matplotlib.pyplot as plt plt.rc('text', usetex=True) fig =plt.figure()txt = fig.suptitle('place holder', fontsize=14)ax = plt.subplot()plt.show() # generate meshl = 6.nx = 301 ny = nxdx = l/nxdy = dxmesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]] # initial conditionx_init, y_init = -2., 2.x_var, y_var = 0.1, 0.1# Replacement for matplotlib.mlab.multivariate_normal# from https://stackoverflow.com/a/67595065/2019542rv = multivariate_normal([x_init, y_init], [[x_var, 0], [0, y_var]])distr = rv.pdf(np.dstack((mesh.x, mesh.y)))phi = CellVariable(name="solution variable", mesh=mesh, value=distr) # function to reshape solution for plottingdef data_process(arg): p1 = arg.value.reshape(mesh.shape) return p1 # baseline diffusion coefficientd0 = 0.1 # positional variablesxy_cell = mesh.cellCenters()xy_face = FaceVariable(mesh=mesh, value=mesh.faceCenters(), rank=1) # internal zero-flux boundariesxx, yy = xy_facemask = ((xx > -1.7) & (xx < -1.5) & (yy > 1.7))mask = FaceVariable(mesh, value=mask)d = FaceVariable(mesh=mesh, value=d0)d[mask.value] = 0. phi[..., (mesh.x > -1.7) & (mesh.x < -1.5) & (mesh.y > 1.7)] = 0 largeValue = 1e+10 # main equationeqX = (TransientTerm(var=phi) == DiffusionTerm(coeff=d, var=phi)) """ SIMULATION """# time stepsdt = 0.004steps = 100 viewer = Viewer(vars=phi, datamin=0, datamax=1.5)viewer.plot()print('step', 0, phi.cellVolumeAverage) # solving transient solutionsfor step in range(steps): eqX.solve(var=phi, dt=dt) viewer.plot() print('step', step, phi.cellVolumeAverage) I used FiPy's standard viewer because your imshow() weren't displaying for me in Jupyter (matplotlib animation updates are really sensitive to versions and exactly what the display environment is). I also zeroed out the variable inside the barrier, but that's easily removed if undesirable. I used FiPy 3.4.3, Python 3.10.6, and petsc4py 3.17.3. — Reply to this email directly, view it on GitHub <#918 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMF556J3JNQ7VISWYJ6Z5LXC23UDANCNFSM6AAAAAAXH5S5SM> . You are receiving this because you authored the thread.Message ID: ***@***.***>

-- Yun Tao Research Faculty Institute of Bioinformatics, University of Georgia

[guyer]

I assume that this approach applies to a convection-diffusion
model as well: there, I should just zero both the convection and diffusion
coefficients at the barrier and not worry about implicit/explicit terms.

Correct