#%matplotlib inline
from IPython.display import display
from sympy.interactive import printing
printing.init_printing(use_latex='mathjax')
from sympy import *

x, y    = symbols("x y")
t       = symbols('t')
e       = symbols("e")
Gamma   = symbols('gamma')
R       = symbols('R')
Cv      = symbols('Cv')
mu      = symbols('mu')
mut     = symbols('mu_t')

# SA variables and constants
chi     = symbols('chi')
fv1     = symbols('fv1')
fn      = symbols('fn')
Pr      = symbols('Pr')
Prt     = symbols('Prt')
sig     = symbols('sigma')
cn1     = symbols('c_n1')
cv1     = symbols('c_v1')
cv2     = symbols('c_v2')
cv3     = symbols('c_v3')
ct3     = symbols('c_t3')
ct4     = symbols('c_t4')
rlim    = symbols('r_lim')
cb1     = symbols('c_b1')   
cb2     = symbols('c_b2')   
cw1     = symbols('c_w1')   
cw2     = symbols('c_w2')   
cw3     = symbols('c_w3')   
K       = symbols('kappa')  
d       = symbols('d')  

k       = symbols('k')
kt      = symbols('k_t')
x0, x1, x2, x3, x4, F = symbols("x0 x1 x2 x3 x4 F")

# intermediary vars for derivatives in the forcing functions
drhodx   = symbols('dRHOdx')
drhody   = symbols('dRHOdy')

dudx     = symbols('dUdx')
dudy     = symbols('dUdy')
d2udxx   = symbols('d2Udxx')
d2udxy   = symbols('d2Udxy')
d2udyy   = symbols('d2Udyy')

dvdx     = symbols('dVdx')
dvdy     = symbols('dVdy')
d2vdxx   = symbols('d2Vdxx')
d2vdxy   = symbols('d2Vdxy')
d2vdyy   = symbols('d2Vdyy')

dEdx     = symbols('dEdx')
dEdy     = symbols('dEdy')
d2Edxx   = symbols('d2Edxx')
d2Edxy   = symbols('d2Edxy')
d2Edyy   = symbols('d2Edyy')

dpdx     = symbols('dPdx')
dpdy     = symbols('dPdy')

dntldx   = symbols('dNTLdx')
dntldy   = symbols('dNTLdy')
d2ntldxx = symbols('d2NTLdxx')
d2ntldxy = symbols('d2NTLdxy')
d2ntldyy = symbols('d2NTLdyy')

# primitive vars
rho      = Function('rho')(x,y)
u        = Function('u')(x,y)
v        = Function('v')(x,y)
w        = Function('w')(x,y)
#E is the total energy per mass
E        = Function('E')(x,y)
ntl      = Function('ntl')(x,y)
p        = Function('p')(x,y)

MS_RHO = Function('MS_RHO')(x,y)
MS_U   = Function('MS_U')(x,y)
MS_V   = Function('MS_V')(x,y)
MS_P   = Function('MS_P')(x,y)
#computing MS_E as a function of MS_P, MS_U, MS_V and their derivatives so it makes shorter expressions
MS_E   = MS_P/((Gamma-1)*MS_RHO)+0.5*(MS_U*MS_U+MS_V*MS_V);

"""
Explanation: Calculation of manufactured solution and forcing functions for the verification of 2D compressible and high-order Euler/NS/RANS-SA solvers via the method of manufactured solutions
by Farshad Navah
farshad.navah at mail.mcgill.ca
source:
https://github.com/fanav/Verification_MMS_Tools
REFERENCE: 
"A comprehensive high-order solver verification methodology for free fluid flows"
Farshad Navah, McGill University; Sivakumaran Nadarajah, McGill University
Read More: https://arxiv.org/abs/1712.09478
arXiv 1712.09478
Preambule
End of explanation
"""


