"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
In [1]:
#%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);
In [2]:
# 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)
Choosing the MS
In [3]:
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
In [4]:
dim = 2 # number of dimensions
if NS<=1:
neq = dim+2
elif NS==2:
neq = dim+3
In [5]:
Cp = R*Gamma/(Gamma-1)
Cv = R/(Gamma-1)
In [6]:
W = Matrix([rho, rho*u, rho*v, rho*E, rho*ntl])
In [7]:
e = simplify(E)-simplify(1/2*(u**2+v**2))
e = simplify(e)
In [8]:
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
In [9]:
#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)
In [10]:
Div_F_inv = subs_derivs_by_dwdx(Div_F_inv)
Div_F_inv = subs_derivs_by_dwdy(Div_F_inv)
In [11]:
# conductivity
k = Cp*mu/Pr
In [12]:
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
In [13]:
# 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))
In [14]:
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
In [15]:
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
In [16]:
if NS>0:
T = e/Cv
dTdx = simplify(diff(T,x))
dTdy = simplify(diff(T,y))
In [17]:
if NS>0:
q = simplify(k_eff*Matrix([dTdx, dTdy]))
In [18]:
#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) ], ])
In [19]:
if NS>0:
vis_F = simplify(vis_F_orig)
In [20]:
#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)
In [21]:
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)
In [22]:
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)
In [23]:
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)
In [24]:
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
In [25]:
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")
In [26]:
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);
In [27]:
#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)
In [28]:
# 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")
In [28]: