Programming This forum is for all programming questions.
The question does not have to be directly related to Linux and any language is fair game. 
Notices 
Welcome to LinuxQuestions.org, a friendly and active Linux Community.
You are currently viewing LQ as a guest. By joining our community you will have the ability to post topics, receive our newsletter, use the advanced search, subscribe to threads and access many other special features. Registration is quick, simple and absolutely free. Join our community today!
Note that registered members see fewer ads, and ContentLink is completely disabled once you log in.
Are you new to LinuxQuestions.org? Visit the following links:
Site Howto 
Site FAQ 
Sitemap 
Register Now
If you have any problems with the registration process or your account login, please contact us. If you need to reset your password, click here.
Having a problem logging in? Please visit this page to clear all LQrelated cookies.

Introduction to Linux  A Hands on Guide
This guide was created as an overview of the Linux Operating System, geared toward new users as an exploration tour and getting started guide, with exercises at the end of each chapter.
For more advanced trainees it can be a desktop reference, and a collection of the base knowledge needed to proceed with system and network administration. This book contains many real life examples derived from the author's experience as a Linux system and network administrator, trainer and consultant. They hope these examples will help you to get a better understanding of the Linux system and that you feel encouraged to try out things on your own.
Click Here to receive this Complete Guide absolutely free. 


03152011, 08:17 PM

#1

LQ Newbie
Registered: Jul 2010
Posts: 17
Rep:

Python numerical integration, not getting any output.
I'm translating a numerical integration program from Fortran to Python, and I thought that I had gotten rid of all the syntax problems. However, my program runs for a long time without any writing anything into the output file. The output file is created, but there just isn't anything in it. The program has been running for over thirty minutes now. Although it may just be slow, I am new to Python and I am worried that I may have made a mistake somewhere, especially with the output code.
There are some variables that are defined, but not really used. Those are for when the code gets expanded.
Here is the Python code.
Code:
#! /usr/bin/env python
# This program solves the parallel coupled electronic dimer
from numpy import *
import math
# System Parameters
# np = number of LRC circuit units in the chain
# dll = inductance
# dmm = mutual inductance
# dkk = proportionality coefficient. dmm = dkk*dll
# dcc = capacitance
# omega0 = natural frequency of a unit. omega0 ** 2 = 1/(LC)
# GG = gain/loss. Sqrt(L/C)/R
# mu = rescaled mutual induction. dmm / dll
# tmax = maximum time evolution
# dt = integration time step
# ntimes = nsnaps of the number of times we are calling the output
# subroutine in order to write the output data
# init = initial excitation
GG = float(raw_input('Enter the value of gain/loss.\n'))
"""np = int(raw_input('Enter the number of LRC circuits.\n'))
print np
"""
np = 2
init = 1
ntimes = 1000
dll = 1.0
dcc = 1.0
dmm = 0.8
mu =0 #dmm / dll
tmax = 250.0
dt = 0.01
m = 0
dt2 = dt/2.0
ifile = tmax/(ntimes*dt)
ydot = zeros([2,2], dtype=float)
#
# Initial Condition
# Temporary workaround; only for np=2. Not yet general.
# y1 and y2 are the first and second circuits.
# The first and second elements of each array correspond to Y[i,1] and Y[i,2] in the Fortran code.
#
y = array([[1,0],[0,0]])
savey = zeros([np,2], dtype=float)
phi = zeros([np,2], dtype=float)
#z = zeros([np,2], dtype=float)
#zt = zeros([np,2], dtype=float)
# Equations of Motion 
# Evaluation of the right hand side of the ODE.
def eqmot(np, t, y, ydot, dmm, dll, dcc, GG, mu):
mu2= mu * mu
var1 = (( GG / sqrt(dll * dcc)) + (mu2 * GG / (dll*dcc)) ** 4) / (1( mu2 / (dll*dcc)**3))
var2 = (1/(dll*dcc))/(1(mu2 / (dll*dcc)**3))
var3 = (2*GG*mu/(dll*dcc))/(1(mu2 / (dll*dcc))**3)
var4 = (mu * math.sqrt(dll*dcc) / (dll*dcc)**2)/(1 (mu2 / (dll*dcc)**3))
ydot[0,0]=y[0,1]
ydot[1,0]=y[1,1]
ydot[0,1]= (var1)*y[0,1]var2*y[0,0]var3*y[1,1]+var4*y[1,0]
ydot[1,1]= var1*y[1,1]var2*y[1,0]+var3*y[0,1]+var4*y[0,0]
return ydot[0,0], ydot[1,0], ydot[0,1], ydot[1][1]
# Runge Kutta 
def runge(t, n, z, zt, dt2, savey, phi, m):
m = m + 1
if m == 1:
return 1, t, n, z, zt, dt2, savey, phi, m
elif m == 2:
for j in range(1,3):
for i in range(1, n + 1):
savey[i1][j1]=z[i1][j1]
phi[i1][j1]=zt[i1][j1]
z[i1][j1]=savey[i1][j1] + dt2 * zt[i1][j1]
t = t + dt2
return 1, t, n, z, zt, dt2, savey, phi, m
elif m == 3:
for j in range(1,3):
for i in range(1,n+1):
phi[i1][j1]=phi[i1][j1] + 2.0 * zt[i1][j1]
z[i1][j1]=savey[i1][j1] + dt2 * zt[i1][j1]
return 1, t, n, z, zt, dt2, savey, phi, m
elif m == 4:
for j in range(1,3):
for i in range(1,n+1):
phi[i1][j1] = phi[i1][j1] + 2.0 * zt[i1][j1]
z[i1][j1] = savey[i1][j1] + 2.0 * dt2 * zt[i1][j1]
t = t + dt2
return 1, t, n, z, zt, dt2, savey, phi, m
elif m == 5:
for j in range(1,3):
for i in range(1,n+1):
z[i1][j1]=savey[i1][j1] + (phi[i1][j1] + zt[i1][j1])*dt2/3.0
m = 0
return 1, t, n, z, zt, dt2, savey, phi, m
else:
print 'M has gone out of its range.'
# crfile 
def crfile(y,ydot,np,ntimes,t,tmax,intex,mu,dcc,dll,dmm):
energyL1 =0.5*dll*y[0][0]**2
energyL2 =0.5*dll*y[1][0]**2
energyC1 =0.5*(dcc*dll*y[0][1]dcc*dmm*y[1][1]**2/dcc)
energyC2 =0.5*(dcc*dll*y[1][1]dcc*dmm*y[0][1]**2/dcc)
energyM =dmm*y[0][0]*y[1][1]
energytot =energyL1+energyL2+energyC1+energyC2+energyM
return t, energyL1, energyL2, energyC1, energyC2, energytot
fout = open('output.txt', 'w')
# Integration 
t = 0.0
for i in range(1,ntimes+1):
for j in range(1, int(ifile+1)):
k = 1
while k == 1:
ydot[0,0], ydot[1,0], ydot[0,1], ydot[1][1] = eqmot(np, t, y , ydot, dmm, dll, dcc, GG, mu)
k, t, np, y, ydot, dt2, savey, phi, m = runge(t, np, y , ydot, dt2, savey, phi, m)
crfile(y,tdot,np,ntimes,t,tmax,i,mu,dcc,dll,dmm)
fout.write(crfile)
f.out.close()
Here is the original Fortran code. There are some errors with variable/array names in the 'crfile' subroutine, but I don't think there are any syntax mistakes.
Code:
Program electro_dimer_evol
!
! This program solves the coupled electronic dimer
! This is for the parallel coupled dimer
!
!  M.C. Zheng, T. Kottos Feb 2011
!
IMPLICIT NONE
INTEGER :: NP, i, ii, j, k, kk, ntimes, ifile, init
INTEGER :: runge, m, nsnaps
REAL(8) :: omega0, GG, mu, tmax
REAL(8) :: dll, dmm, dcc, dkk
REAL(8) :: dt, dt2, t
REAL(8), DIMENSION(:,:), ALLOCATABLE :: Y, YDOT, phi, savey
!SYSTEM PARAMETERS
!
! NP  number of LRC circuits units in the chain
! dll  inductance
! dmm  mutual inductance
! dkk  proportionality coefficient dmm = dkk*dll
! dcc  capacitance
! omega0  natural frequency of a unit: \omega0^2 = 1/LC
! GG  gain/loss: gg=R*SQRT[C/L] use new parameter
! mu  rescaled mutaul induction \mu = MC = dkk/omega0^2, use new para
!
! Tmax  maximum time evolution
! dt  integration time step
! ntimes  nsnaps of the number of times we are calling the output subroutine in order to write the output data
! init  initial excitation
!
Print*, 'Enter: GG'
Read(5,*) GG
Np = 2
init = 1
ntimes = 1000 !nn = Tmax/(dt*ntimes); nn INTEGER
dll = 1.d0
dcc = 1.d0
! dkk = .5d0
dmm = .8d0
! omega0 = 1.d0/dsqrt(dll*dcc)
mu = dmm/dll
! GG = 1.8d0
Tmax = 250.d0
dt = 0.01d0
!
ALLOCATE(Y(Np,2),YDot(Np,2),phi(Np,2),savey(Np,2))
OPEN(2, file = 'IE.dat')
m = 0
dt2 = dt/2.d0
ifile = tmax/(ntimes*dt)
!Initial Condition
CALL INCOND(Np, Y, init)
!LRC Integration
t=0.d0
Do i=1, ntimes
Do j=1, ifile
k = runge(t,NP,Y,YDOT,dt2,savey,phi,m)
Do while (k .EQ. 1)
call eqmot(NP,T,Y,YDOT,dmm,dll,dcc,GG,mu)
k = runge(t,NP,Y,YDOT,dt2,savey,phi,m)
END Do
END Do
call crfile(Y,YDOT,np,ntimes,t,tmax,i,mu,dcc,dll,dmm)
END Do
Close(2)
DEALLOCATE(Y,YDot,phi,savey)
END Program electro_dimer_evol
!
subroutine incond(np,Y,init)
integer np,init,i
real*8 Y(np,2)
do i=1,np
Y(i,1)=0.0
Y(i,2)=0.0
enddo
Y(init,1)=1.d0
return
end
! ======================================================================
SUBROUTINE crfile(Y,YDOT,np,ntimes,t,tmax,intex,mu,dcc,dll,dmm)
implicit none
integer j,NP,ntimes,intex
real*8 Y(NP,2),YDOT(NP,2),tmax
real*8 t,omega0,GG,mu,dcc,dll,dmm
real*8 energyL1, energyL2, energyC1, energyC2, energyM
real*8 energytot
energyL1 = .5d0*dll*Y(1,1)**2
energyL2 = .5d0*dll*Y(2,1)**2
energyC1 = .5d0*(dcc*dll*YDOT(1,1)dcc*dmm*YDOT(2,1))**2/dcc
energyC1 = .5d0*(dcc*dll*YDOT(2,1)dcc*dmm*YDOT(1,2))**2/dcc
energyM = dmm*Y(1,1)*Y(2,2)
energytot = energyL1+energyL2+energyC1+energyC2+energyM
!write out energies
write(2,*)t,energyL1,energyL2,energyC1,energyC2,energytot
return
end
! ======================================================================
SUBROUTINE eqmot(NP,T,Y,YDOT,dmm,dll,dcc,GG,mu)
!c***********************************************************************
!c Routine for evaluating the right hand side of the ODE.
!c The number of variables is 4*N.
!c The correspondences of the varibles are,
!c Y(NP,2)=(Q(np),I(np)); YDOT(NP,2)=(dQ(np)/dt,dI(np)/dt)
!c
!c***********************************************************************
implicit none
integer j,NP
real*8 Y(NP,2),YDOT(NP,2)
REAL*8 dll, dmm, dcc
real*8 t,GG,mu,pi,pi2
real*8 mu2
real*8 var1,var2,var3,var4
!c******************* EQUATIONS OF MOTION *******************************
mu2 = mu*mu
var1= ((GG/sqrt(dll*dcc))+(mu2*GG/(dll*dcc)**4))/(1(mu2/(dll*dcc)**3))
var2= (1/(dll*dcc))/(1(mu2/(dll*dcc)**3))
var3= (2*GG*mu/(dll*dcc))/(1(mu2/(dll*dcc)**3))
var4= (mu*sqrt(dll*dcc)/(dll*dcc)**2)/(1(mu2/(dll*dcc)**3))
YDOT(1,1) = Y(1,2)
YDOT(2,1) = Y(2,2)
YDOT(1,2) = var1*Y(1,2)var2*Y(1,1)var3*Y(2,2)+var4*Y(2,1)
YDOT(2,2) = var1*Y(2,2)var2*Y(2,1)+var3*Y(1,2)+var4*Y(1,1)
RETURN
END
!====================================================================
function runge(t,N,z,zt,dt2,savey,phi,m)
integer :: runge, m, i, j
integer :: N
real*8 dt2,t,phi(N,2),savey(N,2)
real(8), dimension(N,2) :: z,zt
m=m+1
selectcase(m)
case(1)
runge=1
return
case(2)
do j=1,2
do i=1,N
savey(i,j)=z(i,j)
phi(i,j)=zt(i,j)
z(i,j)=savey(i,j)+dt2*zt(i,j)
enddo
enddo
t=t+dt2
runge=1
return
case(3)
do j=1,2
do i=1,N
phi(i,j)=phi(i,j)+2D0*zt(i,j)
z(i,j)=savey(i,j)+dt2*zt(i,j)
enddo
enddo
runge=1
return
case(4)
do j=1,2
do i=1,N
phi(i,j)=phi(i,j)+2D0*zt(i,j)
z(i,j)=savey(i,j)+2D0*dt2*zt(i,j)
enddo
enddo
t=t+dt2
runge=1
return
case(5)
do j=1,2
do i=1,N
z(i,j)=savey(i,j)+(phi(i,j)+zt(i,j))*dt2/3D0
enddo
enddo
m=0
runge=0
return
end select
end function runge



