

In [2]:

    
def gradVec2(u_vec, x, y):
    return sp.Matrix([[diff(u_vec[0], x), diff(u_vec[1],x)], [diff(u_vec[0], y), diff(u_vec[1], y)]])

def divTen2(tensor, x, y):
    return sp.Matrix([diff(tensor[0,0], x) + diff(tensor[1,0], y), diff(tensor[0, 1], x) + diff(tensor[1,1], y)])

def divVec2(u_vec, x, y):
    return diff(u_vec[0], x) + diff(u_vec[1], y)

2/1/17

FEProblemBase::reinitMaterials only calls property computation for Material objects currently active on the given subdomain so that's good. However, it's possible that material objects "active" on the subdomain aren't actually being used in any computing objects like kernels, etc. So we would like to do some additional checking.

Alright, let's say we're computing the residual thread. Then assuming we cannot compute properties in a material in isolation, we would like to do the next best thing: only call computeQpProperties for materials that have actually been asked to supply properties to kernels, dg_kernels, boundary_conditions, and interface_kernels.

So what am I doing as of commit d7dbfe5? I am determining the needed_mat_props through ComputeResidualThread::subdomainChanged() -> FEProblemBase::prepareMaterials. In the latter method, we first ask all materials--if there are any materials active on the block--to update their material property dependencies and then we ask materials on the boundaries of the subdomain to also update their dependencies. Note that this could lead to a boundary material object getting asked to update their material property dependencies twice because we first pass to MaterialWarehouse::updateMatPropDependenceyHelper all material objects as long as there is any one material object active in a block sense on the subdomain, and then we pass active material boundary objects. But this overlap doesn't matter so much because our needed_mat_props is a set, so if we try to insert the same material properties multiple times, it will silently and correctly fail. Note, however, that this could also pass needed_mat_props from material objects not on the current block, so that needs to be changed.

So what happens in MaterialWarehouse::updateMatPropDependencyHelper? We add mp_deps from MaterialPropertyInterface::getMatPropDependencies. However, it should be noted that this is only done for material objects. Is this fine? Well let's figure it out. It returns _material_property_dependencies which is a set added to by calling addMatPropDependency. Now this gets called when the object that inherits from MaterialPropertyInterface calls its own getMaterialProperty method. So I hypothesize that in my simple test, if I ask to have a material property in my kernel object with getMaterialProperty that will not register in any material objects list of _material_property_dependencies and consequently computeQpProperties will never get called. I will test that the next time I sit down at my comp.

2/2/17

Three tests:

Run with one material that doesn't supply any properties. Desired behavior: computeQpProperties does not get called. Expected to Pass. With devel MOOSE: expected to Fail (expected change)
Run two materials, one that supplies properties, another that does not. Desired behavior: computeQpProperties does not get called for the material not supplying properties while the other one does. Expected behavior: both materials compute methods get called. Fail. With devel MOOSE: expected to Fail (expected to not change)
Run with a kernel that uses a material property and an elemental aux variable that does not. Desired behavior: computeQpProperties should get called through the residual and jacobian threads but not through the aux kernel thread. Expected to Pass. With devel MOOSEE: expected to Fail (expected change)

Calls to computeProperties:

ComputeResidualThread
ComputeResidualThread 0th nonlinear residual printed
ComputeJacobianThread
ComputeResidualThread 0th linear residual printed
ComputeResidualThread 1st linear residual printed
ComputeResidualThread
ComputeResidualThread 1st nonlinear residual printed
ComputeElemAuxVarsThread -> Actually this is fine because this is the Aux Kernel that is created for outputting the material property

Number of calls: 8

1-4
5-8 ...
25-28
29-32

Failed tests:

random.material_serial
controls*

Failed but now passing:

element_aux_boundary
bnd_material_test
elem_aux_bc_on_bound
output.boundary
multiplicity
material_point_source_test
line_material_sampler

2/6/17

Ok my new test is failing with threads and I don't really know why. It seems like the number of calls to computing threads should be the same...

Calls to computeProperties:

ComputeResidualThread
ComputeResidualThread 0th nonlinear residual printed
ComputeJacobianThread
ComputeResidualThread 0th linear residual printed
ComputResidualThread 1st linear residual printed
ComputeResidualThread
ComputeResidualThread 1st nonlinear residual printed
ComputeElemAuxVarsThread

Yep so thread computing pattern is the exact same. How about whether the material is the same location in memory every time?

0x7fed90 (1, 2, 8) Increments: 1-4, 5-8, 9-12 -> average of 10.5 which is what is observed in the output
0x810b10 (3, 4, 5, 6, 7

4/28/17

Navier Stokes module development



In [2]:

    
import sympy as sp

sxx, sxy, syx, syy, nx, ny = sp.var('sxx sxy syx syy nx ny')



In [5]:

    
s = sp.Matrix([[sxx, sxy],[syx, syy]])
n = sp.Matrix([nx, ny])



In [6]:

    
s*n









    Out[6]:





Matrix([
[nx*sxx + ny*sxy],
[nx*syx + ny*syy]])



In [14]:

    
prod = n.transpose()*s*n
prod2 = n.transpose()*(s*n)



In [15]:

    
print(prod)
print(prod2)
print(prod==prod2)









    



Matrix([[nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy)]])
Matrix([[nx*(nx*sxx + ny*sxy) + ny*(nx*syx + ny*syy)]])
False



In [13]:

    
prod.shape









    Out[13]:





(1, 1)



In [18]:

    
sp.expand(prod) == sp.expand(prod2)









    Out[18]:





True



In [20]:

    
lhs = n.transpose()*s
print(lhs.shape)









    



(1, 2)



In [21]:

    
rhs = (n.transpose() * s * n) * n.transpose()
print(rhs.shape)









    



(1, 2)



In [62]:

    
rhs2 = (n.transpose()*s) * (n*n.transpose())
print(rhs2)
rhs3 = n.transpose() * (s*n*n.transpose())
print(sp.expand(rhs) == sp.expand(rhs2) == sp.expand(rhs3))









    



Matrix([[nx**2*(nx*sxx + ny*syx) + nx*ny*(nx*sxy + ny*syy), nx*ny*(nx*sxx + ny*syx) + ny**2*(nx*sxy + ny*syy)]])
True



In [58]:

    
print(n*n.transpose())









    



Matrix([[nx**2, nx*ny], [nx*ny, ny**2]])



In [59]:

    
print(n.transpose()*n)









    



Matrix([[nx**2 + ny**2]])



In [48]:

    
print(sp.simplify(lhs))
print(sp.simplify(rhs))









    



Matrix([[nx*sxx + ny*syx, nx*sxy + ny*syy]])
Matrix([[nx*(nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy)), ny*(nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy))]])



In [27]:

    
elml = lhs[0,0]
elmr = rhs[0,0]



In [29]:

    
print(elml.expand())
print(elmr.expand())









    



nx*sxx + ny*syx
nx**3*sxx + nx**2*ny*sxy + nx**2*ny*syx + nx*ny**2*syy



In [43]:

    
elmr.expand()









    Out[43]:





nx**3*sxx + nx**2*ny*sxy + nx**2*ny*syx + nx*ny**2*syy



In [42]:

    
elmr.expand().subs(nx, sp.sqrt(1 - ny**2))









    Out[42]:





ny**2*syy*sqrt(-ny**2 + 1) + ny*sxy*(-ny**2 + 1) + ny*syx*(-ny**2 + 1) + sxx*(-ny**2 + 1)**(3/2)



In [44]:

    
elmr.expand().subs(nx, sp.sqrt(1 - ny**2)).simplify()









    Out[44]:





ny**2*syy*sqrt(-ny**2 + 1) - ny*sxy*(ny**2 - 1) - ny*syx*(ny**2 - 1) + sxx*(-ny**2 + 1)**(3/2)



In [36]:

    
help(expr.replace)









    



Help on method replace in module sympy.core.basic:

replace(query, value, map=False, simultaneous=True, exact=False) method of sympy.core.add.Add instance
    Replace matching subexpressions of ``self`` with ``value``.
    
    If ``map = True`` then also return the mapping {old: new} where ``old``
    was a sub-expression found with query and ``new`` is the replacement
    value for it. If the expression itself doesn't match the query, then
    the returned value will be ``self.xreplace(map)`` otherwise it should
    be ``self.subs(ordered(map.items()))``.
    
    Traverses an expression tree and performs replacement of matching
    subexpressions from the bottom to the top of the tree. The default
    approach is to do the replacement in a simultaneous fashion so
    changes made are targeted only once. If this is not desired or causes
    problems, ``simultaneous`` can be set to False. In addition, if an
    expression containing more than one Wild symbol is being used to match
    subexpressions and  the ``exact`` flag is True, then the match will only
    succeed if non-zero values are received for each Wild that appears in
    the match pattern.
    
    The list of possible combinations of queries and replacement values
    is listed below:
    
    Examples
    ========
    
    Initial setup
    
        >>> from sympy import log, sin, cos, tan, Wild, Mul, Add
        >>> from sympy.abc import x, y
        >>> f = log(sin(x)) + tan(sin(x**2))
    
    1.1. type -> type
        obj.replace(type, newtype)
    
        When object of type ``type`` is found, replace it with the
        result of passing its argument(s) to ``newtype``.
    
        >>> f.replace(sin, cos)
        log(cos(x)) + tan(cos(x**2))
        >>> sin(x).replace(sin, cos, map=True)
        (cos(x), {sin(x): cos(x)})
        >>> (x*y).replace(Mul, Add)
        x + y
    
    1.2. type -> func
        obj.replace(type, func)
    
        When object of type ``type`` is found, apply ``func`` to its
        argument(s). ``func`` must be written to handle the number
        of arguments of ``type``.
    
        >>> f.replace(sin, lambda arg: sin(2*arg))
        log(sin(2*x)) + tan(sin(2*x**2))
        >>> (x*y).replace(Mul, lambda *args: sin(2*Mul(*args)))
        sin(2*x*y)
    
    2.1. pattern -> expr
        obj.replace(pattern(wild), expr(wild))
    
        Replace subexpressions matching ``pattern`` with the expression
        written in terms of the Wild symbols in ``pattern``.
    
        >>> a = Wild('a')
        >>> f.replace(sin(a), tan(a))
        log(tan(x)) + tan(tan(x**2))
        >>> f.replace(sin(a), tan(a/2))
        log(tan(x/2)) + tan(tan(x**2/2))
        >>> f.replace(sin(a), a)
        log(x) + tan(x**2)
        >>> (x*y).replace(a*x, a)
        y
    
        When the default value of False is used with patterns that have
        more than one Wild symbol, non-intuitive results may be obtained:
    
        >>> b = Wild('b')
        >>> (2*x).replace(a*x + b, b - a)
        2/x
    
        For this reason, the ``exact`` option can be used to make the
        replacement only when the match gives non-zero values for all
        Wild symbols:
    
        >>> (2*x + y).replace(a*x + b, b - a, exact=True)
        y - 2
        >>> (2*x).replace(a*x + b, b - a, exact=True)
        2*x
    
    2.2. pattern -> func
        obj.replace(pattern(wild), lambda wild: expr(wild))
    
        All behavior is the same as in 2.1 but now a function in terms of
        pattern variables is used rather than an expression:
    
        >>> f.replace(sin(a), lambda a: sin(2*a))
        log(sin(2*x)) + tan(sin(2*x**2))
    
    3.1. func -> func
        obj.replace(filter, func)
    
        Replace subexpression ``e`` with ``func(e)`` if ``filter(e)``
        is True.
    
        >>> g = 2*sin(x**3)
        >>> g.replace(lambda expr: expr.is_Number, lambda expr: expr**2)
        4*sin(x**9)
    
    The expression itself is also targeted by the query but is done in
    such a fashion that changes are not made twice.
    
        >>> e = x*(x*y + 1)
        >>> e.replace(lambda x: x.is_Mul, lambda x: 2*x)
        2*x*(2*x*y + 1)
    
    See Also
    ========
    subs: substitution of subexpressions as defined by the objects
          themselves.
    xreplace: exact node replacement in expr tree; also capable of
              using matching rules



In [49]:

    
t = lhs - rhs
print(t)









    



Matrix([[nx*sxx - nx*(nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy)) + ny*syx, nx*sxy + ny*syy - ny*(nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy))]])



In [51]:

    
t1 = t[0,0]
t2 = t[0,1]



In [52]:

    
print(t1)
print(t2)









    



nx*sxx - nx*(nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy)) + ny*syx
nx*sxy + ny*syy - ny*(nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy))



In [53]:

    
t1.simplify()









    Out[53]:





nx*sxx - nx*(nx*(nx*sxx + ny*syx) + ny*(nx*sxy + ny*syy)) + ny*syx



In [63]:

    
ddx, ddy, ux, uy = sp.var('ddx ddy ux uy')

grad = sp.Matrix([ddx,ddy])
u = sp.Matrix([ux,uy])
print(grad.shape)

phij,mu = sp.var('phij mu')
uDuxj = sp.Matrix([phij,0])

uDuyj = sp.Matrix([0,phij])



In [65]:

    
grad*u.transpose()









    Out[65]:





Matrix([
[ddx*ux, ddx*uy],
[ddy*ux, ddy*uy]])



In [78]:

    
jacx = n.transpose() * (mu * (grad*uDuxj.transpose() + (grad*uDuxj.transpose()).transpose())) * n
print(jacx)









    



Matrix([[ddy*mu*nx*ny*phij + nx*(2*ddx*mu*nx*phij + ddy*mu*ny*phij)]])



In [79]:

    
sp.expand(jacx[0,0])*nx









    Out[79]:





nx*(2*ddx*mu*nx**2*phij + 2*ddy*mu*nx*ny*phij)



In [81]:

    
jacy = n.transpose() * (mu * (grad*uDuyj.transpose() + (grad*uDuyj.transpose()).transpose())) * n
print(jacy)









    



Matrix([[ddx*mu*nx*ny*phij + ny*(ddx*mu*nx*phij + 2*ddy*mu*ny*phij)]])



In [82]:

    
sp.expand(jacy[0,0])*ny









    Out[82]:





ny*(2*ddx*mu*nx*ny*phij + 2*ddy*mu*ny**2*phij)



In [83]:

    
sp.factor(jacy[0,0])









    Out[83]:





2*mu*ny*phij*(ddx*nx + ddy*ny)



In [87]:

    
print(sp.factor((jacx[0,0]*n.transpose())[0,0]))
print(sp.factor((jacy[0,0]*n.transpose())[0,1]))









    



2*mu*nx**2*phij*(ddx*nx + ddy*ny)
2*mu*ny**2*phij*(ddx*nx + ddy*ny)



In [88]:

    
sJacX = mu * (grad*uDuxj.transpose() + (grad*uDuxj.transpose()).transpose())
sJacY = mu * (grad*uDuyj.transpose() + (grad*uDuyj.transpose()).transpose())
print(sJacX)
print(sJacY)









    



Matrix([[2*ddx*mu*phij, ddy*mu*phij], [ddy*mu*phij, 0]])
Matrix([[0, ddx*mu*phij], [ddx*mu*phij, 2*ddy*mu*phij]])



In [90]:

    
print(sp.factor((n.transpose()*sJacX)[0,0]))
print(sp.factor((n.transpose()*sJacY)[0,1]))









    



mu*phij*(2*ddx*nx + ddy*ny)
mu*phij*(ddx*nx + 2*ddy*ny)



In [85]:

    
jacx.shape









    Out[85]:





(1, 1)

Note that the total stress tensor ($\hat{\sigma}$) is equal to the sum of the deviatoric stress tensor ($\hat{\tau}$) and the stress induced by pressure ($-p\hat{I}$), e.g.

\begin{equation} \hat{\sigma} = \hat{\tau} - p\hat{I} \end{equation}



In [5]:

    
import sympy as sp

sxx, sxy, syx, syy, nx, ny, mu = sp.var('sxx sxy syx syy nx ny mu')
ddx, ddy, ux, uy = sp.var('ddx ddy ux uy')

grad = sp.Matrix([ddx,ddy])
u = sp.Matrix([ux,uy])

phij,mu = sp.var('phij mu')
uDuxj = sp.Matrix([phij,0])
uDuyj = sp.Matrix([0,phij])

rateOfStrain = (grad*u.transpose() + (grad*u.transpose()).transpose()) * 1 / 2
d_rateOfStrain_d_uxj = (grad*uDuxj.transpose() + (grad*uDuxj.transpose()).transpose()) * 1 / 2
d_rateOfStrain_d_uyj = (grad*uDuyj.transpose() + (grad*uDuyj.transpose()).transpose()) * 1 / 2
print(rateOfStrain)
print(d_rateOfStrain_d_uxj)
print(d_rateOfStrain_d_uyj)









    



Matrix([[ddx*ux, ddx*uy/2 + ddy*ux/2], [ddx*uy/2 + ddy*ux/2, ddy*uy]])
Matrix([[ddx*phij, ddy*phij/2], [ddy*phij/2, 0]])
Matrix([[0, ddx*phij/2], [ddx*phij/2, ddy*phij]])



In [7]:

    
tau = rateOfStrain * 2 * mu
d_tau_d_uxj = d_rateOfStrain_d_uxj * 2 * mu
d_tau_d_uyj = d_rateOfStrain_d_uyj * 2 * mu
print(tau)
print(d_tau_d_uxj)
print(d_tau_d_uyj)









    



Matrix([[2*ddx*mu*ux, mu*(ddx*uy + ddy*ux)], [mu*(ddx*uy + ddy*ux), 2*ddy*mu*uy]])
Matrix([[2*ddx*mu*phij, ddy*mu*phij], [ddy*mu*phij, 0]])
Matrix([[0, ddx*mu*phij], [ddx*mu*phij, 2*ddy*mu*phij]])



In [10]:

    
normals = sp.Matrix([nx,ny])
y_component_normal = sp.Matrix([0,ny])
x_component_normal = sp.Matrix([nx,0])
test = sp.var('test')
test_x = sp.Matrix([test,0])
test_y = sp.Matrix([0,test])

This is an example of an off-diagonal jacobian computation: derivative with respect to $x$ while test function corresponds to $y$

Specifically this corresponds to an off-diagonal contribution corresponding to the residual term:

\begin{equation} \vec{n}^T \cdot \hat{\tau} \cdot \vec{v}_y \end{equation}



In [12]:

    
normals.transpose() * d_tau_d_uxj * test_y









    Out[12]:





Matrix([[ddy*mu*nx*phij*test]])

Now let's look at an off diagonal-term for:

\begin{equation} \left(\vec{n}^T \cdot \hat{\tau} \cdot \vec{n} \right) \vec{n}^T \cdot \vec{v}_y \end{equation}



In [14]:

    
sp.factor(normals.transpose() * d_tau_d_uxj * normals * normals.transpose() * test_y)









    Out[14]:





Matrix([[2*mu*nx*ny*phij*test*(ddx*nx + ddy*ny)]])

Hmm...that's not very revealing...this result is completely symmetric...it doesn't tell me what the code implementation should be. Let's try 3D in order to elucidate



In [3]:

    
import sympy as sp

nx, ny, nz, mu, phij, ddx, ddy, ddz, ux, uy, uz = sp.var('nx ny nz mu phij ddx ddy ddz ux uy uz')
grad = sp.Matrix([ddx,ddy,ddz])
u = sp.Matrix([ux, uy, uz])
uDuxj = sp.Matrix([phij,0,0])
uDuyj = sp.Matrix([0,phij,0])
uDuzj = sp.Matrix([0,0,phij])

rateOfStrain = (grad*u.transpose() + (grad*u.transpose()).transpose()) * 1 / 2
d_rateOfStrain_d_uxj = (grad*uDuxj.transpose() + (grad*uDuxj.transpose()).transpose()) * 1 / 2
d_rateOfStrain_d_uyj = (grad*uDuyj.transpose() + (grad*uDuyj.transpose()).transpose()) * 1 / 2
d_rateOfStrain_d_uzj = (grad*uDuzj.transpose() + (grad*uDuzj.transpose()).transpose()) * 1 / 2
print(rateOfStrain)
print(d_rateOfStrain_d_uxj)
print(d_rateOfStrain_d_uyj)
print(d_rateOfStrain_d_uzj)









    



Matrix([[ddx*ux, ddx*uy/2 + ddy*ux/2, ddx*uz/2 + ddz*ux/2], [ddx*uy/2 + ddy*ux/2, ddy*uy, ddy*uz/2 + ddz*uy/2], [ddx*uz/2 + ddz*ux/2, ddy*uz/2 + ddz*uy/2, ddz*uz]])
Matrix([[ddx*phij, ddy*phij/2, ddz*phij/2], [ddy*phij/2, 0, 0], [ddz*phij/2, 0, 0]])
Matrix([[0, ddx*phij/2, 0], [ddx*phij/2, ddy*phij, ddz*phij/2], [0, ddz*phij/2, 0]])
Matrix([[0, 0, ddx*phij/2], [0, 0, ddy*phij/2], [ddx*phij/2, ddy*phij/2, ddz*phij]])



In [4]:

    
tau = rateOfStrain * 2 * mu
d_tau_d_uxj = d_rateOfStrain_d_uxj * 2 * mu
d_tau_d_uyj = d_rateOfStrain_d_uyj * 2 * mu
d_tau_d_uzj = d_rateOfStrain_d_uzj * 2 * mu
print(tau)
print(d_tau_d_uxj)
print(d_tau_d_uyj)
print(d_tau_d_uzj)









    



Matrix([[2*ddx*mu*ux, mu*(ddx*uy + ddy*ux), mu*(ddx*uz + ddz*ux)], [mu*(ddx*uy + ddy*ux), 2*ddy*mu*uy, mu*(ddy*uz + ddz*uy)], [mu*(ddx*uz + ddz*ux), mu*(ddy*uz + ddz*uy), 2*ddz*mu*uz]])
Matrix([[2*ddx*mu*phij, ddy*mu*phij, ddz*mu*phij], [ddy*mu*phij, 0, 0], [ddz*mu*phij, 0, 0]])
Matrix([[0, ddx*mu*phij, 0], [ddx*mu*phij, 2*ddy*mu*phij, ddz*mu*phij], [0, ddz*mu*phij, 0]])
Matrix([[0, 0, ddx*mu*phij], [0, 0, ddy*mu*phij], [ddx*mu*phij, ddy*mu*phij, 2*ddz*mu*phij]])



In [5]:

    
normals = sp.Matrix([nx,ny,nz])
test = sp.var('test')
test_x = sp.Matrix([test,0,0])
test_y = sp.Matrix([0,test,0])
test_z = sp.Matrix([0,0,test])



In [21]:

    
sp.factor(normals.transpose() * d_tau_d_uxj * normals * normals.transpose() * test_y)









    Out[21]:





Matrix([[2*mu*nx*ny*phij*test*(ddx*nx + ddy*ny + ddz*nz)]])



In [22]:

    
sp.factor(normals.transpose() * d_tau_d_uxj * normals * normals.transpose() * test_z)









    Out[22]:





Matrix([[2*mu*nx*nz*phij*test*(ddx*nx + ddy*ny + ddz*nz)]])



In [23]:

    
sp.factor(normals.transpose() * d_tau_d_uyj * normals * normals.transpose() * test_x)









    Out[23]:





Matrix([[2*mu*nx*ny*phij*test*(ddx*nx + ddy*ny + ddz*nz)]])

Alright, it looks like we get the normal components corresponding to residual $i$ and derivative variable $j$!!! Boom!



In [7]:

    
(normals.transpose() * tau)[0]









    Out[7]:





2*ddx*mu*nx*ux + mu*ny*(ddx*uy + ddy*ux) + mu*nz*(ddx*uz + ddz*ux)



In [9]:

    
sp.factor(_)









    Out[9]:





mu*(2*ddx*nx*ux + ddx*ny*uy + ddx*nz*uz + ddy*ny*ux + ddz*nz*ux)

5/3/17



In [30]:

    
from scipy.special import erf
from numpy import exp, sqrt, pi
import numpy as np

def u(x, y, u1, u2, sigma):
    return (u1 + u2) / 2. - (u1 - u2) / 2. * erf(sigma * y / x)

def v(x, y, u1, u2, sigma):
    return (u1 - u2) / (2. * sigma * sqrt(pi)) * exp(-(sigma * y / x)**2)

def p():
    return 0

def k(x, y, k0, sigma):
    return k0 * exp(-(sigma * y / x)**2)

def epsilon(x, y, epsilon0, sigma):
    return epsilon0 / x * exp(-(sigma * y / x)**2)

def muT(x, y, muT0, sigma):
    return muT0 * x * exp(-(sigma * y / x)**2)

def k0(u1, u2, sigma):
    return 343. / 75000. * u1 * (u1 - u2) * sigma / sqrt(pi)

def epsilon0(u1, u2, sigma, Cmu):
    return 343. / 22500. * Cmu * u1 * (u1 - u2)**2 * sigma**2 / pi

def muT0(u1, rho):
    return 343. / 250000. * rho * u1

def Re(rho, u1, L, mu):
    return rho * u1 * L / mu

u1 = 1
u2 = 0
sigma = 13.5
Cmu = 0.9
x = np.arange(10, 100.5, .5)
y = np.arange(-30, 30.5, .5)
x,y = np.meshgrid(x, y)
uplot = u(x, y, u1, u2, sigma)
vplot = v(x, y, u1, u2, sigma)
kplot = k(x, y, k0(u1, u2, sigma), sigma)
epsPlot = epsilon(x, y, epsilon0(u1, u2, sigma, Cmu), sigma)
muTplot = muT(x, y, muT0(u1, 1), sigma)



In [25]:

    
import matplotlib.pyplot as plt

plt.pcolor(x, y, uplot)
plt.colorbar()
plt.show()



In [27]:

    
plt.pcolor(x, y, vplot)
plt.colorbar()
plt.show()



In [28]:

    
plt.pcolor(x, y, kplot)
plt.colorbar()
plt.show()



In [29]:

    
plt.pcolor(x, y, epsPlot)
plt.colorbar()
plt.show()



In [31]:

    
plt.pcolor(x, y, muTplot)
plt.colorbar()
plt.show()



In [1]:

    
import sympy as sp
from sympy import diff

x, y, sigma, Cmu, rho, mu, k0, eps0 = sp.var('x y sigma Cmu rho mu k0 eps0')



In [2]:

    
def gradVec2(u_vec, x, y):
    return sp.Matrix([[diff(u_vec[0], x), diff(u_vec[1],x)], [diff(u_vec[0], y), diff(u_vec[1], y)]])

def divTen2(tensor, x, y):
    return sp.Matrix([diff(tensor[0,0], x) + diff(tensor[1,0], y), diff(tensor[0, 1], x) + diff(tensor[1,1], y)])

def divVec2(u_vec, x, y):
    return diff(u_vec[0], x) + diff(u_vec[1], y)



In [4]:

    
u = (1 - sp.erf(sigma * y / x)) / 2
v = sp.exp(-(sigma * y / x)**2) / 2 / sigma / sp.sqrt(sp.pi)
k = k0 * sp.exp(-(sigma * y / x)**2)
eps = eps0 / x * sp.exp(-(sigma * y / x)**2)
muT = rho * Cmu * k**2 / eps

u_vec = sp.Matrix([u, v])
grad_u_vec = gradVec2(u_vec, x, y)
visc_term = divTen2((mu + muT) * (grad_u_vec + grad_u_vec.transpose()), x, y)



In [5]:

    
print(sp.simplify(divVec2(u_vec, x, y)))



In [6]:

    
visc_term









    Out[6]:





Matrix([
[                   -2*Cmu*k0**2*rho*sigma**2*y*(-sigma*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x) + sigma*y**2*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))*exp(-sigma**2*y**2/x**2)/(eps0*x) + 4*sigma**3*y**3*(Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) + 2*sigma*y*(2*Cmu*k0**2*rho*sigma**2*y**2*exp(-sigma**2*y**2/x**2)/(eps0*x**2) + Cmu*k0**2*rho*exp(-sigma**2*y**2/x**2)/eps0)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2) - 4*sigma*y*(Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) + (Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*(2*sigma**3*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) - 2*sigma**3*y**3*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) + 2*sigma*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))],
[4*Cmu*k0**2*rho*sigma**3*y**2*exp(-2*sigma**2*y**2/x**2)/(sqrt(pi)*eps0*x**3) + 4*sigma**3*y**2*(Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**4) - 2*sigma*(Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2) + (-sigma*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x) + sigma*y**2*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))*(2*Cmu*k0**2*rho*sigma**2*y**2*exp(-sigma**2*y**2/x**2)/(eps0*x**2) + Cmu*k0**2*rho*exp(-sigma**2*y**2/x**2)/eps0) + (Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*(-2*sigma**3*y**2*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**4) + 2*sigma**3*y**4*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**6) + sigma*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2) - 3*sigma*y**2*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**4))]])



In [61]:

    
visc_term.shape









    Out[61]:





(2, 1)



In [7]:

    
momentum_equations = rho * u_vec.transpose() * grad_u_vec - visc_term.transpose()



In [12]:

    
u_eq = momentum_equations[0]
v_eq = momentum_equations[1]



In [16]:

    
sp.simplify(v_eq)









    Out[16]:





(-8*pi**5*Cmu*k0**2*rho*sigma**3*x**3*y**2 + 2*pi**5*Cmu*k0**2*rho*sigma*x*(x**2 - y**2)*(2*sigma**2*y**2 + x**2) - pi**5*eps0*rho*sigma*x**3*y**2*(erf(sigma*y/x) - 1)*exp(sigma**2*y**2/x**2) - pi**(9/2)*eps0*rho*x**4*y - 8*pi**5*sigma**3*x**2*y**2*(Cmu*k0**2*rho*x + eps0*mu*exp(sigma**2*y**2/x**2)) + 4*pi**5*sigma*x**4*(Cmu*k0**2*rho*x + eps0*mu*exp(sigma**2*y**2/x**2)) + 2*pi**5*sigma*(Cmu*k0**2*rho*x + eps0*mu*exp(sigma**2*y**2/x**2))*(-2*sigma**2*y**4 - x**4 + x**2*y**2*(2*sigma**2 + 3)))*exp(-2*sigma**2*y**2/x**2)/(2*pi**(11/2)*eps0*x**6)



In [13]:

    
sp.simplify(u_eq)









    Out[13]:





(-8*pi**3*Cmu*k0**2*rho*sigma**3*x**3*y - 8*pi**3*Cmu*k0**2*rho*sigma**3*x*y**3 - 4*pi**3*eps0*mu*sigma**3*x**2*y*exp(sigma**2*y**2/x**2) - 4*pi**3*eps0*mu*sigma**3*y**3*exp(sigma**2*y**2/x**2) + 4*pi**3*eps0*mu*sigma*x**2*y*exp(sigma**2*y**2/x**2) - pi**3*eps0*rho*sigma*x**3*y*exp(sigma**2*y**2/x**2)*erf(sigma*y/x) + pi**3*eps0*rho*sigma*x**3*y*exp(sigma**2*y**2/x**2) - pi**(5/2)*eps0*rho*x**4)*exp(-2*sigma**2*y**2/x**2)/(2*pi**(7/2)*eps0*x**5)



In [14]:

    
sp.collect(u_eq, x)









    Out[14]:





-4*sigma**3*y**3*(Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) + 4*sigma*y*(Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) - (Cmu*k0**2*rho*x*exp(-sigma**2*y**2/x**2)/eps0 + mu)*(2*sigma**3*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) - 2*sigma**3*y**3*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) + 2*sigma*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3)) + (2*Cmu*k0**2*rho*sigma**2*y*(-sigma*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x) + sigma*y**2*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))*exp(-sigma**2*y**2/x**2)/eps0 - rho*exp(-2*sigma**2*y**2/x**2)/(2*pi))/x + (rho*sigma*y*(-erf(sigma*y/x)/2 + 1/2)*exp(-sigma**2*y**2/x**2)/sqrt(pi) - 2*sigma*y*(2*Cmu*k0**2*rho*sigma**2*y**2*exp(-sigma**2*y**2/x**2)/(eps0*x**2) + Cmu*k0**2*rho*exp(-sigma**2*y**2/x**2)/eps0)*exp(-sigma**2*y**2/x**2)/sqrt(pi))/x**2



In [10]:

    
u_eq = u_eq.subs(k0, sigma / sp.sqrt(sp.pi) * 343 / 75000)
print(u_eq)
u_eq = u_eq.subs(eps0, Cmu * sigma**2 / sp.pi * 343 / 22500)
print(u_eq)









    



117649*Cmu*rho*sigma**4*y*(-sigma*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x) + sigma*y**2*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))*exp(-sigma**2*y**2/x**2)/(2812500000*pi*eps0*x) + rho*sigma*y*(-erf(sigma*y/x)/2 + 1/2)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2) - rho*exp(-2*sigma**2*y**2/x**2)/(2*pi*x) - 4*sigma**3*y**3*(117649*Cmu*rho*sigma**2*x*exp(-sigma**2*y**2/x**2)/(5625000000*pi*eps0) + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) - 2*sigma*y*(117649*Cmu*rho*sigma**4*y**2*exp(-sigma**2*y**2/x**2)/(2812500000*pi*eps0*x**2) + 117649*Cmu*rho*sigma**2*exp(-sigma**2*y**2/x**2)/(5625000000*pi*eps0))*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2) + 4*sigma*y*(117649*Cmu*rho*sigma**2*x*exp(-sigma**2*y**2/x**2)/(5625000000*pi*eps0) + mu)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) - (117649*Cmu*rho*sigma**2*x*exp(-sigma**2*y**2/x**2)/(5625000000*pi*eps0) + mu)*(2*sigma**3*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) - 2*sigma**3*y**3*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) + 2*sigma*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))
343*rho*sigma**2*y*(-sigma*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x) + sigma*y**2*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))*exp(-sigma**2*y**2/x**2)/(125000*x) + rho*sigma*y*(-erf(sigma*y/x)/2 + 1/2)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2) - rho*exp(-2*sigma**2*y**2/x**2)/(2*pi*x) - 4*sigma**3*y**3*(mu + 343*rho*x*exp(-sigma**2*y**2/x**2)/250000)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) - 2*sigma*y*(343*rho*sigma**2*y**2*exp(-sigma**2*y**2/x**2)/(125000*x**2) + 343*rho*exp(-sigma**2*y**2/x**2)/250000)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2) + 4*sigma*y*(mu + 343*rho*x*exp(-sigma**2*y**2/x**2)/250000)*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) - (mu + 343*rho*x*exp(-sigma**2*y**2/x**2)/250000)*(2*sigma**3*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3) - 2*sigma**3*y**3*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**5) + 2*sigma*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**3))