# Function for cleaning up the C output
def Cout_cleanup(s):         
    s = s.replace("rho(x, y)","RHO")
    s = s.replace("u(x, y)","U")
    s = s.replace("v(x, y)","V")
    s = s.replace("W(x, y)","W")    
    s = s.replace("E(x, y)","E")
    s = s.replace("p(x, y)","P")    
    s = s.replace("ntl(x, y)","NTL")    
         
    s = s.replace("M_PI","piMS")
    s = s.replace("gamma","Gamma")
    s = s.replace("Min","MIN")  
    
    s = s.replace("Derivative(MS_RHO(x, y), x)",   "dRHOdx")
    s = s.replace("Derivative(MS_RHO(x, y), x, x)","d2RHOdxx")
    s = s.replace("Derivative(MS_RHO(x, y), y)",   "dRHOdy")    
    s = s.replace("Derivative(MS_RHO(x, y), y, y)","d2RHOdyy")    
    s = s.replace("MS_RHO(x, y)","RHO")     
    
    s = s.replace("Derivative(MS_P(x, y), x)",     "dPdx")
    s = s.replace("Derivative(MS_P(x, y), x, x)",  "d2Pdxx")
    s = s.replace("Derivative(MS_P(x, y), y)",     "dPdy")    
    s = s.replace("Derivative(MS_P(x, y), y, y)",  "d2Pdyy")    
    s = s.replace("MS_P(x, y)","P")    
    
    s = s.replace("Derivative(MS_U(x, y), x)",     "dUdx")    
    s = s.replace("Derivative(MS_U(x, y), x, x)",  "d2Udxx")    
    s = s.replace("Derivative(MS_U(x, y), y)",     "dUdy")    
    s = s.replace("Derivative(MS_U(x, y), y, y)",  "d2Udyy")    
    s = s.replace("MS_U(x, y)","U")        
    
    s = s.replace("Derivative(MS_V(x, y), x)",     "dVdx")    
    s = s.replace("Derivative(MS_V(x, y), x, x)",  "d2Vdxx")        
    s = s.replace("Derivative(MS_V(x, y), y)",     "dVdy")    
    s = s.replace("Derivative(MS_V(x, y), y, y)",  "d2Vdyy")            
    s = s.replace("MS_V(x, y)","V")         
    return(s)

# Functions to substitute derivative of primitive variables by intermediate variable so the final expression is simpler
#                                                                         (broken to 2 expressions dwdx and S(w,dwdx))

def subs_derivs_by_dwdx(var):
    var = var.subs(diff(rho,x)  ,drhodx) 
    var = var.subs(diff(u  ,x)  ,dudx)    
    var = var.subs(diff(v  ,x)  ,dvdx)    
    var = var.subs(diff(E  ,x)  ,dEdx)
    var = var.subs(diff(p  ,x)  ,dpdx)    
    var = var.subs(diff(ntl,x)  ,dntldx)    
    return(var)

def subs_derivs_by_dwdy(var):
    var = var.subs(diff(rho,y)  ,drhody)    
    var = var.subs(diff(u  ,y)  ,dudy)    
    var = var.subs(diff(v  ,y)  ,dvdy)    
    var = var.subs(diff(E  ,y)  ,dEdy)
    var = var.subs(diff(p  ,y)  ,dpdy)    
    var = var.subs(diff(ntl,y)  ,dntldy)     
    return(var)

def subs_derivs_by_d2wdxy(var):
    var = var.subs(diff(u  ,x,x),d2udxx)
    var = var.subs(diff(u  ,x,y),d2udxy)
    var = var.subs(diff(v  ,x,x),d2vdxx)
    var = var.subs(diff(v  ,x,y),d2vdxy)
    var = var.subs(diff(E  ,x,x),d2Edxx)
    var = var.subs(diff(E  ,x,y),d2Edxy)
    var = var.subs(diff(ntl,x,x),d2ntldxx)
    var = var.subs(diff(ntl,x,y),d2ntldxy)
    
    var = var.subs(diff(u  ,y,x),d2udxy)
    var = var.subs(diff(u  ,y,y),d2udyy)
    var = var.subs(diff(v  ,y,x),d2vdxy)
    var = var.subs(diff(v  ,y,y),d2vdyy)
    var = var.subs(diff(E  ,y,x),d2Edxy)
    var = var.subs(diff(E  ,y,y),d2Edyy)
    var = var.subs(diff(ntl,y,x),d2ntldxy)    
    var = var.subs(diff(ntl,y,y),d2ntldyy)    
    return(var)


"""
Explanation: Useful Functions
End of explanation
"""


MS_no = 5

#1 : Inviscid flows in subsonic regime - MS-1
if MS_no==1: 
    NS = 0   
