Skip to content

Instantly share code, notes, and snippets.

@traviscj

traviscj/mac22_integrator.py

Created December 8, 2009 04:23
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save traviscj/251398 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save traviscj/251398 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. traviscj created this gist Dec 8, 2009.
    123 changes: 123 additions & 0 deletions mac22_integrator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,123 @@
    #!/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')