03162011, 01:22 AM

#2

Senior Member
Registered: Sep 2006
Location: Liverpool  England
Distribution: slackware64 13.37 and current, Dragonfly BSD
Posts: 1,810

The usual way to approach these things is to insert print statements to output values from your loops so you can tell whats happening in your code. Python is very forgiving and just because the program may run doesn't mean you have no syntax errors! Ther errrors may be in a part of your code that is never being reached.For example here:
Code:
crfile(y,tdot,np,ntimes,t,tmax,i,mu,dcc,dll,dmm) # returned results not being stored anywhere
fout.write(crfile) # trying to write a function to a file
f.out.close() # possibly meant to be fout.close()
If the file write was being reached it would cause a TypeError. These are just initial obervations without going into the code in any detail. Start be putting in print statements in your loops (for and while) to output variables and create counters to see where the code is up to. This will show where it's getting stuck in loops. It looks like you are getting stuck in the inner while loop. Fix these then other errors will become apparent. Also it makes things more readable if you organize your code so you have for example any module imports at the top (as you already have), then function definitions then the main variable definitions then the actual program code getting input and looping when the functions are called. Your code mixes main code, variable definitions and function definitions making it very hard to read. Firstly though put some output so you can start to debug.



03162011, 06:54 AM

#3

LQ Newbie
Registered: Jul 2010
Posts: 17
Original Poster
Rep:

Thank you. I shall try these things right away.
Edit: At least one of the problems seems to be that k does not change.
Last edited by C.L.; 03162011 at 09:51 AM.



Thread Tools 
Search this Thread 


Posting Rules

You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
HTML code is Off



All times are GMT 5. The time now is 01:28 PM.

LinuxQuestions.org is looking for people interested in writing
Editorials, Articles, Reviews, and more. If you'd like to contribute
content, let us know.

Latest Threads
LQ News