#2 : Inviscid flows in supersonic regime - MS-2    
if MS_no==2: 
    NS = 0    
#3 : Laminar flows - MS-3    
if MS_no==3: 
    NS = 1 
#4 : Turbulent flows - MS-4 (original SA)    
if MS_no==4: 
    NS      =  2
    SA      = +1
    SA_Sbar = +1
    lam_sup =  1
#5 : Turbulent flows - MS-5 (negative SA)    
if MS_no==5: 
    NS      =  2    
    SA      = -1    
    SA_Sbar = +1
    lam_sup =  1
        
# Legend:    
# ----------------------------------------------------------------------------
# NS::      (NS = 0) Euler || (NS = 1) NS || (NS = 2) RANS
# SA::      (SA = -1) Negative SA and ntl<0 || (SA = +1) Original SA and 0<=ntl
# lam_sup:: (lam_sup = 0) No laminar suppression term || (lam_sup = 1) Laminar suppression considered (ft2)
# SA_Sbar:: (SA_Sbar = -1) Negative SA  || (SA_Sbar = +1): Original SA


dim  = 2 # number of dimensions
if NS<=1:
    neq = dim+2
elif NS==2:
    neq = dim+3

"""
Explanation: SWITCHES
Choosing the MS
End of explanation
"""


Cp = R*Gamma/(Gamma-1)
Cv = R/(Gamma-1)

W = Matrix([rho, rho*u, rho*v, rho*E, rho*ntl])

e = simplify(E)-simplify(1/2*(u**2+v**2))
e = simplify(e)

"""
Explanation: Some fundamental definitions
End of explanation
"""


inv_F_orig=Matrix([ [rho*u,       rho*v       ], 
                    [rho*u**2+p,  rho*u*v     ], 
                    [rho*v*u,     rho*v**2+p  ],
                    [(rho*E+p)*u, (rho*E+p)*v ],
                    [(rho*ntl)*u, (rho*ntl)*v ],])
inv_F = inv_F_orig

"""
Explanation: Inviscid Fluxes
End of explanation
"""


#Divergence is not yet developed in Sympy. I have to compute it the long way:
Div_F_inv    = zeros(neq,1)
for i in range(0,neq):
    Div_F_inv[i] = diff(inv_F_orig[i,0],x)+diff(inv_F_orig[i,1],y)

Div_F_inv    = simplify(Div_F_inv)

Div_F_inv = subs_derivs_by_dwdx(Div_F_inv)
Div_F_inv = subs_derivs_by_dwdy(Div_F_inv)

"""
Explanation: Computing the Inviscid Flux divergence:
End of explanation
"""


# conductivity
k  = Cp*mu/Pr

"""
Explanation: Viscous Fluxes
End of explanation
"""


if (NS>=2):
    nu   = (1/W[0])*mu
    chi  = ntl/nu
    fv1  = pow(chi,3)/(pow(chi,3)+pow(cv1,3))
    mut  = rho*ntl*fv1
    kt   = Cp*mut/Prt

"""
Explanation: Spalart-Allmaras variables
End of explanation
"""


# Depending on the sign of ntl, there are two possibilities:
if (NS==2):
    if (SA>0):
        fn = 1
    else:
    #for ntl<0
        fn = (cn1+pow(chi,3))/(cn1-pow(chi,3))

"""
Explanation: negative SA's 'Fn' function
End of explanation
"""


mu_eff = 0
k_eff  = 0
# if NS
if (NS==1):
    mu_eff = mu
    k_eff  = k 

# if RANS-SA:
elif (NS==2):
    # negative SA decouples RANS and SA one-way
    if(SA<0):
        mu_eff = mu
        k_eff  = k
    else:
        mu_eff  = mu + mut
        k_eff   = k  + kt

"""
Explanation: Effective viscosity and conductivity (depend on the PDE being solved)
End of explanation
"""


Tau = zeros(dim,dim)
if NS>0:
    Tau = 2*mu_eff*Matrix([ [diff(u,x)-1/3*(diff(u,x)+diff(v,y)), 1/2*(diff(v,x)+diff(u,y))            ],
                          [1/2*(diff(u,y)+diff(v,x))            , diff(v,y)-1/3*(diff(u,x)+diff(v,y))  ]])
    Tau = simplify(Tau)
    Tau = nsimplify(Tau)
    Tau

