LinuxQuestions.org

LinuxQuestions.org (/questions/)
-   Programming (http://www.linuxquestions.org/questions/programming-9/)
-   -   Python -numerical integration, not getting any output. (http://www.linuxquestions.org/questions/programming-9/python-numerical-integration-not-getting-any-output-868802/)

C.L. 03-15-2011 08:17 PM

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


bgeddy 03-16-2011 01:22 AM

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.

C.L. 03-16-2011 06:54 AM

Thank you. I shall try these things right away.

Edit: At least one of the problems seems to be that k does not change.


All times are GMT -5. The time now is 09:21 AM.