Ondrej Certik wrote:
> On Mon, Aug 17, 2009 at 11:06 AM, William Stein<wst...@gmail.com> wrote:
>
>> On Mon, Aug 17, 2009 at 6:27 AM, Mani chandra<mchan...@iitk.ac.in> wrote:
>>
>>> William Stein wrote:
>>>
>>>> On Mon, Aug 17, 2009 at 2:26 AM, Juan Jose
>>>> Garcia-Ripoll<juanjose.garciarip...@googlemail.com> wrote:
>>>>
>>>>
>>>>> On Mon, Aug 17, 2009 at 10:57 AM, William Stein<wst...@gmail.com> wrote:
>>>>>
>>>>>
>>>>>> Note that Sage's Maxima uses ECL. So the basic question is, how can
>>>>>> we increase the memory that Maxima + ECL can use?
>>>>>>
>>>>>>
>>>>> The limits ECL has are by default too small for big applications, but
>>>>> it is intentionally done so. However, changing them is pretty easy:
>>>>> add a call to ext:set-limit in any file of maxima that forms part of
>>>>> the final executable.
>>>>>
>>>>> The different memory limits that can be independently controlled are
>>>>> listed here
>>>>> http://ecls.sourceforge.net/new-manual/re34.html
>>>>> So for instance, in your case, which hits the dynamically allocated
>>>>> memory limit, you might add the following
>>>>> (ext:set-limit 'ext:heap-size (* 1024 1024 1024))
>>>>> to maxima/src/ecl-port.lisp in order get 1GB memory limit.
>>>>>
>>>>> Juanjo
>>>>>
>>>>>
>>>> Thanks!!! I've made this trac #6772:
>>>>
>>>> http://trac.sagemath.org/sage_trac/ticket/6772
>>>>
>>>> -- William
>>>>
>>>>
>>>>
>>> Hi,
>>>
>>> I'm glad that this issue has finally been taken up. But as an interm
>>> measure, how does one pass this command through a sage code? I tried
>>> maxima.eval_String, but it doesn't seem to work.
>>>
>> My understanding is that one has to modify the maxima source code
>> itself and recompile maxima.
>>
>>
>>> In an unrelated note, I rewrote my code trying to use a standalone
>>> installation of sympy and it takes forever to compute what maxima in
>>> sage did in a few mins, so I guess sympy is nowhere close to "industrial
>>> strength".
>>>
>>> Mani chandra
>>>
>>>
>> You should post your program to the sympy list and complain. I'm sure
>> they would love to fix whatever is making it so slow. A lot of people
>> are actively working on sympy this summer.
>>
>
> Could you please post your program to us? I just spent 3 months
> working on MHD in Los Alamos and as a SymPy developer I would love to
> see your real life application and what makes it slow and fix it.
>
> Thanks,
> Ondrej
>
> >
>
>
Hi,
I've attached the code. It basically constructs a set of coupled
equations for each velocity and magnetic field mode, which I can then
analyse for various non-linear phenomenon. I'm not fully done with the
code yet, however most of it is there. The time consuming step is the
select_mode() function, which does an inner product of an expression
with a chosen mode and the picks up the co-efficients of that mode.
Regards,
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
-~----------~----~----~----~------~----~------~--~---
from sympy import *
modes = []; N = 1
# 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."