"""
Explanation: The Shear Stress definition
End of explanation
"""


if NS>0:
    T = e/Cv
    dTdx = simplify(diff(T,x))
    dTdy = simplify(diff(T,y))

"""
Explanation: The Temperature and its derivatives
End of explanation
"""


if NS>0:
    q = simplify(k_eff*Matrix([dTdx, dTdy]))

"""
Explanation: Heat Flux
End of explanation
"""


#Viscous_Flux
if NS>0:
    vis_F_orig = -Matrix([ [0                                     , 0                            ], 
                           [Tau[0,0]                              , Tau[0,1]                     ], 
                           [Tau[1,0]                              , Tau[1,1]                     ],
                           [u*Tau[0,0]+v*Tau[0,1]+q[0]            , u*Tau[1,0]+v*Tau[1,1]+q[1]   ],
                           [(1/sig)*(mu+W[4]*fn)*diff(ntl,x)      ,(1/sig)*(mu+W[4]*fn)*diff(ntl,y)  ], ]) 

if NS>0:
    vis_F = simplify(vis_F_orig)

"""
Explanation: Computing the viscous Flux
End of explanation
"""


#Compute the divergence the long way:
if NS>0:
    Div_F_vis    = zeros(neq,1)
    for i in range(0,neq):
        Div_F_vis[i] = diff(vis_F_orig[i,0],x)+diff(vis_F_orig[i,1],y)
    Div_F_vis    = simplify(Div_F_vis)

if NS>0:
    Div_F_vis = subs_derivs_by_d2wdxy(Div_F_vis)

    Div_F_vis = subs_derivs_by_dwdx(Div_F_vis)
    Div_F_vis = subs_derivs_by_dwdy(Div_F_vis)

"""
Explanation: Computing the Viscous Flux divergence:
End of explanation
"""


def SA_Sources(NS,SA,SA_Sbar,lam_sup):
    if NS==2:
        # interm vars
        dtp2   = d*d
        oosig  = 1/sig      
        Ktp2   = K*K
        cw3tp6 = pow(cw3,6)    
        chitp2 = pow(chi,2)        
        chitp3 = pow(chi,3)

        fv2  = 1-(chi/(1+chi*fv1))
        S    = abs((dvdx-dudy))
        Sbar = fv2*ntl/((Ktp2)*(dtp2))
        #Laminar suppression and trip terms:
        if lam_sup!=0:
            ft2  = ct3*exp(-ct4*chitp2)
            Trp  = 0
            print ("Laminar suppression is considered, but no trip term")
        else:
            ft2  = 0
            Trp  = 0

    # The modified vorticity: Stl (S_tilda)
    # the condition on the S_bar (and modified vorticity) is not directly dependent on the sign of ntl
    # it is rather based on the numerical condition (Sbar>=-cv2*S).

        if SA_Sbar>0:
            Stl = S + Sbar;
        elif SA_Sbar<0:
            Stl = S + (S*(cv2*cv2*S+cv3*Sbar))/(S*(cv3-2*cv2)-Sbar)  

        # The Production term
        if(SA>0):
            P = cb1*Stl*ntl*(1-ft2)        
        elif SA<0:
            P = cb1*S*ntl*(1-ct3)

        # The Destruction term
        if SA>0:
            r    = Min(ntl/(Stl*Ktp2*dtp2),rlim)      
            rtp6 = pow(r,6)
            g    = r + cw2*(rtp6-r)
            gtp6 = pow(g,6)
            fw   = g*pow((1+cw3tp6)/(gtp6+cw3tp6),1/6)
            D    = (cw1*fw-cb1*ft2/(Ktp2))*(ntl*ntl)/(dtp2)
        # Negative SA
        elif SA<0:
            D    = -cw1*(ntl*ntl)/(dtp2)  

        # fn function is defined earlier    

        # The "Distribution" term
        trm1 = cb2/sig*(dntldx*dntldx+dntldy*dntldy);

        # The density counterpart of the conservative form
        trm2 = oosig*(nu+ntl*fn)*(drhodx*dntldx+drhody*dntldy);    

        #src  = rho*(-P+D-Trp)-rho*trm1+trm2;    
        S_SA_src_product = nsimplify(  rho*(-P)   )
        S_SA_src_destruc = nsimplify(  rho*(+D)   )
        S_SA_src_distrib = nsimplify( -rho*(trm1) )
        S_SA_src_conserv = nsimplify(  trm2       )

        S_SA = S_SA_src_product + S_SA_src_destruc + S_SA_src_distrib + S_SA_src_conserv
    else:
        S_SA = 0

    S_SA  = nsimplify(S_SA)  
    
    return(S_SA_src_product,S_SA_src_destruc,S_SA_src_distrib,S_SA_src_conserv,S_SA)    

if NS==2:
    S_SA_src_product, S_SA_src_destruc, S_SA_src_distrib, S_SA_src_conserv, S_SA = SA_Sources(NS,SA,SA_Sbar,lam_sup)

"""
Explanation: SA source terms
End of explanation
"""


S_inv = zeros(neq,1)
S_vis = zeros(neq,1)
S     = zeros(neq,1)

for i in range(0,neq):
    S_inv[i] = Div_F_inv[i]
    if NS>0:
        S_vis[i] = Div_F_vis[i] 
        
    S[i] = S_inv[i]+S_vis[i]   
    
    if NS==2 and i==neq-1:
        S[i] = S[i]+S_SA

"""
Explanation: Computing the Forcing Functions for the MMS (balancing discrete residuals)
End of explanation
"""


s_inv = ccode(S_inv[0],assign_to='FFQ_RHO_INV')
s_inv = Cout_cleanup(s_inv)
print (len(s_inv),"\n")
print (s_inv,"\n\n")

s_vis = ccode(S_vis[0],assign_to='FFQ_RHO_VIS')
s_vis = Cout_cleanup(s_vis)
print (len(s_vis),"\n")
print (s_vis,"\n\n")

print (100*"-")

s_inv = ccode(S_inv[1],assign_to='FFQ_U_INV')
s_inv = Cout_cleanup(s_inv)
print (len(s_inv),"\n")
print (s_inv,"\n\n")

s_vis = ccode(S_vis[1],assign_to='FFQ_U_VIS')
s_vis = Cout_cleanup(s_vis)
print (len(s_vis),"\n")
print (s_vis,"\n\n")

print (100*"-")

s_inv = ccode(S_inv[2],assign_to='FFQ_V_INV')
s_inv = Cout_cleanup(s_inv)
print (len(s_inv),"\n")
print (s_inv,"\n\n")

s_vis = ccode(S_vis[2],assign_to='FFQ_V_VIS')
s_vis = Cout_cleanup(s_vis)
print (len(s_vis),"\n")
print (s_vis,"\n\n")

print (100*"-")

s_inv = ccode(S_inv[3],assign_to='FFQ_E_INV')
s_inv = Cout_cleanup(s_inv)
print (len(s_inv),"\n")
print (s_inv,"\n\n")

s_vis = ccode(S_vis[3],assign_to='FFQ_E_VIS')
s_vis = Cout_cleanup(s_vis)
print (len(s_vis),"\n")
print (s_vis,"\n\n")

print (100*"-")

if NS==2:
    s_inv = ccode(S_inv[4],assign_to='FFQ_NTL_INV')
    s_inv = Cout_cleanup(s_inv)    
    print (len(s_inv),"\n")
    print (s_inv,"\n\n")

if NS==2:
    s_vis = ccode(S_vis[4],assign_to='FFQ_NTL_VIS')
    s_vis = Cout_cleanup(s_vis)    
    print (len(s_vis),"\n")
    print (s_vis,"\n\n")

if NS==2:
    s_vis = ccode(S_SA_src_product,assign_to='FFQ_NTL_SRC_PRODUCT')
    s_vis = Cout_cleanup(s_vis)    
    print (len(s_vis),"\n")
    print (s_vis,"\n\n")
    
if NS==2:
    s_vis = ccode(S_SA_src_destruc,assign_to='FFQ_NTL_SRC_DESTRUCT')
    s_vis = Cout_cleanup(s_vis)    
    print (len(s_vis),"\n")
    print (s_vis,"\n\n")
    
if NS==2:
    s_vis = ccode(S_SA_src_distrib,assign_to='FFQ_NTL_SRC_DISTRIB')
    s_vis = Cout_cleanup(s_vis)    
    print (len(s_vis),"\n")
    print (s_vis,"\n\n")
    
