How to declare two dimensional boundaries using Fipy

Question

I am trying to reproduce Matlab example https://se.mathworks.com/help/matlab/math/solve-system-of-pdes.html using FiPy. However, I am unable to reproduce the Matlab results (I am newbie when using Fipy). I think the problem is probably in function or in boundary condition declarations. Also, I did not find Fipy setValue property for variable with respect of t for initial conditions, so I override it with very clumsy loop structure. Here is my test code:

from fipy import *
from fipy import Grid2D, CellVariable,TransientTerm, DiffusionTerm
from fipy.tools import numerix
import matplotlib.pyplot as plt
import numpy as np

from sys import stdout
from builtins import range

import plotly.graph_objects as go

nx = 100
ny = 100
Lx = 1
Ly = 2

mesh = Grid2D(nx=nx, ny=ny,Lx=Lx,Ly=Ly)
u1 = CellVariable(name = "solution variable u1", hasOld=True, mesh = mesh, value = 0.)
u2 = CellVariable(name = "solution variable u2", hasOld=True, mesh = mesh, value = 0.)

z1=np.zeros((nx,ny))
z2=np.zeros((nx,ny))

u_mask1=np.zeros(nx*ny)
u_mask2=np.zeros(nx*ny)

#initial conditions
#compute mask

z1[:,0]=1
z2[:,0]=0

for i in range(nx):      #x
    for j in range(ny):  #t
        if i==0:
            position=j
        else:
            position=i*ny+j-1
        u_mask1[position]=z1[i,j]
        u_mask2[position]=z2[i,j]

u1.setValue(1, where=u_mask1)
#u2.setValue(0, where=u_mask2) #unnecessary

#constrains
u1.constrain(1.0, where=mesh.facesRight) #Direction?
u2.constrain(0.0, where=mesh.facesLeft)  #Direction?
u1.grad.dot([1.0, 0.0]).constrain(0.0, where=mesh.facesLeft) 
u2.grad.dot([1.0, 0.0]).constrain(0.0, where=mesh.facesRight)


eqn0 = TransientTerm(var=u1) == DiffusionTerm(0.024, var=u1) - (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
eqn1 = TransientTerm(var=u2) == DiffusionTerm(0.170, var=u2) + (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
    
eqn = eqn0 & eqn1
#compute equations
max_range=1000

for t in range(0,max_range):
    stdout.write("\r%3.2f" % (t/max_range*100))
    stdout.flush()
    u1.updateOld()
    u2.updateOld()
    eqn.sweep(dt=0.002)
    #eqn.solve(dt=1.e-3)

z1=np.zeros((nx,ny))
z2=np.zeros((nx,ny))
for i in range(nx):      #x
    for j in range(ny):  #t
        if i==0:
            position=j
        else:
            position=i*ny+j-1
        z1[i,j]=u1[position]
        z2[i,j]=u2[position]



x, y = np.linspace(0, Lx, nx), np.linspace(0, Ly, ny)

fig = go.Figure(data=[go.Surface(z=z1, x=x, y=y)])
fig.update_layout(scene = dict(
                    xaxis_title='X',
                    yaxis_title='T',
                    zaxis_title='Z'),
                    title='u1',
                    width=600,
                    height=600,
                    margin=dict(r=20, b=10, l=10, t=10)
                  )
fig.show()

fig = go.Figure(data=[go.Surface(z=z2, x=x, y=y)])
fig.update_layout(scene = dict(
                    xaxis_title='X',
                    yaxis_title='T',
                    zaxis_title='Z'),
                    title='u2',
                    width=600,
                    height=600,
                    margin=dict(r=20, b=10, l=10, t=10)
                  )
fig.show()

Answer 1

I would implement in FiPy as:

from fipy import Grid1D, CellVariable,TransientTerm, DiffusionTerm, Viewer
from fipy.tools import numerix

nx = 100
Lx = 1.

mesh = Grid1D(nx=nx, Lx=Lx)
u1 = CellVariable(name = "solution variable u1", hasOld=True, mesh = mesh, value = 0.)
u2 = CellVariable(name = "solution variable u2", hasOld=True, mesh = mesh, value = 0.)

#initial conditions

u1.setValue(1)
#u2.setValue(0) #unnecessary

#boundary conditions

# u1.faceGrad.constrain([0.], where=mesh.facesLeft)  # no-flux is default BC
u2.constrain(0.0, where=mesh.facesLeft)  #Direction?
# u2.faceGrad.constrain([0.], where=mesh.facesLeft)  # no-flux is default BC
u1.constrain(1.0, where=mesh.facesRight) #Direction?


eqn0 = TransientTerm(var=u1) == DiffusionTerm(0.024, var=u1) - (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
eqn1 = TransientTerm(var=u2) == DiffusionTerm(0.170, var=u2) + (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
    
eqn = eqn0 & eqn1
#compute equations

viewer = Viewer(vars=(u1, u2), datamin=0., datamax=1.)
viewer.plot()

elapsed = 0.
dt = 0.002
while elapsed <= 2.:
    print(elapsed)
    u1.updateOld()
    u2.updateOld()
    eqn.sweep(dt=dt)
    viewer.plot()
    elapsed += dt

The main differences are:

• FiPy has no concept of defining a CellVariable in time. A CellVariable is defined in space and can evolve in time.
• I removed all of the plotly code. You can copy the 1D FiPy results at each timestep into an array and render them with matplotlib or plotly if you like.
• I put the initial and boundary conditions in the form FiPy is looking for.
  • Note: The default boundary condition for FiPy is no-flux; for this problem, that is equivalent to the normal gradient being zero.

As observed in the Matlab example, the initial time steps need to be much smaller than later on. You can use the steppyngstounes package to reduce the total number of sweeps by about a third and significantly improve the solution convergence at early times:

from steppyngstounes import PIDStepper

u1.updateOld()
u2.updateOld()

for step in PIDStepper(start=0., stop=2., size=0.0002):
    res = 1e100
    for sweep in range(3):
        res = eqn.sweep(dt=step.size)
    if step.succeeded(error=res / 0.01):
        u1.updateOld()
        u2.updateOld()
        viewer.plot()
    else:
        u1.value = u1.old
        u2.value = u2.old