In [11]:

    
sp.simplify(u_eq)









    Out[11]:





(-125000*pi**3*mu*sigma**3*x**2*y*exp(sigma**2*y**2/x**2) - 125000*pi**3*mu*sigma**3*y**3*exp(sigma**2*y**2/x**2) + 125000*pi**3*mu*sigma*x**2*y*exp(sigma**2*y**2/x**2) - 343*pi**3*rho*sigma**3*x**3*y - 343*pi**3*rho*sigma**3*x*y**3 - 31250*pi**3*rho*sigma*x**3*y*exp(sigma**2*y**2/x**2)*erf(sigma*y/x) + 31250*pi**3*rho*sigma*x**3*y*exp(sigma**2*y**2/x**2) - 31250*pi**(5/2)*rho*x**4)*exp(-2*sigma**2*y**2/x**2)/(62500*pi**(7/2)*x**5)



In [52]:

    
grad_u_vec = sp.Matrix([[diff(u, x), diff(v, x)], [diff(u, y), diff(v, y)]])



In [57]:

    
grad_u_vec









    Out[57]:





Matrix([
[sigma*y*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x**2), -y*exp(sigma*y/x)/(2*sqrt(pi)*x**2)],
[    -sigma*exp(-sigma**2*y**2/x**2)/(sqrt(pi)*x),       exp(sigma*y/x)/(2*sqrt(pi)*x)]])



In [40]:

    
clear(pi)



In [49]:

    
from sympy.physics.vector import ReferenceFrame
R = ReferenceFrame('R')

u = (1 - sp.erf(sigma * R[1] / R[0])) / 2
v = sp.exp(sigma * R[1] / R[0]) / 2 / sigma / sp.sqrt(pi)
k = k0 * sp.exp(-(sigma * R[1] / R[0])**2)
eps = eps0 / R[0] * sp.exp(-(sigma * R[1] / R[0])**2)
muT = rho * Cmu * k**2 / eps



In [50]:

    
u_vec[0]









    Out[50]:





-erf(sigma*y/x)/2 + 1/2



In [ ]:

    
grad_u_vec = gradVec2(u_vec)



In [18]:

    
from scipy.special import erf

erf(2)









    Out[18]:





0.99532226501895271



In [19]:

    
erf(-1)









    Out[19]:





-0.84270079294971478



In [20]:

    
erf(.99)









    Out[20]:





0.83850806955536983



In [47]:

    
from numpy import pi, sqrt, exp

def d_erf(x):
    return 2. / sqrt(pi) * exp(-x**2)

def d_half_erf(x):
    return 2. / sqrt(pi) * exp(-(0.5*x)**2) * 0.5



In [48]:

    
d_half_erf(-2)









    Out[48]:





0.20755374871029736



In [22]:

    
d_erf(-1)









    Out[22]:





0.41510749742059472



In [23]:

    
print(pi)









    



3.141592653589793



In [24]:

    
d_erf(0)









    Out[24]:





1.1283791670955126



In [25]:

    
d_erf(1)









    Out[25]:





0.41510749742059472



In [27]:

    
import numpy as np



In [28]:

    
libmesh = np.loadtxt("/home/lindsayad/projects/moose/libmesh/contrib/fparser/examples/first_orig.dat")



In [30]:

    
libmesh.shape









    Out[30]:





(200, 2)



In [42]:

    
xl = libmesh[:,0]
ypl = libmesh[:,1]



In [45]:

    
xt = np.arange(-1,1,.01)
yt = d_erf(xt)



In [33]:

    
import matplotlib.pyplot as plt



In [43]:

    
plt.close()
plt.plot(xl, ypl, label="libmesh")
plt.plot(xt, yt, label='true')
plt.legend()
plt.show()



In [46]:

    
plt.close()
plt.plot(xl, yt / ypl)
plt.show()

5/10/17

Ok, the analytic_turbulence problem sucks. Even if I start with Dirichlet boundary conditions on all boundaries and initial conditions representing the supposed analytic solution, and then run a transient simulation, the solution evolves away from the supposed analytic solution. backwards_step_adaptive.i runs to completion, but that's for a relatively low inlet velocity.

Getting some pretty good results now also with backwards_step_adaptive_inlet_v_100.i which I wasn't a few days before. This could perhaps be due to the introduction of the SUPG terms. Convergence becomes a little slow at longer time steps, perhaps because of incomplete Jacobian implementation? Or poor relative scaling of the variables? Results for kin actually don't look too far off from the results in the Kuzmin paper. This simulation uses a Reynolds number of 100, which is still pretty small! Next effort will be with the Reynolds number in the Kuzmin paper of 47,625.

It's something I've observed over the years that decreasing element size can lead to decreasing solver convergence. Note that I'm not talking about convergence to the true solution. I wish I could find a good piece of literature discussing this phenomenon. There are just so many things to consider about a finite element solve; it can be fun at times and frustrating at others.

5/11/17

Ok, going to do some methods of manufactured solutions!



In [1]:

    
from sympy import *



In [8]:

    
x, y, L = var('x y L')



In [10]:

    
from random import randint, random, uniform
for i in range(30):
    print('%.2f' % uniform(.1, .99))



In [64]:

    
from random import randint, random, uniform

def sym_func(x, y, L):
    return round(uniform(.1, .99),1) + round(uniform(.1, .99),1) * sin(round(uniform(.1, .99),1) * pi * x / L) \
                            + round(uniform(.1, .99),1) * sin(round(uniform(.1, .99),1) * pi * y / L) \
                            + round(uniform(.1, .99),1) * sin(round(uniform(.1, .99),1) * pi * x * y / L)



In [65]:

    
u = sym_func(x, y, 1)
v = sym_func(x, y, 1)
p = sym_func(x, y, 1)
k = sym_func(x, y, 1)
eps = sym_func(x, y, 1)

print(u, v, p, k, eps, sep="\n")









    



0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5
0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3
0.5*sin(0.5*pi*x) + 1.0*sin(0.3*pi*y) + 0.5*sin(0.2*pi*x*y) + 0.5
0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4
0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5



In [113]:

    
import sympy as sp

def gradVec2(u_vec, x, y):
    return sp.Matrix([[diff(u_vec[0], x), diff(u_vec[1],x)], [diff(u_vec[0], y), diff(u_vec[1], y)]])

def divTen2(tensor, x, y):
    return sp.Matrix([diff(tensor[0,0], x) + diff(tensor[1,0], y), diff(tensor[0, 1], x) + diff(tensor[1,1], y)])

def divVec2(u_vec, x, y):
    return diff(u_vec[0], x) + diff(u_vec[1], y)

def gradScalar2(u, x, y):
    return sp.Matrix([diff(u, x), diff(u,y)])

def strain_rate(u_vec, x, y):
    return gradVec2(u_vec, x, y) + gradVec2(u_vec, x, y).transpose()

def strain_rate_squared_2(u_vec, x, y):
    tensor = gradVec2(u_vec, x, y) + gradVec2(u_vec, x, y).transpose()
    rv = 0
    for i in range(2):
        for j in range(2):
            rv += tensor[i, j] * tensor[i, j]
    return rv

def laplace2(u, x, y):
    return diff(diff(u, x), x) + diff(diff(u, y), y)

Momentum equations



In [94]:

    
pnew = Integer(0)
type(pnew)









    Out[94]:





sympy.core.numbers.Zero

Traction Form



In [104]:

    
cmu = 0.09

uvec = sp.Matrix([u, v])

mu, rho = var('mu rho')

visc_term = (-mu * divTen2(gradVec2(uvec, x, y) + gradVec2(uvec, x, y).transpose(), x, y)).transpose()
conv_term = rho * uvec.transpose() * gradVec2(uvec, x, y)
pressure_term = gradScalar2(p, x, y).transpose()
turbulent_visc_term = -(divTen2(rho * cmu * k**2 / eps * (gradVec2(uvec, x, y) + gradVec2(uvec, x, y).transpose()), x, y)).transpose()
# print(visc_term.shape, conv_term.shape, pressure_term.shape, sep="\n")
source = conv_term + visc_term + pressure_term + turbulent_visc_term
print(source[0])
print(source[1])









    



