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
|