source: trunk/anuga_work/development/2010-projects/anuga_1d/sww-sudi/run-discrepancy.py

import os, time, csv, pprint
from scipy.special import jn
from scipy import sin, cos, sqrt, linspace, pi, zeros
from rootsearch import *
from bisect_function import *
from Numeric import zeros,Float,dot
from gaussPivot import *
from config import g
from analytical_prescription import *


def j0(x):
    return jn(0.0, x)

def j1(x):
    return jn(1.0, x)

def j2(x):
    return jn(2.0, x)

def j3(x):
    return jn(3.0, x)

def jm1(x):
    return jn(-1.0, x)

def jm2(x):
    return jn(-2.0, x)

def bed(x):
    return x-1.0


def wu_at_O_by_formal_expansion(t):
    W1 = j0(4*pi/T)*cos(2*pi*t/T)
    U1 = -j1(4*pi/T)*sin(2*pi*t/T)

    gu = -2*pi/T*j0(4*pi/T)*sin(2*pi*t/T)
    gw = -2*pi/T*j1(4*pi/T)*cos(2*pi*t/T)
    fu = gw
    fw = gu + j1(4*pi/T)*sin(2*pi*t/T)
    W2 = gu*U1 + gw*W1 - 0.5*U1**2
    U2 = fu*U1 + fw*W1

    w = W1*A + W2*A*A
    u = U1*A + U2*A*A
    return w,u

def wu_at_O_by_fixed_point_iteration(t):
    w1 = A*j0(4.0*pi/T)*cos(2.0*pi*t/T)
    u1 = -1.0*A*j1(4.0*pi/T)*sin(2.0*pi*t/T)

    w2 = -0.5*u1**2.0 + A*j0(4.0*pi*(w1+1.0)**0.5/T)*cos(2.0*pi*(u1+t)/T)
    u2 = -1.0*A*j1(4.0*pi*(w1+1.0)**0.5/T)*(w1+1.0)**-0.5*sin(2.0*pi*(u1+t)/T)

    #w3 = -0.5*u2**2.0 + A*j0(4.0*pi*(w2+1.0)**0.5/T)*cos(2.0*pi*(u2+t)/T)
    #u3 = -1.0*A*j1(4.0*pi*(w2+1.0)**0.5/T)*(w2+1.0)**-0.5*sin(2.0*pi*(u2+t)/T)

    #w4 = -0.5*u3**2.0 + A*j0(4.0*pi*(w3+1.0)**0.5/T)*cos(2.0*pi*(u3+t)/T)
    #u4 = -1.0*A*j1(4.0*pi*(w3+1.0)**0.5/T)*(w3+1.0)**-0.5*sin(2.0*pi*(u3+t)/T)

    #w5 = -0.5*u4**2.0 + A*j0(4.0*pi*(w4+1.0)**0.5/T)*cos(2.0*pi*(u4+t)/T)
    #u5 = -1.0*A*j1(4.0*pi*(w4+1.0)**0.5/T)*(w4+1.0)**-0.5*sin(2.0*pi*(u4+t)/T)    
    return w2,u2


Time = linspace(0.0,T,1000)
N_T = len(Time)


Stage = zeros(N_T, Float)    
Veloc = zeros(N_T, Float)
for i in range(N_T):
    t=Time[i]
    #print "Here.................................t=",t
    zet, vel = prescribe(0.0,t)
    Stage[i] = zet
    Veloc[i] = vel


Stage_sudi = zeros(N_T, Float)   
Veloc_sudi = zeros(N_T, Float)
for i in range(N_T):
    t=Time[i]
    Stage_sudi[i],Veloc_sudi[i] = wu_at_O_by_formal_expansion(t)#wu_at_O_by_fixed_point_iteration(t)


num=len(Stage)
error_w=(1.0/num)*sum(abs(Stage-Stage_sudi))*h_0
error_u=(1.0/num)*sum(abs(Veloc-Veloc_sudi))*sqrt(g*h_0)
print "error_w=", error_w
print "error_u=", error_u

#######Numerical C-G at centroids########
table = zeros((num,3),Float)
for r in range(num):
    for c in range(3):
        if c==0:
            table[r][c] = Time[r]
        elif c==1:
            table[r][c] = Stage[r]
        else:
            table[r][c] = Veloc[r]

outname = "%s%f%s" %("Exact_at_O_", Tp, ".csv")                
outfile = open(outname, 'w')
writer = csv.writer(outfile)
for row in table:
    writer.writerow(row)        
outfile.close()


#######Numerical stage and velocity########
table = zeros((num,3),Float)
for r in range(num):
    for c in range(3):
        if c==0:
            table[r][c] = Time[r]
        elif c==1:
            table[r][c] = Stage_sudi[r]
        else:
            table[r][c] = Veloc_sudi[r]

outname = "%s%f%s" %("Numeric_w2u2_FE_at_O_", Tp, ".csv")             
outfile = open(outname, 'w')
writer = csv.writer(outfile)
for row in table:
    writer.writerow(row)        
outfile.close()

from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot

hold(False)
clf()
plot1 = subplot(211)
plot(Time/T,Stage*h_0,'b-', Time/T,Stage_sudi*h_0,'k--')
#xlabel('t/T')
ylabel('Stage')
#plot1.set_xlim([0.000,0.030])
#plot1.set_ylim([0.980,1.005]) #([-9.0e-3,9.0e-3])
legend(('C-G', 'M-R'),
       'lower left', shadow=False)

plot2 = subplot(212)
plot(Time/T,Veloc*sqrt(g*h_0),'b-', Time/T,Veloc_sudi*sqrt(g*h_0),'k--')
xlabel('t/T')
ylabel('Velocity')
#plot1.set_xlim([0.0,1.1])
#plot2.set_ylim([-0.05,0.05]) #([-1.0e-12,1.0e-12])
legend(('C-G', 'M-R'),
       'upper right', shadow=False)


#filename = "discrepancy-closer"
#filename += str(i)
#filename += ".eps"
#savefig(filename)
#show()

#plot(Time,Vel_at_O)
show()