-mu*(-0.028*pi**2*x**2*sin(0.2*pi*x*y) - 0.018*pi**2*x*y*sin(0.3*pi*x*y) - 0.056*pi**2*y**2*sin(0.2*pi*x*y) - 0.2*pi**2*sin(0.5*pi*x) - 0.4*pi**2*sin(pi*y) + 0.06*pi*cos(0.3*pi*x*y)) + rho*(0.14*pi*x*cos(0.2*pi*x*y) + 0.4*pi*cos(pi*y))*(0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3) - 0.09*rho*(0.56*pi*x*cos(0.4*pi*x*y) + 1.26*pi*cos(0.7*pi*y))*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(-0.48*pi*x*cos(0.6*pi*x*y) - 0.81*pi*cos(0.9*pi*y))*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 + rho*(0.14*pi*y*cos(0.2*pi*x*y) + 0.2*pi*cos(0.5*pi*x))*(0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5) - 0.09*rho*(0.28*pi*y*cos(0.2*pi*x*y) + 0.4*pi*cos(0.5*pi*x))*(0.56*pi*y*cos(0.4*pi*x*y) + 0.56*pi*cos(0.7*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(0.28*pi*y*cos(0.2*pi*x*y) + 0.4*pi*cos(0.5*pi*x))*(-0.48*pi*y*cos(0.6*pi*x*y) - 0.18*pi*cos(0.3*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 - 0.09*rho*(-0.056*pi**2*y**2*sin(0.2*pi*x*y) - 0.2*pi**2*sin(0.5*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(-0.028*pi**2*x**2*sin(0.2*pi*x*y) - 0.018*pi**2*x*y*sin(0.3*pi*x*y) - 0.4*pi**2*sin(pi*y) + 0.06*pi*cos(0.3*pi*x*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.1*pi*y*cos(0.2*pi*x*y) + 0.25*pi*cos(0.5*pi*x)
-mu*(-0.036*pi**2*x**2*sin(0.3*pi*x*y) - 0.028*pi**2*x*y*sin(0.2*pi*x*y) - 0.018*pi**2*y**2*sin(0.3*pi*x*y) - 0.384*pi**2*sin(0.8*pi*x) - 0.054*pi**2*sin(0.3*pi*y) + 0.14*pi*cos(0.2*pi*x*y)) + rho*(0.06*pi*x*cos(0.3*pi*x*y) + 0.09*pi*cos(0.3*pi*y))*(0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3) - 0.09*rho*(0.12*pi*x*cos(0.3*pi*x*y) + 0.18*pi*cos(0.3*pi*y))*(0.56*pi*x*cos(0.4*pi*x*y) + 1.26*pi*cos(0.7*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(0.12*pi*x*cos(0.3*pi*x*y) + 0.18*pi*cos(0.3*pi*y))*(-0.48*pi*x*cos(0.6*pi*x*y) - 0.81*pi*cos(0.9*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 + rho*(0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x))*(0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5) - 0.09*rho*(0.56*pi*y*cos(0.4*pi*x*y) + 0.56*pi*cos(0.7*pi*x))*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(-0.48*pi*y*cos(0.6*pi*x*y) - 0.18*pi*cos(0.3*pi*x))*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 - 0.09*rho*(-0.036*pi**2*x**2*sin(0.3*pi*x*y) - 0.054*pi**2*sin(0.3*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(-0.028*pi**2*x*y*sin(0.2*pi*x*y) - 0.018*pi**2*y**2*sin(0.3*pi*x*y) - 0.384*pi**2*sin(0.8*pi*x) + 0.14*pi*cos(0.2*pi*x*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.1*pi*x*cos(0.2*pi*x*y) + 0.3*pi*cos(0.3*pi*y)

Laplace Form



In [105]:

    
cmu = 0.09

uvec = sp.Matrix([u, v])

mu, rho = var('mu rho')

visc_term = (-mu * divTen2(gradVec2(uvec, x, y), x, y)).transpose()
conv_term = rho * uvec.transpose() * gradVec2(uvec, x, y)
pressure_term = gradScalar2(p, x, y).transpose()
turbulent_visc_term = -(divTen2(rho * cmu * k**2 / eps * (gradVec2(uvec, x, y)), x, y)).transpose()
# print(visc_term.shape, conv_term.shape, pressure_term.shape, sep="\n")
source = conv_term + visc_term + pressure_term + turbulent_visc_term
print(source[0])
print(source[1])









    



-mu*(-0.028*pi**2*x**2*sin(0.2*pi*x*y) - 0.028*pi**2*y**2*sin(0.2*pi*x*y) - 0.1*pi**2*sin(0.5*pi*x) - 0.4*pi**2*sin(pi*y)) - 0.09*rho*(0.14*pi*x*cos(0.2*pi*x*y) + 0.4*pi*cos(pi*y))*(0.56*pi*x*cos(0.4*pi*x*y) + 1.26*pi*cos(0.7*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(0.14*pi*x*cos(0.2*pi*x*y) + 0.4*pi*cos(pi*y))*(-0.48*pi*x*cos(0.6*pi*x*y) - 0.81*pi*cos(0.9*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 + rho*(0.14*pi*x*cos(0.2*pi*x*y) + 0.4*pi*cos(pi*y))*(0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3) - 0.09*rho*(0.14*pi*y*cos(0.2*pi*x*y) + 0.2*pi*cos(0.5*pi*x))*(0.56*pi*y*cos(0.4*pi*x*y) + 0.56*pi*cos(0.7*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(0.14*pi*y*cos(0.2*pi*x*y) + 0.2*pi*cos(0.5*pi*x))*(-0.48*pi*y*cos(0.6*pi*x*y) - 0.18*pi*cos(0.3*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 + rho*(0.14*pi*y*cos(0.2*pi*x*y) + 0.2*pi*cos(0.5*pi*x))*(0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5) - 0.09*rho*(-0.028*pi**2*x**2*sin(0.2*pi*x*y) - 0.4*pi**2*sin(pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(-0.028*pi**2*y**2*sin(0.2*pi*x*y) - 0.1*pi**2*sin(0.5*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.1*pi*y*cos(0.2*pi*x*y) + 0.25*pi*cos(0.5*pi*x)
-mu*(-0.018*pi**2*x**2*sin(0.3*pi*x*y) - 0.018*pi**2*y**2*sin(0.3*pi*x*y) - 0.384*pi**2*sin(0.8*pi*x) - 0.027*pi**2*sin(0.3*pi*y)) - 0.09*rho*(0.06*pi*x*cos(0.3*pi*x*y) + 0.09*pi*cos(0.3*pi*y))*(0.56*pi*x*cos(0.4*pi*x*y) + 1.26*pi*cos(0.7*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(0.06*pi*x*cos(0.3*pi*x*y) + 0.09*pi*cos(0.3*pi*y))*(-0.48*pi*x*cos(0.6*pi*x*y) - 0.81*pi*cos(0.9*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 + rho*(0.06*pi*x*cos(0.3*pi*x*y) + 0.09*pi*cos(0.3*pi*y))*(0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3) - 0.09*rho*(0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x))*(0.56*pi*y*cos(0.4*pi*x*y) + 0.56*pi*cos(0.7*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x))*(-0.48*pi*y*cos(0.6*pi*x*y) - 0.18*pi*cos(0.3*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2 + rho*(0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x))*(0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5) - 0.09*rho*(-0.018*pi**2*x**2*sin(0.3*pi*x*y) - 0.027*pi**2*sin(0.3*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.09*rho*(-0.018*pi**2*y**2*sin(0.3*pi*x*y) - 0.384*pi**2*sin(0.8*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.1*pi*x*cos(0.2*pi*x*y) + 0.3*pi*cos(0.3*pi*y)

Pressure equation



In [73]:

    
-divVec2(uvec, x, y)









    Out[73]:





-0.06*pi*x*cos(0.3*pi*x*y) - 0.14*pi*y*cos(0.2*pi*x*y) - 0.2*pi*cos(0.5*pi*x) - 0.09*pi*cos(0.3*pi*y)

Or testing with a simple diffusion term



In [96]:

    
diff_term = -laplace2(p, x, y)
print(diff_term)









    



0.02*pi**2*x**2*sin(0.2*pi*x*y) + 0.02*pi**2*y**2*sin(0.2*pi*x*y) + 0.125*pi**2*sin(0.5*pi*x) + 0.09*pi**2*sin(0.3*pi*y)

Turbulent kinetic energy equation



In [99]:

    
cmu = 0.09
sigk = 1.
sigeps = 1.3
c1eps = 1.44
c2eps = 1.92

conv_term = rho * uvec.transpose() * gradScalar2(k, x, y)
diff_term = - divVec2((mu + rho * cmu * k**2 / eps / sigk) * gradScalar2(k, x, y), x, y)
creation_term = - rho * cmu * k**2 / eps / 2 * strain_rate_squared_2(uvec, x, y)
destruction_term = rho * eps

terms = [conv_term[0,0], diff_term, creation_term, destruction_term]
L = 0
for term in terms:
    L += term
print(L)









    



rho*(0.28*pi*x*cos(0.4*pi*x*y) + 0.63*pi*cos(0.7*pi*y))*(0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3) + rho*(0.28*pi*y*cos(0.4*pi*x*y) + 0.28*pi*cos(0.7*pi*x))*(0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5) - 0.045*rho*((0.12*pi*x*cos(0.3*pi*x*y) + 0.18*pi*cos(0.3*pi*y))**2 + (0.28*pi*y*cos(0.2*pi*x*y) + 0.4*pi*cos(0.5*pi*x))**2 + 2*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))**2)*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + rho*(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - (mu + 0.09*rho*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5))*(-0.112*pi**2*x**2*sin(0.4*pi*x*y) - 0.441*pi**2*sin(0.7*pi*y)) - (mu + 0.09*rho*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5))*(-0.112*pi**2*y**2*sin(0.4*pi*x*y) - 0.196*pi**2*sin(0.7*pi*x)) - (0.28*pi*x*cos(0.4*pi*x*y) + 0.63*pi*cos(0.7*pi*y))*(0.09*rho*(0.56*pi*x*cos(0.4*pi*x*y) + 1.26*pi*cos(0.7*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.09*rho*(-0.48*pi*x*cos(0.6*pi*x*y) - 0.81*pi*cos(0.9*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2) - (0.28*pi*y*cos(0.4*pi*x*y) + 0.28*pi*cos(0.7*pi*x))*(0.09*rho*(0.56*pi*y*cos(0.4*pi*x*y) + 0.56*pi*cos(0.7*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.09*rho*(-0.48*pi*y*cos(0.6*pi*x*y) - 0.18*pi*cos(0.3*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2)

Turbulent dissipation



In [100]:

    
cmu = 0.09
sigk = 1.
sigeps = 1.3
c1eps = 1.44
c2eps = 1.92

conv_term = rho * uvec.transpose() * gradScalar2(eps, x, y)
diff_term = - divVec2((mu + rho * cmu * k**2 / eps / sigeps) * gradScalar2(eps, x, y), x, y)
creation_term = - rho * c1eps * cmu * k / 2 * strain_rate_squared_2(uvec, x, y)
destruction_term = rho * c2eps * eps**2 / k

terms = [conv_term[0,0], diff_term, creation_term, destruction_term]
L = 0
for term in terms:
    L += term
print(L)









    



rho*(0.48*pi*x*cos(0.6*pi*x*y) + 0.81*pi*cos(0.9*pi*y))*(0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3) + rho*(0.48*pi*y*cos(0.6*pi*x*y) + 0.18*pi*cos(0.3*pi*x))*(0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5) - 0.0648*rho*((0.12*pi*x*cos(0.3*pi*x*y) + 0.18*pi*cos(0.3*pi*y))**2 + (0.28*pi*y*cos(0.2*pi*x*y) + 0.4*pi*cos(0.5*pi*x))**2 + 2*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))**2)*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4) + 1.92*rho*(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2/(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4) - (mu + 0.0692307692307692*rho*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5))*(-0.288*pi**2*x**2*sin(0.6*pi*x*y) - 0.729*pi**2*sin(0.9*pi*y)) - (mu + 0.0692307692307692*rho*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5))*(-0.288*pi**2*y**2*sin(0.6*pi*x*y) - 0.054*pi**2*sin(0.3*pi*x)) - (0.48*pi*x*cos(0.6*pi*x*y) + 0.81*pi*cos(0.9*pi*y))*(0.0692307692307692*rho*(0.56*pi*x*cos(0.4*pi*x*y) + 1.26*pi*cos(0.7*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.0692307692307692*rho*(-0.48*pi*x*cos(0.6*pi*x*y) - 0.81*pi*cos(0.9*pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2) - (0.48*pi*y*cos(0.6*pi*x*y) + 0.18*pi*cos(0.3*pi*x))*(0.0692307692307692*rho*(0.56*pi*y*cos(0.4*pi*x*y) + 0.56*pi*cos(0.7*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) + 0.0692307692307692*rho*(-0.48*pi*y*cos(0.6*pi*x*y) - 0.18*pi*cos(0.3*pi*x))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)**2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)**2)

Simple diffusion



In [93]:

    
diff_term = -laplace2(u, x, y)
print(diff_term)









    



0.028*pi**2*x**2*sin(0.2*pi*x*y) + 0.028*pi**2*y**2*sin(0.2*pi*x*y) + 0.1*pi**2*sin(0.5*pi*x) + 0.4*pi**2*sin(pi*y)



In [69]:

    
def z(func, xh, yh):
    u = np.zeros(xh.shape)
    for i in range(0,xh.shape[0]):
        for j in range(0,xh.shape[1]):
            u[i][j] = func.subs({x:xh[i][j], y:yh[i][j]}).evalf()
#             print(func.subs({x:xh[i][j], y:yh[i][j]}).evalf())
    return u
            

xnum = np.arange(0, 1.01, .05)
ynum = np.arange(0, 1.01, .05)
xgrid, ygrid = np.meshgrid(xnum, ynum)
uh = z(u, xgrid, ygrid)

vh = z(v, xgrid, ygrid)
ph = z(p, xgrid, ygrid)
kh = z(k, xgrid, ygrid)
epsh = z(eps, xgrid, ygrid)



In [72]:

    
import matplotlib.pyplot as plt

plot_funcs = [uh, vh, ph, kh, epsh]

for func in plot_funcs:
    plt.pcolor(xgrid, ygrid, func, cmap='coolwarm')
    cbar = plt.colorbar()
    plt.show()



In [112]:

    
f, g = symbols('f g', cls=Function)
f(x,y).diff(x)









    Out[112]:





$$\frac{\partial}{\partial x} f{\left (x,y \right )}$$



In [123]:

    
vx, vy = symbols('v_x v_y', cls=Function)
mu, x, y = var('mu x y')
nx, ny = var('n_x n_y')
n = sp.Matrix([nx, ny])
v_vec = sp.Matrix([vx(x, y), vy(x, y)])
sigma = strain_rate(v_vec, x, y)

tw = n.transpose() * sigma - n.transpose() * sigma * n * n.transpose()
tw[0]









    Out[123]:





$$- n_{x} \left(n_{x} \left(2 n_{x} \frac{\partial}{\partial x} \operatorname{v_{x}}{\left (x,y \right )} + n_{y} \left(\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{y}}{\left (x,y \right )}\right)\right) + n_{y} \left(n_{x} \left(\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{y}}{\left (x,y \right )}\right) + 2 n_{y} \frac{\partial}{\partial y} \operatorname{v_{y}}{\left (x,y \right )}\right)\right) + 2 n_{x} \frac{\partial}{\partial x} \operatorname{v_{x}}{\left (x,y \right )} + n_{y} \left(\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{y}}{\left (x,y \right )}\right)$$



In [124]:

    
tw[1]









    Out[124]:





$$n_{x} \left(\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{y}}{\left (x,y \right )}\right) - n_{y} \left(n_{x} \left(2 n_{x} \frac{\partial}{\partial x} \operatorname{v_{x}}{\left (x,y \right )} + n_{y} \left(\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{y}}{\left (x,y \right )}\right)\right) + n_{y} \left(n_{x} \left(\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} + \frac{\partial}{\partial x} \operatorname{v_{y}}{\left (x,y \right )}\right) + 2 n_{y} \frac{\partial}{\partial y} \operatorname{v_{y}}{\left (x,y \right )}\right)\right) + 2 n_{y} \frac{\partial}{\partial y} \operatorname{v_{y}}{\left (x,y \right )}$$



In [129]:

    
vx, vy = symbols('v_x v_y', cls=Function)
mu, x, y = var('mu x y')
nx, ny = var('n_x n_y')
n = sp.Matrix([Integer(0), Integer(1)])
v_vec = sp.Matrix([vx(x, y), 0])
sigma = strain_rate(v_vec, x, y)

tw = n.transpose() * sigma - n.transpose() * sigma * n * n.transpose()
tw[0]









    Out[129]:





$$\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )}$$

Verified RANS kernels

INSK
INSEpsilon
INSMomentumTurbulentViscosityTractionForm
INSMomentumTurbulentViscosityLaplaceForm
INSMomentumShearStressWallFunction with |u|/yStarPlus branch of uTau

with exp_form = false.



In [130]:

    
u









    Out[130]:





$$0.4 \sin{\left (0.5 \pi x \right )} + 0.4 \sin{\left (\pi y \right )} + 0.7 \sin{\left (0.2 \pi x y \right )} + 0.5$$



In [113]:

    
import sympy as sp

def gradVec2(u_vec, x, y):
    return sp.Matrix([[diff(u_vec[0], x), diff(u_vec[1],x)], [diff(u_vec[0], y), diff(u_vec[1], y)]])

def divTen2(tensor, x, y):
    return sp.Matrix([diff(tensor[0,0], x) + diff(tensor[1,0], y), diff(tensor[0, 1], x) + diff(tensor[1,1], y)])

def divVec2(u_vec, x, y):
    return diff(u_vec[0], x) + diff(u_vec[1], y)

def gradScalar2(u, x, y):
    return sp.Matrix([diff(u, x), diff(u,y)])

def strain_rate(u_vec, x, y):
    return gradVec2(u_vec, x, y) + gradVec2(u_vec, x, y).transpose()

def strain_rate_squared_2(u_vec, x, y):
    tensor = gradVec2(u_vec, x, y) + gradVec2(u_vec, x, y).transpose()
    rv = 0
    for i in range(2):
        for j in range(2):
            rv += tensor[i, j] * tensor[i, j]
    return rv

def laplace2(u, x, y):
    return diff(diff(u, x), x) + diff(diff(u, y), y)



In [210]:

    
def L_momentum_traction(uvec, k, eps, x, y):
    cmu = 0.09
    mu, rho = sp.var('mu rho')
    visc_term = (-mu * divTen2(gradVec2(uvec, x, y) + gradVec2(uvec, x, y).transpose(), x, y)).transpose()
    conv_term = rho * uvec.transpose() * gradVec2(uvec, x, y)
    pressure_term = gradScalar2(p, x, y).transpose()
    turbulent_visc_term = -(divTen2(rho * cmu * k**2 / eps * (gradVec2(uvec, x, y) + gradVec2(uvec, x, y).transpose()), x, y)).transpose()
    # print(visc_term.shape, conv_term.shape, pressure_term.shape, sep="\n")
    source = conv_term + visc_term + pressure_term + turbulent_visc_term
    return source

def bc_terms_momentum_traction(uvec, nvec, k, eps, x, y):
    cmu = 0.09
    mu, rho = sp.var('mu rho')
    visc_term = (-mu * nvec.transpose() * (gradVec2(uvec, x, y) + gradVec2(uvec, x, y).transpose(), x, y)).transpose()
    turbulent_visc_term = -(nvec.transpose() * (rho * cmu * k**2 / eps * (gradVec2(uvec, x, y) + gradVec2(uvec, x, y).transpose()), x, y)).transpose()
    return visc_term + turbulent_visc_term
    
def L_momentum_laplace(uvec, k, eps, x, y):
    cmu = 0.09
    mu, rho = var('mu rho')
    visc_term = (-mu * divTen2(gradVec2(uvec, x, y), x, y)).transpose()
    conv_term = rho * uvec.transpose() * gradVec2(uvec, x, y)
    pressure_term = gradScalar2(p, x, y).transpose()
    turbulent_visc_term = -(divTen2(rho * cmu * k**2 / eps * (gradVec2(uvec, x, y)), x, y)).transpose()
    # print(visc_term.shape, conv_term.shape, pressure_term.shape, sep="\n")
    source = conv_term + visc_term + pressure_term + turbulent_visc_term
    return source

def L_pressure(uvec, x, y):
    return -divVec2(uvec, x, y)

def L_kin(uvec, k, eps, x, y):
    cmu = 0.09
    sigk = 1.
    sigeps = 1.3
    c1eps = 1.44
    c2eps = 1.92
    conv_term = rho * uvec.transpose() * gradScalar2(k, x, y)
    diff_term = - divVec2((mu + rho * cmu * k**2 / eps / sigk) * gradScalar2(k, x, y), x, y)
    creation_term = - rho * cmu * k**2 / eps / 2 * strain_rate_squared_2(uvec, x, y)
    destruction_term = rho * eps
    terms = [conv_term[0,0], diff_term, creation_term, destruction_term]
    L = 0
    for term in terms:
        L += term
    return L

def L_eps(uvec, k, eps, x, y):
    cmu = 0.09
    sigk = 1.
    sigeps = 1.3
    c1eps = 1.44
    c2eps = 1.92
    conv_term = rho * uvec.transpose() * gradScalar2(eps, x, y)
    diff_term = - divVec2((mu + rho * cmu * k**2 / eps / sigeps) * gradScalar2(eps, x, y), x, y)
    creation_term = - rho * c1eps * cmu * k / 2 * strain_rate_squared_2(uvec, x, y)
    destruction_term = rho * c2eps * eps**2 / k
    terms = [conv_term[0,0], diff_term, creation_term, destruction_term]
    L = 0
    for term in terms:
        L += term
    return L



In [159]:

    
def prep_moose_input(sym_expr):
    rep1 = re.sub(r'\*\*',r'^',str(sym_expr))
    rep2 = re.sub(r'mu',r'${mu}',rep1)
    rep3 = re.sub(r'rho',r'${rho}',rep2)
    return rep3

def write_all_functions():
    target = open('/home/lindsayad/python/mms_input.txt','w')
    target.write("[Functions]" + "\n")
    target.write("  [./u_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_momentum_traction(uVecNew, kinNew, epsilonNew, x, y)[0]) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./v_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_momentum_traction(uVecNew, kinNew, epsilonNew, x, y)[1]) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./p_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_pressure(uVecNew, x, y)) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./kin_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_kin(uVecNew, kinNew, epsilonNew, x, y)) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./epsilon_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_eps(uVecNew, kinNew, epsilonNew, x, y)) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./u_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(uNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./v_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(vNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./p_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(pNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./kin_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(kinNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./epsilon_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(epsilonNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("[]" + "\n")
    target.close()
    
def write_reduced_functions():
    target = open('/home/lindsayad/python/mms_input.txt','w')
    target.write("[Functions]" + "\n")
    target.write("  [./u_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_momentum_traction(uVecNew, kinNew, epsilonNew, x, y)[0]) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./kin_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_kin(uVecNew, kinNew, epsilonNew, x, y)) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./epsilon_source_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + prep_moose_input(L_eps(uVecNew, kinNew, epsilonNew, x, y)) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./u_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(uNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./kin_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(kinNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("  [./epsilon_func]" + "\n")
    target.write("    type = ParsedFunction" + "\n")
    target.write("    value = '" + str(epsilonNew) + "'" + "\n")
    target.write("  [../]" + "\n")
    target.write("[]" + "\n")
    target.close()



In [205]:

    
yStarPlus = 11.06

# uNew = yStarPlus**2 / y + u * (y - 1.) * 200
# uNew = u * (y - 1.) * 200
# vNew = Integer(0)
# vNew = v * (y - 1.5) * 200

# pNew = Integer(0)

# # Converges
# uNew = u
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Does not converge
# uNew = u * 200
# vNew = v * 200
# pNew = p * 200
# kinNew = k * 200
# epsilonNew = eps * 200

# # Converges
# uNew = u * (y - 1.)
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Converges
# uNew = u * (y - 1.) + 1. / y
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Does not converge
# uNew = u * (y - 1.) + yStarPlus**2 / y
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Converges
# uNew = u * (y - 1.) + 1.1 / y
# vNew = 0
# pNew = 0
# kinNew = k
# epsilonNew = eps

# Want to test natural boundary condition
uNew = 0.5 + sin(pi * x / 2) + sin(pi * y / 2)
vNew = 0
pNew = 0
kinNew = k
epsilonNew = eps

uVecNew = sp.Matrix([uNew, vNew])

write_reduced_functions()



In [169]:

    
print(u)
print(v)
print(p)
print(k)
print(eps)









    



0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5
0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3
0.5*sin(0.5*pi*x) + 1.0*sin(0.3*pi*y) + 0.5*sin(0.2*pi*x*y) + 0.5
0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4
0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5

Ok, so apparently just scaling every variable's manufactured solution by 200 causes MOOSE convergence issues. Sigh



In [ ]:

    
vx, vy = symbols('v_x v_y', cls=Function)
mu, x, y = var('mu x y')
nx, ny = var('n_x n_y')
n = sp.Matrix([Integer(0), Integer(1)])
v_vec = sp.Matrix([vx(x, y), 0])
kinFunc, epsFunc = symbols('kinFunc epsFunc', cls=Function)



In [225]:

    
blah = bc_terms_momentum_traction(uVecNew, n, kinNew, epsilonNew, x, y)



In [218]:

    
type(blah)









    Out[218]:





sympy.matrices.dense.MutableDenseMatrix



In [226]:

    
blah[0]









    Out[226]:





$$\left ( \right )$$



In [207]:

    
sigma = mu * strain_rate(v_vec, x, y)

tw = n.transpose() * sigma - n.transpose() * sigma * n * n.transpose()
tw[0]









    Out[207]:





$$\mu \frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )}$$



In [204]:

    
(n.transpose() * sigma)[0]









    Out[204]:





$$\mu \frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )}$$



In [221]:

    
cmu = 0.09
mu, rho = sp.var('mu rho')
visc_term = (-mu * n.transpose() * (gradVec2(v_vec, x, y) + gradVec2(v_vec, x, y).transpose(), x, y)).transpose()
turbulent_visc_term = -(n.transpose() * (rho * cmu * k**2 / eps * (gradVec2(v_vec, x, y) + gradVec2(v_vec, x, y).transpose()), x, y)).transpose()



In [222]:

    
visc_term









    Out[222]:





$$\left[\begin{matrix}\left ( \right )\\- \mu \left ( \left[\begin{matrix}2 \frac{\partial}{\partial x} \operatorname{v_{x}}{\left (x,y \right )} & \frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )}\\\frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} & 0\end{matrix}\right], \quad x, \quad y\right )\end{matrix}\right]$$



In [223]:

    
print(visc_term)









    



Matrix([[()], [-mu*(Matrix([
[2*Derivative(v_x(x, y), x), Derivative(v_x(x, y), y)],
[  Derivative(v_x(x, y), y),                        0]]), x, y)]])



In [ ]:

    
yStarPlus = 11.06

# uNew = yStarPlus**2 / y + u * (y - 1.) * 200
# uNew = u * (y - 1.) * 200
# vNew = Integer(0)
# vNew = v * (y - 1.5) * 200

# pNew = Integer(0)

# # Converges
# uNew = u
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Does not converge
# uNew = u * 200
# vNew = v * 200
# pNew = p * 200
# kinNew = k * 200
# epsilonNew = eps * 200

# # Converges
# uNew = u * (y - 1.)
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Converges
# uNew = u * (y - 1.) + 1. / y
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Does not converge
# uNew = u * (y - 1.) + yStarPlus**2 / y
# vNew = v
# pNew = p
# kinNew = k
# epsilonNew = eps

# # Converges
# uNew = u * (y - 1.) + 1.1 / y
# vNew = 0
# pNew = 0
# kinNew = k
# epsilonNew = eps



In [1]:

    
from moose_calc_routines import *
from sympy import *

init_printing()

yStarPlus = 1.1
vx, vy = symbols('v_x v_y', cls=Function, positive=True, real=True)
mu, x, y = var('mu x y', real=True, positive=True)
nx, ny = var('n_x n_y')
n = sp.Matrix([Integer(0), Integer(1)])
v_vec = sp.Matrix([vx(x, y), 0])
kinFunc, epsFunc = symbols('k_f \epsilon_f', cls=Function)

u = sym_func(x, y, 1)
v = sym_func(x, y, 1)
p = sym_func(x, y, 1)
k = sym_func(x, y, 1)
eps = sym_func(x, y, 1)

# Want to test wall function bc
uNew = mu * yStarPlus**2 / y
# uNew = 0.5 + sin(pi * x / 2) + sin(pi * y / 2) + mu * yStarPlus**2 / y
vNew = 0
pNew = 0
kinNew = k
epsilonNew = eps

# # Want to test natural boundary condition
# uNew = 0.5 + sin(pi * x / 2) + sin(pi * y / 2)
# vNew = 0
# pNew = 0
# kinNew = k
# epsilonNew = eps

uVecNew = sp.Matrix([uNew, vNew])

numeric = bc_terms_momentum_traction(uVecNew, n, kinNew, epsilonNew, x, y, symbolic=False)
numeric_wall_function = wall_function_momentum_traction(uVecNew, n, kinNew, epsilonNew, x, y, "kin", symbolic=False)
symbolic = bc_terms_momentum_traction(v_vec, n, kinFunc(x, y), epsFunc(x, y), x, y, symbolic=True)
wall_function = wall_function_momentum_traction(v_vec, n, kinFunc(x, y), epsFunc(x, y), x, y, "kin", symbolic=True)



In [4]:

    
write_reduced_functions(uVecNew, kinNew, epsilonNew, x, y)



In [6]:

    
expr = numeric[0] - numeric_wall_function[0]

expr.subs(y, 1).collect('mu')









    Out[6]:





$$- 1.21 \mu^{2} + 1.21 \mu \mu$$



In [2]:

    
expr = symbolic[0] - wall_function[0]
# print(expr)
# expr

newexp = expr.subs(vx(x, y), vx(y)).subs(Abs(vx(y)), vx(y))
# print(newexp)
newexp









    Out[2]:





$$- \frac{c_{\mu}^{0.25}}{y_{\mu}} \sqrt{\operatorname{k_{f}}{\left (x,y \right )}} \operatorname{v_{x}}{\left (y \right )} - \mu \frac{d}{d y} \operatorname{v_{x}}{\left (y \right )}$$



In [45]:

    
dsolve(newexp, vx(y))









    Out[45]:





$$\operatorname{v_{x}}{\left (y \right )} = \frac{\mu y_{\mu}^{2}}{C_{1} y_{\mu}^{2} + y}$$



In [39]:

    
symbolic[0]









    Out[39]:





$$- \frac{c_{\mu} \rho \operatorname{k_{f}}^{2}{\left (x,y \right )}}{\epsilon_{f}{\left (x,y \right )}} \frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} - \mu \frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )}$$



In [40]:

    
wall_function[0]









    Out[40]:





$$- \frac{c_{\mu} \rho \operatorname{k_{f}}^{2}{\left (x,y \right )}}{\epsilon_{f}{\left (x,y \right )}} \frac{\partial}{\partial y} \operatorname{v_{x}}{\left (x,y \right )} + \frac{1}{y_{\mu}^{2}} \operatorname{v_{x}}{\left (x,y \right )} \left|{\operatorname{v_{x}}{\left (x,y \right )}}\right|$$



In [68]:

    
from moose_calc_routines import *
from sympy import *
import sympy as sp
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
init_printing()

x, y = var('x y')

# # INS turbulence
# u = 0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5
# v = 0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3
# p = 0.5*sin(0.5*pi*x) + 1.0*sin(0.3*pi*y) + 0.5*sin(0.2*pi*x*y) + 0.5
# k = 0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4
# eps = 0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5
# uvec = sp.Matrix([u, v])
# n = sp.Matrix([Integer(0), Integer(1)])

# INS only
u = 0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5
v = 0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3
p = 0.5*sin(0.5*pi*x) + 1.0*sin(0.3*pi*y) + 0.5*sin(0.2*pi*x*y) + 0.5
uvec = sp.Matrix([u, v])
nvec = sp.Matrix([Integer(0), Integer(1)])
nvecs = {'left' : sp.Matrix([-1, 0]), 'top' : sp.Matrix([0, 1]), \
         'right' : sp.Matrix([1, 0]), 'bottom' : sp.Matrix([0, -1])}

source = {bnd_name : 
          prep_moose_input(-bc_terms_momentum_traction_no_turbulence(uvec, nvec, p, x, y, parts=True)[0])
          for bnd_name, nvec in nvecs.items()}



In [75]:

    
source









    Out[75]:





{'bottom': '-${mu}*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))',
 'left': '-${mu}*(0.28*pi*y*cos(0.2*pi*x*y) + 0.4*pi*cos(0.5*pi*x)) + 0.5*sin(0.5*pi*x) + 1.0*sin(0.3*pi*y) + 0.5*sin(0.2*pi*x*y) + 0.5',
 'right': '${mu}*(0.28*pi*y*cos(0.2*pi*x*y) + 0.4*pi*cos(0.5*pi*x)) - 0.5*sin(0.5*pi*x) - 1.0*sin(0.3*pi*y) - 0.5*sin(0.2*pi*x*y) - 0.5',
 'top': '${mu}*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))'}



In [3]:

    
surface_terms = bc_terms_momentum_traction_no_turbulence(uvec, nvec, p, x, y, parts=True)
tested_bc = no_bc_bc(uvec, nvec, p, x, y, parts=True)
needed_func = tested_bc - surface_terms
print(prep_moose_input(needed_func[0]))
print(prep_moose_input(needed_func[1]))









    



> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(172)no_bc_bc()
-> if parts:
(Pdb) c
0
0



In [1]:

    
surface_terms = bc_terms_momentum_traction(uvec, n, p, k, eps, x, y, symbolic=False, parts=True)
tested_bc = wall_function_momentum_traction(uvec, n, p, k, eps, x, y, "kin", symbolic=False, parts=True)









    



> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(127)wall_function_momentum_traction()
-> if symbolic:
(Pdb) n
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(131)wall_function_momentum_traction()
-> cmu = 0.09
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(132)wall_function_momentum_traction()
-> yStarPlus = 1.1
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(133)wall_function_momentum_traction()
-> if tau_type == "vel":
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(136)wall_function_momentum_traction()
-> elif tau_type == "kin":
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(137)wall_function_momentum_traction()
-> uTau = cmu**.25 * sp.sqrt(k)
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(141)wall_function_momentum_traction()
-> mu, rho = sp.var('mu rho')
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(142)wall_function_momentum_traction()
-> normal_stress_term = (-nvec.transpose() * mu * strain_rate(uvec, x, y) * nvec * nvec.transpose()).transpose()
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(143)wall_function_momentum_traction()
-> tangential_stress_term = uTau / yStarPlus * uvec
(Pdb) p normal_stress_term
Matrix([
[                                                      0],
[-mu*(0.12*pi*x*cos(0.3*pi*x*y) + 0.18*pi*cos(0.3*pi*y))]])
(Pdb) n
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(144)wall_function_momentum_traction()
-> muT = rho * cmu * k * k / eps
(Pdb) 
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(146)wall_function_momentum_traction()
-> turbulent_stress_term = (-nvec.transpose() * strain_rate(uvec, x, y)).transpose()
(Pdb) s
--Call--
> /home/lindsayad/miniconda3/lib/python3.6/site-packages/sympy/matrices/matrices.py(355)transpose()
-> def transpose(self):
(Pdb) return
--Return--
> /home/lindsayad/miniconda3/lib/python3.6/site-packages/sympy/matrices/matrices.py(356)transpose()->Matrix([[0, 1]])
-> return self._eval_transpose()
(Pdb) s
--Call--
> /home/lindsayad/miniconda3/lib/python3.6/site-packages/sympy/matrices/matrices.py(595)__neg__()
-> def __neg__(self):
(Pdb) return
--Return--
> /home/lindsayad/miniconda3/lib/python3.6/site-packages/sympy/matrices/matrices.py(596)__neg__()->Matrix([[0, -1]])
-> return -1*self
(Pdb) s
--Call--
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(26)strain_rate()
-> def strain_rate(u_vec, x, y):
(Pdb) n
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(27)strain_rate()
-> return gradVec2(u_vec, x, y) + gradVec2(u_vec, x, y).transpose()
(Pdb) return
--Return--
> /home/lindsayad/jupyter_notebooks/moose_calc_routines.py(27)strain_rate()->Matrix([
[   ...s(0.3*pi*y)]])
-> return gradVec2(u_vec, x, y) + gradVec2(u_vec, x, y).transpose()
(Pdb) quit






    



---------------------------------------------------------------------------
BdbQuit                                   Traceback (most recent call last)
<ipython-input-1-01e77c4bb29f> in <module>()
     17 
     18 surface_terms = bc_terms_momentum_traction(uvec, n, p, k, eps, x, y, symbolic=False, parts=True)
---> 19 tested_bc = wall_function_momentum_traction(uvec, n, p, k, eps, x, y, "kin", symbolic=False, parts=True)

/home/lindsayad/jupyter_notebooks/moose_calc_routines.py in wall_function_momentum_traction(uvec, nvec, p, k, eps, x, y, tau_type, symbolic, parts)
    144         uTau = uvec_norm / yStarPlus
    145     elif tau_type == "kin":
--> 146         uTau = cmu**.25 * sp.sqrt(k)
    147     else:
    148         raise ValueError("Must either pass 'vel' or 'kin' for tau_type")

/home/lindsayad/jupyter_notebooks/moose_calc_routines.py in strain_rate(u_vec, x, y)
     25 
     26 def strain_rate(u_vec, x, y):
---> 27     return gradVec2(u_vec, x, y) + gradVec2(u_vec, x, y).transpose()
     28 
     29 def strain_rate_squared_2(u_vec, x, y):

/home/lindsayad/miniconda3/lib/python3.6/bdb.py in trace_dispatch(self, frame, event, arg)
     50             return self.dispatch_call(frame, arg)
     51         if event == 'return':
---> 52             return self.dispatch_return(frame, arg)
     53         if event == 'exception':
     54             return self.dispatch_exception(frame, arg)

/home/lindsayad/miniconda3/lib/python3.6/bdb.py in dispatch_return(self, frame, arg)
     94             finally:
     95                 self.frame_returning = None
---> 96             if self.quitting: raise BdbQuit
     97             # The user issued a 'next' or 'until' command.
     98             if self.stopframe is frame and self.stoplineno != -1:

BdbQuit:



In [2]:

    
needed_func = tested_bc - surface_terms
needed_func
print(prep_moose_input(needed_func[0]))









    Out[2]:





$$\left[\begin{matrix}\mu \left(0.14 \pi x \cos{\left (0.2 \pi x y \right )} + 0.06 \pi y \cos{\left (0.3 \pi x y \right )} + 0.48 \pi \cos{\left (0.8 \pi x \right )} + 0.4 \pi \cos{\left (\pi y \right )}\right) + \frac{0.09 \rho \left(0.4 \sin{\left (0.7 \pi x \right )} + 0.9 \sin{\left (0.7 \pi y \right )} + 0.7 \sin{\left (0.4 \pi x y \right )} + 0.4\right)^{2}}{0.6 \sin{\left (0.3 \pi x \right )} + 0.9 \sin{\left (0.9 \pi y \right )} + 0.8 \sin{\left (0.6 \pi x y \right )} + 0.5} \left(0.14 \pi x \cos{\left (0.2 \pi x y \right )} + 0.06 \pi y \cos{\left (0.3 \pi x y \right )} + 0.48 \pi \cos{\left (0.8 \pi x \right )} + 0.4 \pi \cos{\left (\pi y \right )}\right) - 0.14 \pi x \cos{\left (0.2 \pi x y \right )} - 0.06 \pi y \cos{\left (0.3 \pi x y \right )} + 0.497929597731969 \left(0.4 \sin{\left (0.5 \pi x \right )} + 0.4 \sin{\left (\pi y \right )} + 0.7 \sin{\left (0.2 \pi x y \right )} + 0.5\right) \sqrt{0.4 \sin{\left (0.7 \pi x \right )} + 0.9 \sin{\left (0.7 \pi y \right )} + 0.7 \sin{\left (0.4 \pi x y \right )} + 0.4} - 0.48 \pi \cos{\left (0.8 \pi x \right )} - 0.4 \pi \cos{\left (\pi y \right )}\\\frac{0.09 \rho \left(0.12 \pi x \cos{\left (0.3 \pi x y \right )} + 0.18 \pi \cos{\left (0.3 \pi y \right )}\right) \left(0.4 \sin{\left (0.7 \pi x \right )} + 0.9 \sin{\left (0.7 \pi y \right )} + 0.7 \sin{\left (0.4 \pi x y \right )} + 0.4\right)^{2}}{0.6 \sin{\left (0.3 \pi x \right )} + 0.9 \sin{\left (0.9 \pi y \right )} + 0.8 \sin{\left (0.6 \pi x y \right )} + 0.5} - 0.12 \pi x \cos{\left (0.3 \pi x y \right )} + 0.497929597731969 \sqrt{0.4 \sin{\left (0.7 \pi x \right )} + 0.9 \sin{\left (0.7 \pi y \right )} + 0.7 \sin{\left (0.4 \pi x y \right )} + 0.4} \left(0.6 \sin{\left (0.8 \pi x \right )} + 0.3 \sin{\left (0.3 \pi y \right )} + 0.2 \sin{\left (0.3 \pi x y \right )} + 0.3\right) - 0.18 \pi \cos{\left (0.3 \pi y \right )}\end{matrix}\right]$$






    



${mu}*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y)) + 0.09*${rho}*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)^2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5) - 0.14*pi*x*cos(0.2*pi*x*y) - 0.06*pi*y*cos(0.3*pi*x*y) + 0.497929597731969*(0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5)*sqrt(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4) - 0.48*pi*cos(0.8*pi*x) - 0.4*pi*cos(pi*y)



In [3]:

    
print(prep_moose_input(-surface_terms[0]))









    



${mu}*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y)) + 0.09*${rho}*(0.14*pi*x*cos(0.2*pi*x*y) + 0.06*pi*y*cos(0.3*pi*x*y) + 0.48*pi*cos(0.8*pi*x) + 0.4*pi*cos(pi*y))*(0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4)^2/(0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5)



In [3]:

    
wall_function_momentum_traction(uvec, n, p, k, eps, x, y, "kin", symbolic=True, parts=True)









    Out[3]:





$$\left[\begin{matrix}- 0.14 \pi x \cos{\left (0.2 \pi x y \right )} - 0.06 \pi y \cos{\left (0.3 \pi x y \right )} - 0.48 \pi \cos{\left (0.8 \pi x \right )} - 0.4 \pi \cos{\left (\pi y \right )}\\- \mu \left(0.12 \pi x \cos{\left (0.3 \pi x y \right )} + 0.18 \pi \cos{\left (0.3 \pi y \right )}\right) - 0.12 \pi x \cos{\left (0.3 \pi x y \right )} + 0.5 \sin{\left (0.5 \pi x \right )} + 1.0 \sin{\left (0.3 \pi y \right )} + 0.5 \sin{\left (0.2 \pi x y \right )} - 0.18 \pi \cos{\left (0.3 \pi y \right )} + 0.5\end{matrix}\right]$$



In [1]:

    
surface_terms = bc_terms_diffusion(u, n, x, y)
tested_bc = vacuum(u, n)



In [2]:

    
surface_terms
tested_bc









    Out[2]:





$$- 0.14 \pi x \cos{\left (0.2 \pi x y \right )} - 0.4 \pi \cos{\left (\pi y \right )}$$






    Out[2]:





$$0.2 \sin{\left (0.5 \pi x \right )} + 0.2 \sin{\left (\pi y \right )} + 0.35 \sin{\left (0.2 \pi x y \right )} + 0.25$$



In [3]:

    
needed_func = tested_bc - surface_terms



In [6]:

    
needed_func
print(needed_func)









    Out[6]:





$$0.14 \pi x \cos{\left (0.2 \pi x y \right )} + 0.2 \sin{\left (0.5 \pi x \right )} + 0.2 \sin{\left (\pi y \right )} + 0.35 \sin{\left (0.2 \pi x y \right )} + 0.4 \pi \cos{\left (\pi y \right )} + 0.25$$






    



0.14*pi*x*cos(0.2*pi*x*y) + 0.2*sin(0.5*pi*x) + 0.2*sin(pi*y) + 0.35*sin(0.2*pi*x*y) + 0.4*pi*cos(pi*y) + 0.25

5/16/17



In [1]:

    
from moose_calc_routines import *
from sympy import *
import sympy as sp
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
init_printing()

x, y = var('x y')

# INS turbulence
u = 0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5
v = 0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3
p = 0.5*sin(0.5*pi*x) + 1.0*sin(0.3*pi*y) + 0.5*sin(0.2*pi*x*y) + 0.5
k = 0.4*sin(0.7*pi*x) + 0.9*sin(0.7*pi*y) + 0.7*sin(0.4*pi*x*y) + 0.4
eps = 0.6*sin(0.3*pi*x) + 0.9*sin(0.9*pi*y) + 0.8*sin(0.6*pi*x*y) + 0.5

# # INS only
# u = 0.4*sin(0.5*pi*x) + 0.4*sin(pi*y) + 0.7*sin(0.2*pi*x*y) + 0.5
# v = 0.6*sin(0.8*pi*x) + 0.3*sin(0.3*pi*y) + 0.2*sin(0.3*pi*x*y) + 0.3
# p = 0.5*sin(0.5*pi*x) + 1.0*sin(0.3*pi*y) + 0.5*sin(0.2*pi*x*y) + 0.5

uvec = sp.Matrix([u, v])
nvecs = {'left' : sp.Matrix([-1, 0]), 'top' : sp.Matrix([0, 1]), \
         'right' : sp.Matrix([1, 0]), 'bottom' : sp.Matrix([0, -1])}

source = {bnd_name : 
          prep_moose_input(#ins_epsilon_wall_function_bc(nvec, k, eps, x, y)
                          -bc_terms_eps(nvec, k, eps, x, y)[0,0])
          for bnd_name, nvec in nvecs.items()}
# anti_bounds = {'left' : 'top right bottom', 'top' : 'right bottom left',
#                'right' : 'bottom left top', 'bottom' : 'left top right'}
anti_bounds = {'left' : 'top right bottom left', 'top' : 'right bottom left top',
               'right' : 'bottom left top right', 'bottom' : 'left top right bottom'}
h_list = ['5', '10']
base = "k_epsilon_general_bc"
h_array = np.array([.2, .1])
volume_source = {'u' : prep_moose_input(L_momentum_traction(uvec, p, k, eps, x, y)[0]),
                'v' : prep_moose_input(L_momentum_traction(uvec, p, k, eps, x, y)[1]),
                'p' : prep_moose_input(L_pressure(uvec, x, y)),
                'k' : prep_moose_input(L_kin(uvec, k, eps, x, y)),
                'eps' : prep_moose_input(L_eps(uvec, k , eps, x, y))}
diri_func = {'u' : u, 'v' : v, 'p' : p, 'k' : k, 'eps' : eps}



In [1]:

    
a_string = "b"
a_string += "c"
a_string









    Out[1]:





'bc'



In [2]:

    
"a" + None









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-b1f8c2ce73bc> in <module>()
----> 1 "a" + None

TypeError: must be str, not NoneType



In [6]:

    
optional_save_string="epsilon_wall_func_natural"

plot_order_accuracy('left', h_array, base, optional_save_string=optional_save_string)
plot_order_accuracy('right', h_array, base, optional_save_string=optional_save_string)
plot_order_accuracy('top', h_array, base, optional_save_string=optional_save_string)
plot_order_accuracy('bottom', h_array, base, optional_save_string=optional_save_string)

Tasks accomplished today:

Showed formal order accuracy of natural boundary condition using MMS with pure navier stokes
Showed grid convergence with "kinetic" branch of INSMomentumShearStressWallFunctionBC but with accuracy order between formal order and formal order - 1 for u, v, and p for top and bottom boundaries. Unable to solve for left and right boundaries. Formal order accuracy for $\epsilon$ and k for solved cases.
Showed grid convergence with "velocity" branch of INSMomentumShearStressWallFunctionBC but with accuracy order between formal order and formal order - 1 for top and bottom boundaries; formal order accuracy for $\epsilon$ and k for top and bottom boundaries. Between formal order - 1 and formal order - 2 for p, u, and v for left boundary; between formal order and formal order - 1 for $\epsilon$ and k for left boundary. Unable to solve for right boundary.
Demonstrated that just by introducing a small error in the MOOSE code (multiplying a term by 1.1), we can destroy grid convergence by two orders. This makes me feel better about the fact that we're not achieving the exact formal order of accuracy with the INSMomentumShearStressWallFunctionBC but we're still within an order of the formal order.

Natural results suggest the moose python calculation routine for the integrated by parts terms is wrong!



In [4]:

    
string = "Functions" + str('u')
string









    Out[4]:





'Functionsu'



In [5]:

    
import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 2.1, .1)
u = 2*x**2 + 4*x
v = 3*x**2 + 2*x + 1
plt.plot(x, u)
plt.plot(x, v)
plt.show()



In [6]:

    
x = np.arange(0, 2.1, .1)
u = 1 - 2*x + 2*x**2
v = x**2
plt.plot(x, u)
plt.plot(x, v)
plt.show()



In [ ]: