Derivation of a symmetric stencil of
$$b = \nabla\cdot(A\nabla_\perp f)+Bf$$using a forward stencil on $\nabla\cdot(A\nabla_\perp f)$, and a backward stencil on $\nabla_\perp f$.
The stencil is made symmetric by mulitplying with $J(x,y)$. See ForwardsBackwardsNonSymmetric.ipynb for details.
In [1]:
from IPython.display import display
from sympy import init_printing
from sympy import symbols, expand, together, as_finite_diff, collect
from sympy import Function, Eq, Subs, Derivative
from collections import deque
init_printing()
In [2]:
def fromFunctionToGrid(expr, sym):
"""
Change from @(x,z) to @_xz, where @ represents a function
Input:
expr - The expression to change
sym - symbols('@_xz, @_xp1z, @_xm1z, @_xzp1, @_xzm1')
xp1 = x+hx
zm1 = z-hz
etc.
"""
curFun = str(syms[0]).split('_')[0]
for sym in syms:
curSuffix = str(sym).split('_')[1]
if curSuffix == 'xz':
expr = expr.subs(Function(curFun)(x,z), sym)
expr = expr.subs(Function(curFun)(x,y), sym)
elif curSuffix == 'xp1z':
expr = expr.subs(Subs(Function(curFun)(x,z), x, x+hx).doit(), sym)
expr = expr.subs(Subs(Function(curFun)(x,y), x, x+hx).doit(), sym)
elif curSuffix == 'xm1z':
expr = expr.subs(Subs(Function(curFun)(x,z), x, x-hx).doit(), sym)
expr = expr.subs(Subs(Function(curFun)(x,y), x, x-hx).doit(), sym)
elif curSuffix == 'xzp1':
expr = expr.subs(Subs(Function(curFun)(x,z), z, z+hz).doit(), sym)
expr = expr.subs(Subs(Function(curFun)(x,y), z, z+hz).doit(), sym)
elif curSuffix == 'xzm1':
expr = expr.subs(Subs(Function(curFun)(x,z), z, z-hz).doit(), sym)
expr = expr.subs(Subs(Function(curFun)(x,y), z, z-hz).doit(), sym)
return expr
In [3]:
x, y, z, hx, hz = symbols('x, y, z, h_x, h_z')
hx, hz = symbols('h_x, h_z', positive=True)
f = Function('f')(x, z)
A = Function('A')(x, z)
B = Function('B')(x, z)
gxx = Function('g^x^x')(x, y)
gzz = Function('g^z^z')(x, y)
J = Function('J')(x, y)
# Dummy function
g = Function('g')(x,z)
# Stencils
backwardX = [x-hx, x]
forwardX = [x, x+hx]
backwardZ = [z-hz, z]
forwardZ = [z, z+hz]
We are here discretizing the equation
$$ b = \nabla\cdot(A\nabla_\perp f)+Bf \simeq \frac{1}{J}\partial_x \left(JAg^{xx}\partial_x f\right) + \frac{1}{J}\partial_z \left(JAg^{zz}\partial_z f\right) + Bf$$where the derivatives in $y$ has been assumed small in non-orthogonal grids.
We will let $T$ denote "term", the superscript $^F$ denote a forward stencil, and the superscript $^B$ denote a backward stencil.
NOTE:
sympy
has a built in function as_finite_diff
, which could do the derivation easy for us. However it fails if
Subs
object (for example unevaluated derivatives calculated at a point)We therefore do this in a sligthly tedious way.
In [14]:
fx = f.diff(x)
fxB = as_finite_diff(fx, forwardX)
display(Eq(symbols('f_x'), fx))
display(Eq(symbols('f_x^F'), together(fxB)))
In [5]:
term1 = as_finite_diff(Derivative(J*A*gxx*g, x), backwardX)
display(Eq(symbols('T_1^F'), term1))
We now back substitute $g\to \partial_x f$
In [6]:
term1 = term1.subs(Subs(g,x,x-hx).doit(), Subs(fxB,x,x-hx).doit())
term1 = term1.subs(g, fxB)
display(Eq(symbols('T_1^F'), term1))
In [7]:
fz = f.diff(z)
fzB = as_finite_diff(fz, forwardZ)
display(Eq(symbols('f_z'), fz))
display(Eq(symbols('f_z^B'), together(fzB)))
In [8]:
term2 = as_finite_diff(Derivative(J*A*gzz*g, z), backwardZ)
display(Eq(symbols('T_2^F'), term2))
In [9]:
term2 = term2.subs(Subs(g,z,z-hz).doit(), Subs(fzB,z,z-hz).doit())
term2 = term2.subs(g, fzB)
display(Eq(symbols('T_2'), term2))
In [10]:
term3 = J*B*f
display(Eq(symbols('T_3^F'), term3))
In [11]:
Jb = term1 + term2 + term3
display(Eq(symbols('Jb'), Jb))
In [12]:
# Converting to grid syntax
functions = ['f', 'A', 'J', 'g^x^x', 'g^z^z', 'B']
for func in functions:
curStr = '{0}_xz, {0}_xp1z, {0}_xm1z, {0}_xzp1, {0}_xzm1'.format(func)
syms = symbols(curStr)
Jb = fromFunctionToGrid(Jb, syms)
In [13]:
# We must expand before we collect
Jb = collect(expand(Jb), symbols('f_xz, f_xp1z, f_xm1z, f_xzp1, f_xzm1'), exact=True)
display(Eq(symbols('Jb'),Jb))