LinuxQuestions.org
Share your knowledge at the LQ Wiki.
Go Back   LinuxQuestions.org > Forums > Non-*NIX Forums > Programming
User Name
Password
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



Reply
 
Search this Thread
Old 03-15-2011, 09:17 PM   #1
C.L.
LQ Newbie
 
Registered: Jul 2010
Posts: 17

Rep: Reputation: 1
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[i-1][j-1]=z[i-1][j-1]
                phi[i-1][j-1]=zt[i-1][j-1]
                z[i-1][j-1]=savey[i-1][j-1] + dt2 * zt[i-1][j-1]
        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[i-1][j-1]=phi[i-1][j-1] + 2.0 * zt[i-1][j-1]
                z[i-1][j-1]=savey[i-1][j-1] + dt2 * zt[i-1][j-1]
        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[i-1][j-1] = phi[i-1][j-1] + 2.0 * zt[i-1][j-1]
                z[i-1][j-1] = savey[i-1][j-1] + 2.0 * dt2 * zt[i-1][j-1]
        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[i-1][j-1]=savey[i-1][j-1] + (phi[i-1][j-1] + zt[i-1][j-1])*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
 
Old 03-16-2011, 02:22 AM   #2
bgeddy
Senior Member
 
Registered: Sep 2006
Location: Liverpool - England
Distribution: slackware64 13.37 and -current, Dragonfly BSD
Posts: 1,810

Rep: Reputation: 227Reputation: 227Reputation: 227
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.
 
Old 03-16-2011, 07:54 AM   #3
C.L.
LQ Newbie
 
Registered: Jul 2010
Posts: 17

Original Poster
Rep: Reputation: 1
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.; 03-16-2011 at 10:51 AM.
 
  


Reply

Tags
fortran, python


Thread Tools Search this Thread
Search this Thread:

Advanced Search

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Python output MTK358 Programming 4 06-14-2010 09:02 PM
Some questions about Python; speed, system integration Christian271 Programming 2 05-24-2010 05:32 PM
Colored output in Python MTK358 Programming 2 12-24-2009 09:01 PM
Python serial output seems strange edM Programming 2 04-27-2009 05:54 AM
tkinter integration to python 2.3.4 gaddargarson Programming 1 03-09-2005 05:58 AM


All times are GMT -5. The time now is 04:26 PM.

Main Menu
Advertisement
My LQ
Write for LQ
LinuxQuestions.org is looking for people interested in writing Editorials, Articles, Reviews, and more. If you'd like to contribute content, let us know.
Main Menu
Syndicate
RSS1  Latest Threads
RSS1  LQ News
Twitter: @linuxquestions
identi.ca: @linuxquestions
Facebook: linuxquestions Google+: linuxquestions
Open Source Consulting | Domain Registration