if NS==2:
    s_vis = ccode(S_SA_src_conserv,assign_to='FFQ_NTL_SRC_CONSERV')
    s_vis = Cout_cleanup(s_vis)    
    print (len(s_vis),"\n")
    print (s_vis,"\n\n")   

"""
Explanation: Printing the forcing function in C
Split source terms into inviscid, viscous and SA source term components
End of explanation
"""


print ("The Manufactured Solution Number -",MS_no, "- is activated")

if MS_no in range(1,6):

    rho_0   = symbols('rho_0')  
    rho_x   = symbols('rho_x')  
    rho_y   = symbols('rho_y')  
    rho_xy  = symbols('rho_xy')     
    u_0     = symbols('u_0')    
    u_x     = symbols('u_x')    
    u_y     = symbols('u_y')    
    u_xy    = symbols('u_xy') 
    v_0     = symbols('v_0')    
    v_x     = symbols('v_x')    
    v_y     = symbols('v_y')  
    v_xy    = symbols('v_xy') 
    p_0     = symbols('p_0')    
    p_x     = symbols('p_x')    
    p_y     = symbols('p_y')    
    p_xy    = symbols('p_xy')

    a_rhox  = symbols('a_rhox') 
    a_rhoy  = symbols('a_rhoy') 
    a_rhoxy = symbols('a_rhoxy') 
    a_ux    = symbols('a_ux')   
    a_uy    = symbols('a_uy')   
    a_uxy   = symbols('a_uxy')   
    a_vx    = symbols('a_vx')   
    a_vy    = symbols('a_vy')   
    a_vxy   = symbols('a_vxy')   
    a_px    = symbols('a_px')   
    a_py    = symbols('a_py')   
    a_pxy   = symbols('a_pxy')       
    L       = symbols('L')

    # SA
    nu_sa_0  =  symbols('nu_sa_0')
    nu_sa_x  =  symbols('nu_sa_x')
    nu_sa_y  =  symbols('nu_sa_y')
    nu_sa_xy =  symbols('nu_sa_xy')
    nu_sa_t  =  symbols('nu_sa_t')
    a_nusay  =  symbols('a_nusay')
    a_nusax  =  symbols('a_nusax')
    a_nusaxy =  symbols('a_nusaxy')
    a_nusat  =  symbols('a_nusat')

    MS_RHO = rho_0   + rho_x   * sin(a_rhox  * pi * x / L)  + rho_y   * cos(a_rhoy  * pi * y / L)   + rho_xy   * cos(a_rhoxy * pi * x / L)*cos(a_rhoxy * pi * y / L);
    MS_U   = u_0     + u_x     * sin(a_ux    * pi * x / L)  + u_y     * cos(a_uy    * pi * y / L)   + u_xy     * cos(a_uxy   * pi * x / L)*cos(a_uxy   * pi * y / L);      
    MS_V   = v_0     + v_x     * cos(a_vx    * pi * x / L)  + v_y     * sin(a_vy    * pi * y / L)   + v_xy     * cos(a_vxy   * pi * x / L)*cos(a_vxy   * pi * y / L);
    MS_P   = p_0     + p_x     * cos(a_px    * pi * x / L)  + p_y     * sin(a_py    * pi * y / L)   + p_xy     * cos(a_pxy   * pi * x / L)*cos(a_pxy   * pi * y / L);
    MS_NTL = nu_sa_0 + nu_sa_x * cos(a_nusax * pi * x / L)  + nu_sa_y * cos(a_nusay * pi * y / L)   + nu_sa_xy * cos(a_nusaxy* pi * x / L)*cos(a_nusaxy* pi * y / L);

"""
Explanation: Manufactured Solutions
Trigonometric MS
End of explanation
"""


#Compute the MS derivatives

dEdx     = diff(MS_E  ,x)
dEdy     = diff(MS_E  ,y)
d2Edxx   = diff(MS_E  ,x,x)
d2Edyy   = diff(MS_E  ,y,y)

dRHOdx   = diff(MS_RHO,x)
dRHOdy   = diff(MS_RHO,y)
d2RHOdxx = diff(MS_RHO,x,x)
d2RHOdyy = diff(MS_RHO,y,y)

