

In [1]:

    
from sympy import init_session
init_session()









    



IPython console for SymPy 1.0 (Python 3.5.3-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/

Linear axisymmetric conservative interpolation

We want to construct a linear polynomal through that points $r_i$, $r_{i+1}$ that gives the correct averages, $r_i$, and $r_{i+1}$ when integrated over the volume, e.g.

$$\frac{1}{r_i \Delta r} \int_{r_{i-1/2}}^{r_{i+1/2}} f(r) r dr = f_i$$

There are 2 unknowns and two constraints

Define the polynomial



In [3]:

    
a, b = symbols("a b")
r, r0 = symbols("r r0")
f = a*(r-r0) + b
f









    Out[3]:





$$a \left(r - r_{0}\right) + b$$



In [4]:

    
dr = symbols("\Delta")

constraints

Define the constraint equations---here we set them up construct $A$ and $B$ integrals over the 2 control volumes



In [7]:

    
f0, fp = symbols("f_i f_{i+1}")

rm12 = r0 - Rational(1,2)*dr
rp12 = r0 + Rational(1,2)*dr
rp32 = r0 + Rational(3,2)*dr

r1 = r0 + dr

interfaces



In [6]:

    
rm12, rp12, rp32









    Out[6]:





$$\left ( - \frac{\Delta}{2} + r_{0}, \quad \frac{\Delta}{2} + r_{0}, \quad \frac{3 \Delta}{2} + r_{0}\right )$$

centers



In [8]:

    
r0, r1









    Out[8]:





$$\left ( r_{0}, \quad \Delta + r_{0}\right )$$



In [10]:

    
A = simplify(integrate(f*r/(r0*dr), (r, rm12, rp12)))
B = simplify(integrate(f*r/(r1*dr), (r, rp12, rp32)))

The analytic forms of the integrals



In [11]:

    
A, B









    Out[11]:





$$\left ( \frac{\Delta^{2} a}{12 r_{0}} + b, \quad \frac{1}{\Delta + r_{0}} \left(\frac{13 a}{12} \Delta^{2} + \Delta a r_{0} + \Delta b + b r_{0}\right)\right )$$

Our linear system is now:

$$A = f_i$$$$B = f_{i+1}$$

Now find the coefficients of the polynomial



In [12]:

    
coeffs = solve([A-f0, B-fp], [a,b])
coeffs









    Out[12]:





$$\left \{ a : \frac{12 r_{0} \left(\Delta + r_{0}\right) \left(f_{i} - f_{i+1}\right)}{\Delta \left(\Delta \left(\Delta + r_{0}\right) - r_{0} \left(13 \Delta + 12 r_{0}\right)\right)}, \quad b : \frac{\Delta f_{i+1} \left(\Delta + r_{0}\right) - f_{i} r_{0} \left(13 \Delta + 12 r_{0}\right)}{\Delta \left(\Delta + r_{0}\right) - r_{0} \left(13 \Delta + 12 r_{0}\right)}\right \}$$

And in pretty form, here's the polynomial



In [14]:

    
l = f.subs(a,coeffs[a]).subs(b,coeffs[b])
l









    Out[14]:





$$\frac{\Delta f_{i+1} \left(\Delta + r_{0}\right) - f_{i} r_{0} \left(13 \Delta + 12 r_{0}\right)}{\Delta \left(\Delta + r_{0}\right) - r_{0} \left(13 \Delta + 12 r_{0}\right)} + \frac{12 r_{0} \left(\Delta + r_{0}\right) \left(f_{i} - f_{i+1}\right) \left(r - r_{0}\right)}{\Delta \left(\Delta \left(\Delta + r_{0}\right) - r_{0} \left(13 \Delta + 12 r_{0}\right)\right)}$$

Now, evaluate it on the interface



In [16]:

    
a = l.subs(r, rp12)



In [17]:

    
simplify(a)









    Out[17]:





$$\frac{\Delta f_{i+1} \left(\Delta + r_{0}\right) - f_{i} r_{0} \left(13 \Delta + 12 r_{0}\right) + 6 r_{0} \left(\Delta + r_{0}\right) \left(f_{i} - f_{i+1}\right)}{\Delta \left(\Delta + r_{0}\right) - r_{0} \left(13 \Delta + 12 r_{0}\right)}$$



In [ ]: