mac22_integrator.py

#!/usr/bin/python
from numpy import *
import sys

def forward22(uvec,vvec,xvec,tval,f,i):
	return f(uvec[i+1],vvec[i+1],xvec[i+1],tval)-f(uvec[i],vvec[i],xvec[i],tval)
def backward22(uvec,vvec,xvec,tval,f,i):
	return f(uvec[i],vvec[i],xvec[i],tval)-f(uvec[i-1],vvec[i-1],xvec[i-1],tval)
def mac22bf(uvec,vvec,dt,dx,fu,fv,gu,gv,xvec,t):
	# predict: Backward! (needs numerical BC treatment at j=N)
	n = len(uvec)
	(uhat,usol,vhat,vsol) = (zeros((n,1)),zeros((n,1)),zeros((n,1)),zeros((n,1)))
	#usol = zeros((n,1))
	#vhat = zeros((n,1))
	#vsol = zeros((n,1))
	for i in range(1,n):
		uhat[i]=uvec[i] + dt/dx*(backward22(uvec,vvec,xvec,t,fu,i)) + dt*gu(xvec[i],t)
		vhat[i]=vvec[i] + dt/dx*(backward22(uvec,vvec,xvec,t,fv,i)) + dt*gv(xvec[i],t)
	# Extrapolate the fluxes, sir!
        Fum1 = 2*fu(uvec[0],vvec[0],xvec[0],t)-fu(uvec[1],vvec[1],xvec[1],t);
        Fvm1 = 2*fv(uvec[0],vvec[0],xvec[0],t)-fv(uvec[1],vvec[1],xvec[1],t);
        usol[0] = uvec[0]+ dt/dx*(Fum1 - fu(uhat[0],vhat[0],xvec[0],t))
        vsol[0] = vvec[0]+ dt/dx*(Fvm1 - fv(uhat[0],vhat[0],xvec[0],t))

        Funp1 = 2*fu(uvec[n-1],vvec[n-1],xvec[n-1],t)-fu(uvec[n-2],vvec[n-2],xvec[n-2],t);
        Fvnp1 = 2*fv(uvec[n-1],vvec[n-1],xvec[n-1],t)-fv(uvec[n-2],vvec[n-2],xvec[n-2],t);
        uhat[n-1]=.5*(uvec[n-1]+uvec[n-1]+dt/dx*(Funp1-fu(uvec[n-2],vvec[n-2],xvec[n-2],t)));
        vhat[n-1]=.5*(vvec[n-1]+vvec[n-1]+dt/dx*(Fvnp1-fv(uvec[n-2],vvec[n-2],xvec[n-2],t)));
	for i in range(0,n-1):
		usol[i] = .5*(uvec[i]+uhat[i] + dt/dx*(forward22(uhat,vhat,xvec,t,fu,i)) + dt*gu(xvec[i],t))
		vsol[i] = .5*(vvec[i]+vhat[i] + dt/dx*(forward22(uhat,vhat,xvec,t,fv,i)) + dt*gv(xvec[i],t))
        # characteristic boundary conditions.
        (usol[0],vsol[0]) = (.5*(usol[0]+vsol[0]),.5*(usol[0]+vsol[0]))
        (usol[n-1],vsol[n-1])=(.5*(usol[n-1]+vsol[n-1]),-.5*(usol[n-1]+vsol[n-1]))

	return [usol, vsol]
def mac22fb(uvec,vvec,dt,dx,fu,fv,gu,gv,xvec,t):
	# predict: Foward! (needs numerical BC treatment at j=N)
	n=len(uvec)
	(uhat,usol,vhat,vsol) = (zeros((n,1)),zeros((n,1)),zeros((n,1)),zeros((n,1)))
	for i in range(0,n-1):
		uhat[i]=uvec[i] + dt/dx*(forward22(uvec,vvec,xvec,t,fu,i)) + dt*gu(xvec[i],t)
		vhat[i]=vvec[i] + dt/dx*(forward22(uvec,vvec,xvec,t,fv,i)) + dt*gv(xvec[i],t)
	# Extrapolate the fluxes, sir!
	Funp1 = 2*fu(uvec[n-1],vvec[n-1],xvec[n-1],t)-fu(uvec[n-2],vvec[n-2],xvec[n-2],t);
	Fvnp1 = 2*fv(uvec[n-1],vvec[n-1],xvec[n-1],t)-fv(uvec[n-2],vvec[n-2],xvec[n-2],t);
	uhat[n-1]=uvec[n-1]+dt/dx*(Funp1-fu(uvec[n-2],vvec[n-2],xvec[n-2],t));
	vhat[n-1]=vvec[n-1]+dt/dx*(Fvnp1-fv(uvec[n-2],vvec[n-2],xvec[n-2],t));

	Fum1 = 2*fu(uhat[0],vhat[0],xvec[0],t)-fu(uhat[1],vhat[1],xvec[1],t);
	Fvm1 = 2*fv(uhat[0],vhat[0],xvec[0],t)-fv(uhat[1],vhat[1],xvec[1],t);
	usol[0] = .5*(uvec[0]+uhat[0] + dt/dx*(Fum1 - fu(uhat[0],vhat[0],xvec[0],t)))
	vsol[0] = .5*(vvec[0]+vhat[0] + dt/dx*(Fvm1 - fv(uhat[0],vhat[0],xvec[0],t)))

	for i in range(1,n):
		usol[i] = .5*(uvec[i]+uhat[i] + dt/dx*(backward22(uhat,vhat,xvec,t,fu,i)) + dt*gu(xvec[i],t))
		vsol[i] = .5*(vvec[i]+vhat[i] + dt/dx*(backward22(uhat,vhat,xvec,t,fv,i)) + dt*gv(xvec[i],t))
	
	# characteristic boundary conditions.
	(usol[0],vsol[0]) = (.5*(usol[0]+vsol[0]),.5*(usol[0]+vsol[0]))
	(usol[n-1],vsol[n-1])=(.5*(usol[n-1]+vsol[n-1]),-.5*(usol[n-1]+vsol[n-1]))

	return [usol, vsol]
def mac22(FU,FV,GU,GV, tMax, N, cMax, fName_append):
	CFL=.95; #cMax=3; N = 1000

	x = linspace(-30*pi,30*pi,N);
	dx=x[1]-x[0]
	dt = CFL*dx/cMax
	t = linspace(0,tMax,tMax/dt);M = len(t)
	uvec=zeros((N,1))
	vvec=zeros((N,1))
	usol=zeros((M,N))
	vsol=zeros((M,N))

	Mstep = M/40;
	for n in range(0,M-1):
		if n%2==1:
			onestep = mac22fb(usol[n,:],vsol[n,:],dt,dx,FU,FV,GU,GV,x,t[n])
		else:	onestep = mac22bf(usol[n,:],vsol[n,:],dt,dx,FU,FV,GU,GV,x,t[n])
		usol[n+1,:] = onestep[0].T
		vsol[n+1,:] = onestep[1].T
		if n>Mstep:
			Mstep += M/40
			sys.stdout.write(".")
			sys.stdout.flush()
	print("done")
	savetxt('usol_%s.csv'%fName_append, usol, fmt='%.6f', delimiter=';')
	savetxt('vsol_%s.csv'%fName_append, vsol, fmt='%.6f', delimiter=';')
def conservlaw(uval,vval, x,t, mode,F):
	if mode == 1:
		c1=c2=1
		return F(c1,c2,uval,vval)
	elif mode==2:
		if abs(x) < 10*pi:
			c1=c2=1
			return F(c1,c2,uval,vval)
		else:
			c1=c2=3
			return F(c1,c2,uval,vval)
	elif mode==3:
		if abs(x) < 10*pi:
			c1=c2=1
			return F(c1,c2,uval,vval)
		else:
			c1=c2=1.001
			return F(c1,c2,uval,vval)
### Conti
sigma=2;
MODE=int(sys.argv[1])

FU = lambda uval,vval, x, t: conservlaw(uval,vval, x, t, MODE, lambda c1,c2,uval,vval: (c1-c2)/2*uval + (c1+c2)/2*vval)
GU = lambda x,t: exp(-sigma*x**2/2)*exp(-sigma**2*t**2/2)
GV = lambda x,t: 0
if MODE==1:
	FV = lambda uval,vval, x, t: conservlaw(uval,vval, x, t, MODE, lambda c1,c2,uval,vval: (c1+c2)/2*uval + (c1-c2)/2*vval)
	mac22(FU,FV,GU,GV, 100, 1000, 1, 'mode1')
if MODE==2:
	FV = lambda uval,vval, x, t: conservlaw(uval,vval, x, t, MODE, lambda c1,c2,uval,vval: uval + (c1-c2)/2*vval)
	mac22(FU,FV,GU,GV, 60, 3000, 3, 'mode2')
if MODE==3:
	FV = lambda uval,vval, x, t: conservlaw(uval,vval, x, t, MODE, lambda c1,c2,uval,vval: uval + (c1-c2)/2*vval)
	mac22(FU,FV,GU,GV, 60, 3000, 1.001, 'mode3')