dUdx     = diff(MS_U  ,x)
dUdy     = diff(MS_U  ,y)
d2Udxx   = diff(MS_U  ,x,x)
d2Udxy   = diff(MS_U  ,x,y)
d2Udyy   = diff(MS_U  ,y,y)

dVdx     = diff(MS_V  ,x)
dVdy     = diff(MS_V  ,y)
d2Vdxx   = diff(MS_V  ,x,x)
d2Vdxy   = diff(MS_V  ,x,y)
d2Vdyy   = diff(MS_V  ,y,y)

dPdx     = diff(MS_P  ,x)
dPdy     = diff(MS_P  ,y)
d2Pdxx   = diff(MS_P  ,x,x)
d2Pdyy   = diff(MS_P  ,y,y)

if NS==2:
    dNTLdx   = diff(MS_NTL,x)
    dNTLdy   = diff(MS_NTL,y)
    d2NTLdxx = diff(MS_NTL,x,x)    
    d2NTLdyy = diff(MS_NTL,y,y)

# State variables
s = ccode(MS_RHO,assign_to='RHO')
s = Cout_cleanup(s)
print (s,"\n\n")

# First degree derivatives
s = ccode(dRHOdx,assign_to='dRHOdx')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dRHOdy,assign_to='dRHOdy')
s = Cout_cleanup(s)
print (s,"\n\n")

if (NS>0):
    s = ccode(d2RHOdxx,assign_to='d2RHOdxx')
    s = Cout_cleanup(s)
    #print len(s),"\n"
    print (s,"\n\n")

    s = ccode(d2RHOdyy,assign_to='d2RHOdyy')
    s = Cout_cleanup(s)
    #print len(s),"\n"
    print (s,"\n\n")

s = ccode(MS_U,assign_to='U')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dUdx,assign_to='dUdx')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dUdy,assign_to='dUdy')
s = Cout_cleanup(s)
print (s,"\n\n")

if (NS>0):
    s = ccode(d2Udxx,assign_to='d2Udxx')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(d2Udxy,assign_to='d2Udxy')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(d2Udyy,assign_to='d2Udyy')
    s = Cout_cleanup(s)
    print (s,"\n\n")

s = ccode(MS_V,assign_to='V')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dVdx,assign_to='dVdx')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dVdy,assign_to='dVdy')
s = Cout_cleanup(s)
print (s,"\n\n")

if (NS>0):
    s = ccode(d2Vdxx,assign_to='d2Vdxx')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(d2Vdxy,assign_to='d2Vdxy')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(d2Vdyy,assign_to='d2Vdyy')
    s = Cout_cleanup(s)
    print (s,"\n\n")

s = ccode(MS_E,assign_to='E')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dEdx,assign_to='dEdx')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dEdy,assign_to='dEdy')
s = Cout_cleanup(s)
print (s,"\n\n")

if (NS>0):
    s = ccode(d2Edxx,assign_to='d2Edxx')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(d2Edyy,assign_to='d2Edyy')
    s = Cout_cleanup(s)
    print (s,"\n\n")

if NS==2:
    s = ccode(MS_NTL,assign_to='NTL')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(dNTLdx,assign_to='dNTLdx')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(dNTLdy,assign_to='dNTLdy')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(d2NTLdxx,assign_to='d2NTLdxx')
    s = Cout_cleanup(s)
    print (s,"\n\n")

    s = ccode(d2NTLdyy,assign_to='d2NTLdyy')
    s = Cout_cleanup(s)
    print (s,"\n\n")
    
s = ccode(MS_P,assign_to='P')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dPdx,assign_to='dPdx')
s = Cout_cleanup(s)
print (s,"\n\n")

s = ccode(dPdy,assign_to='dPdy')
s = Cout_cleanup(s)
print (s,"\n\n")

if (NS>0):
    s = ccode(d2Pdxx,assign_to='d2Pdxx')
    s = Cout_cleanup(s)
    #print len(s),"\n"
    print (s,"\n\n")

    s = ccode(d2Pdyy,assign_to='d2Pdyy')
    s = Cout_cleanup(s)
    #print len(s),"\n"
    print (s,"\n\n")    

"""
Explanation: Printing code in C for the derivatives of the primitive variables
End of explanation
"""


