Hi,
I wrote a code in sage to construct a Low Dimensional model for the
Magnetohydrodynamic equations to study dynamos etc... But the program
crashes at different points on different machines with the following error:
Traceback (most recent call last):
File "./LowD_model.py", line 158, in <module>
du_x_dt.append(select_mode(RHS_Ux, l, m, n) )
File "./LowD_model.py", line 28, in select_mode
val = (func*exp(-I*(l*x + m*y +
n*z))).integrate(x,0,2*pi).integrate(y,0,2*pi).integrate(z,0,2*pi)
File "expression.pyx", line 5700, in
sage.symbolic.expression.Expression.integral
(sage/symbolic/expression.cpp:24436)
File
"/opt/sage/sage/local/lib/python2.6/site-packages/sage/calculus/calculus.py",
line 566, in integral
result = expression._maxima_().integrate(v, a, b)
File
"/opt/sage/sage/local/lib/python2.6/site-packages/sage/interfaces/maxima.py",
line 2003, in integral
return I(var, min, max)
File
"/opt/sage/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py",
line 1382, in __call__
return self._obj.parent().function_call(self._name, [self._obj] +
list(args), kwds)
File
"/opt/sage/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py",
line 1290, in function_call
return self.new(s)
File
"/opt/sage/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py",
line 1086, in new
return self(code)
File
"/opt/sage/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py",
line 1021, in __call__
return cls(self, x, name=name)
File
"/opt/sage/sage/local/lib/python2.6/site-packages/sage/interfaces/expect.py",
line 1425, in __init__
raise TypeError, x
TypeError: Error executing code in Maxima
CODE:
sage2436 : integrate(sage2432,sage2433,sage2434,sage2435)$
Maxima ERROR:
Maxima encountered a Lisp error:
Memory limit reached. Please jump to an outer point or quit program.
Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.
I'm attaching the entire code but the function where it crashed is the
follows:
def select_mode(func,l,m,n):
val = (func*exp(-I*(l*x + m*y +
n*z))).integrate(x,0,2*pi).integrate(y,0,2*pi).integrate(z,0,2*pi)
return val/(8*pi**3)
Also, I'd like to know if the symbolic computation backend for SAGE is
Maxima, Pynac or Sympy? Or is it a combination of all the three? Will
SAGE move to a single backend in the future? If this bug can't be fixed,
is it suggested that I rewrite the code using the Sympy library?
P.S If anyone is running the code, then please note that it takes a LOT
of time.
Thanking you,
Mani chandra
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
sage-support-unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/sage-support
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
import sys
from sage.all import *
modes = []; N = 2
# Choose the modes for the model
for i in range(-N,N+1):
for j in range(-N,N+1):
for k in range(-N,N+1):
modes.append([i,j,k])
print "Number of modes = ", len(modes)
nu, eta = var('nu, eta')
x, y, z = var('x, y, z')
ko = var('ko'); ko = 2
# Taylor-Green forcing
f_x = sin(ko*x)*cos(ko*y)*cos(ko*z)
f_y = -cos(ko*x)*sin(ko*y)*cos(ko*z)
f_z = 0
def basis(l,m,n):
return exp(I*(l*x + m*y + n*z))
def select_mode(func,l,m,n):
val = (func*exp(-I*(l*x + m*y + n*z))).integrate(x,0,2*pi).integrate(y,0,2*pi).integrate(z,0,2*pi)
return val/(8*pi**3)
u_x = []; b_x = []
u_y = []; b_y = []
u_z = []; b_z = []
p = [];
U_x = var('U_x'); U_x = 0
U_y = var('U_y'); U_y = 0
U_z = var('U_z'); U_z = 0
B_x = var('B_x'); B_x = 0
B_y = var('B_y'); B_y = 0
B_z = var('B_z'); B_z = 0
P = var('P'); P = 0
du_x_dt = []; db_x_dt = []
du_y_dt = []; db_y_dt = []
du_z_dt = []; db_z_dt = []
for i in range(len(modes)):
u_x.append(var('u_x' + str(i)))
u_y.append(var('u_y' + str(i)))
u_z.append(var('u_z' + str(i)))
b_x.append(var('b_x' + str(i)))
b_y.append(var('b_y' + str(i)))
b_z.append(var('b_z' + str(i)))
p.append(var('p' + str(i)))
# Solve u_z in terms of u_x and u_y using the divergence free condition.
for i in range(len(modes)):
l = modes[i][0]
m = modes[i][1]
n = modes[i][2]
if (n!=0):
u_z[i] = (-l*u_x[i] - m*u_y[i])/n
b_z[i] = (-l*b_x[i] - m*b_y[i])/n
elif (m!=0 and n==0):
u_z[i] = 0
b_z[i] = 0
u_y[i] = -l*u_x[i]/m
b_y[i] = -l*b_x[i]/m
elif (m==0 and n==0):
u_x[i] = 0
b_x[i] = 0
u_y[i] = 0
b_y[i] = 0
u_z[i] = 0
b_z[i] = 0
for i in range(len(modes)):
l = modes[i][0]
m = modes[i][1]
n = modes[i][2]
U_x = U_x + u_x[i]*basis(l,m,n)
U_y = U_y + u_y[i]*basis(l,m,n)
U_z = U_z + u_z[i]*basis(l,m,n)
B_x = B_x + b_x[i]*basis(l,m,n)
B_y = B_y + b_y[i]*basis(l,m,n)
B_z = B_z + b_z[i]*basis(l,m,n)
print "Computing Laplacian..."
print " "
laplacian_Ux = diff(U_x, x, 2) + diff(U_x, y, 2) + diff(U_x, z, 2)
laplacian_Uy = diff(U_y, x, 2) + diff(U_y, y, 2) + diff(U_y, z, 2)
laplacian_Uz = diff(U_z, x, 2) + diff(U_z, y, 2) + diff(U_z, z, 2)
laplacian_Bx = diff(B_x, x, 2) + diff(B_x, y, 2) + diff(B_x, z, 2)
laplacian_By = diff(B_y, x, 2) + diff(B_y, y, 2) + diff(B_y, z, 2)
laplacian_Bz = diff(B_z, x, 2) + diff(B_z, y, 2) + diff(B_z, z, 2)
print "Laplacian computation complete."
print " "
print "Computing non-linear terms..."
print " "
nlin_Ux_u = U_x*U_x.diff(x) + U_y*U_x.diff(y) + U_z*U_x.diff(z)
nlin_Ux_b = B_x*B_x.diff(x) + B_y*B_x.diff(y) + B_z*B_x.diff(z)
nlin_Uy_u = U_x*U_y.diff(x) + U_y*U_y.diff(y) + U_z*U_y.diff(z)
nlin_Uy_b = B_x*B_y.diff(x) + B_y*B_y.diff(y) + B_z*B_y.diff(z)
nlin_Uz_u = U_x*U_z.diff(x) + U_y*U_z.diff(y) + U_z*U_z.diff(z)
nlin_Uz_b = B_x*B_z.diff(x) + B_y*B_z.diff(y) + B_z*B_z.diff(z)
nlin_Bx_u = B_x*U_x.diff(x) + B_y*U_x.diff(y) + B_z*U_x.diff(z)
nlin_Bx_b = U_x*B_x.diff(x) + U_y*B_x.diff(y) + U_z*B_x.diff(z)
nlin_By_u = B_x*U_y.diff(x) + B_y*U_y.diff(y) + B_z*U_y.diff(z)
nlin_By_b = U_x*B_y.diff(x) + U_y*B_y.diff(y) + U_z*B_y.diff(z)
nlin_Bz_u = B_x*U_z.diff(x) + B_y*U_z.diff(y) + B_z*U_z.diff(z)
nlin_Bz_b = U_x*B_z.diff(x) + U_y*B_z.diff(y) + U_z*B_z.diff(z)
print "Computation of non-linear terms complete."
print " "
# Solve the Poisson equation for pressure.
print "Computing pressure..."
print " "
div_nlin_U = nlin_Ux_u.diff(x) + nlin_Uy_u.diff(y) + nlin_Uz_u.diff(z)
for i in range(len(modes)):
l = modes[i][0]
m = modes[i][0]
n = modes[i][0]
k_sqr = l**2 + m**2 + n**2
print "Computing pressure mode", modes[i]
if (k_sqr!=0):
p[i] = -select_mode(div_nlin_U, l, m, n)/(l**2 + m**2 + n**2)
else:
p[i] = 0
P = P + p[i]
print " "
print "Pressure computation complete."
print " "
RHS_Ux = -nlin_Ux_u + nlin_Ux_b + nu*laplacian_Ux - P.diff(x) + f_x
RHS_Uy = -nlin_Uy_u + nlin_Uy_b + nu*laplacian_Uy - P.diff(y) + f_y
RHS_Uz = -nlin_Uz_u + nlin_Uz_b + nu*laplacian_Uz - P.diff(z) + f_z
RHS_Bx = -nlin_Bx_b + nlin_Bx_u + eta*laplacian_Bx
RHS_By = -nlin_By_b + nlin_By_u + eta*laplacian_By
RHS_Bz = -nlin_Bz_b + nlin_Bz_u + eta*laplacian_Bz
for i in range(len(modes)):
l = modes[i][0]
m = modes[i][1]
n = modes[i][2]
print "Extracting mode ", modes[i]
du_x_dt.append(select_mode(RHS_Ux, l, m, n) )
du_y_dt.append(select_mode(RHS_Uy, l, m, n) )
du_z_dt.append(select_mode(RHS_Uz, l, m, n) )
db_x_dt.append(select_mode(RHS_Bx, l, m, n) )
db_y_dt.append(select_mode(RHS_By, l, m, n) )
db_z_dt.append(select_mode(RHS_Bz, l, m, n) )
print "Low dimensional model construction complete."
