source: anuga_work/production/sydney_2006/smf_volume_conservation.py @ 4030

"""
Determining volume conservation for SMF slump formulation.
Jane Sexton, GA 2006
"""

from math import sqrt, atan, degrees, exp, cosh, pi, radians
from Numeric import ones
from pylab import *
#import scipy
#from scipy.special import erf

ion()
hold(False)


x = arange(-20,25,0.1)
y = arange(-10,10,0.1)
X,Y = meshgrid(x,y)

# Grilli and Watts 2005 example
dx = 2.0
x0 = 10.0
y0 = 0.0
w = 2.0
kappa = 3
g = 9.8
d = 0.259
T = 0.052
theta = 15.0 #acos(radians(-d/T))  
lam0 = 5.0 #sqrt(g*d)
    
def func(x,y,kappad):
    
    numerator = exp(-((x-x0)/lam0)**2) - kappad*exp(-((x-dx-x0)/lam0)**2)
    denominator = cosh(kappa*(y-y0)/(w+lam0))**2
    return -numerator/denominator

eta1 = func(X,0,0.83)
eta2 = func(X,0,1)

#import sys; sys.exit()
#eta = func(X,Y,kappad)
#im1 = imshow(eta, cmap=cm.jet, interpolation='nearest')
#plot(X/lam0,eta)
#axis([0, 5, -2, 2])

figure(21)
plot(X/lam0,eta1,'b.-',X/lam0,eta2,'g.-')
title('Surface elevation of SMF: kappa = 0.83 and 1')
xlabel('x/lambda_0')
ylabel('eta(x,y=0)')
axis([0, 5, -2, 2])
savefig('check_fig2')

def set_integral_zero(x,dx):
    numerator = erf( (x-x0)/sqrt(lam0) )
    denominator = erf( (x-dx-x0)/sqrt(lam0) )
    return numerator/denominator

def set_integral_zero_2(x,x0):
    numerator = erf( (x-x0)/sqrt(lam0) )
    xg = d/tan(radians(theta)) + T/radians(sin(theta))
    #print 'xg', xg
    denominator = erf( (x-xg-2*x0)/sqrt(lam0) )
    return numerator/denominator

n = 100
kappad_vec1 = []
dx_vec1 = []
step = 0.001
dx = -step #+ 0.016
for i in range(n):
    dx += step
    kappad1 = float(max(set_integral_zero(x,dx)))
    dx_vec1.append(dx)
    kappad_vec1.append(kappad1)
    #print 'dx: %0.5f kappad: %f' % (dx, kappad1)

figure(22)
plot(dx_vec1,kappad_vec,'.-')
title('Volume conservation')
xlabel('delta x')
ylabel('kappa_d')
savefig('Volume_conservation')
#axis([0,0.1,0,80])

n = 60
kappad_vec2 = []
x0_vec2 = []
step = 0.1
x0 = -step
for i in range(n):
    x0 += step
    kappad2 = float(max(set_integral_zero_2(x,x0)))
    x0_vec2.append(x0)
    kappad_vec2.append(kappad2)
    print 'dx: %0.5f kappad: %f' % (x0, kappad2)


figure(23)
plot(x0_vec2,kappad_vec2,'.-')
title('Volume conservation - using xg')
xlabel('x0')
ylabel('kappa_d')
savefig('Volume_conservation_xg')