Untitled

mail@pastecode.io avatar
unknown
python
2 years ago
9.4 kB
1
Indexable
Never
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator
from numpy.core.function_base import linspace
import scipy.sparse as ssp
from scipy.sparse.linalg import spsolve
from scipy.linalg import circulant

def matrix(N):
    dx = 1/(N+1)
    matrix = (1/dx**2)*(np.diag(np.ones(N-1), -1) + np.diag(-2*np.ones(N)) + np.diag(np.ones(N-1), 1))
    return matrix

def eulerstep(Tdx, uold, dt):
    unew = uold + dt*np.dot(Tdx,uold)
    return unew

def trapztep(Tdx, uold, dt):
    unew = np.linalg.solve(np.eye(len(Tdx)) - np.multiply(Tdx, dt/2), uold + dt/2*np.dot(Tdx,uold))
    return unew

def LW_Matrix(N, amu):
    # row = np.zeros(shape=(1,N))
    # row[0][0] = 1- amu**2
    # row[0][1] = (amu/2)*(amu+1)
    # row[0][-1] = (amu/2)*(amu-1)
    # return ssp.csc_matrix(circulant(row)).todense()
    Tdx = (np.diag((amu/2)*(1+amu)*np.ones(N-1), -1) + np.diag((1-np.power(amu,2))*np.ones(N)) + np.diag((amu/2)*(amu-1)*np.ones(N-1), 1))
    Tdx[0][-1] = (amu/2)*(1+amu)
    Tdx[-1][0] = (amu/2)*(amu-1)
    return Tdx

def BTdx(d, dx, sparse ): # without coefficients (dt/(2*dx)) as in own derivation.
    N = int(1/dx)
    TdxB = (d/dx)*( np.diag(np.ones(N-1),-1) + -2*np.diag(np.ones(N)) + np.diag(np.ones(N-1),1) 
                    )
    TdxB[0,-1] = TdxB[1,0]
    TdxB[-1,0] = TdxB[0,1]
    return TdxB if sparse == False else ssp.csr_matrix(TdxB)

def BurgerMatrix(dx, dt, TdxB, sparse):
    N = int(1/dx)
    sub = dt - 1
    main = (2*dx/dt) - (dt/dx)
    sup = dt + 1 
    Tdx_LW = (dt/(2*dx))*(     sub* np.diag(np.ones(N-1),-1)
                             + main*np.diag(np.ones(N))
                             + sup*np.diag(np.ones(N-1),1)
                                                                 )
    Tdx_LW += (dt/(2*dx))*TdxB
    Tdx_LW[0,-1] = Tdx_LW[1,0] # VÄND PÅ MINUSTECKEN OM FEL
    Tdx_LW[-1,0] = Tdx_LW[0,1]
    return Tdx_LW if sparse == False else ssp.csr_matrix(Tdx_LW)

def BurgerStep(uold, TdxB, burger, d, dx, dt, N):
    TdxB[0,-1] = TdxB[1,0]
    TdxB[-1,0] = TdxB[0,1]
    if(ssp.issparse(TdxB) and ssp.issparse(burger)):
        unew = spsolve(    (ssp.eye(N) -(d*dt/(2*dx))*TdxB ) , np.dot(burger, ssp.csr_matrix(uold) )       ).toarray()
    else:
        unew = np.linalg.solve( np.eye(N) - np.multiply(TdxB, (d*dt)/(2*dx)) , np.dot(burger, uold)    )
    return unew


def LaxWen(u, M):
    unew = np.dot(M, u)
    return unew

def advectionsSolver(N , M, g, t0, tf):
    dx = 1/N
    dt = (tf-t0)/M
    # a = dx/dt
    a = 0.9*dx/dt
    amu = a*dt/dx 

    uold = [g(i/N, i/M) for i in range(1, N+1)]
    m = []
    m.append(np.array([0, *uold, 0]))
    xx = np.linspace(0,1,N+2)
    tt = np.linspace(t0,tf,M+1)
    T,X = np.meshgrid(tt, xx)
    lax = LW_Matrix(N, amu = amu)
    for _ in range(1, len(tt)):
        uold = LaxWen(uold, lax)
        m.append(np.array([0, *uold, 0]))
    m = np.array(m)
    return m, T, X

def L2(m, tend, taskParam):
    errV = []
    N = len(m)
    t = np.linspace(0,tend,N)
    for i in range(0,N):
        errV.append(np.sqrt(sum(k*k/(N) for k in m[i])))
    plt.plot(t, errV)
    plt.xlabel('$dt$')
    plt.ylabel('$err$')
    plt.title('Error for Custom g(x), with $a \mu = $ ' + taskParam)
    plt.show()

def step(choice, Tdx, uold, dt):
    if choice == 0:
        return eulerstep(Tdx, uold, dt)
    if choice == 1:
        return trapztep(Tdx, uold, dt)

def plot(m, T, X, task, **param):
    fig = plt.figure()
    ax = plt.axes(projection = '3d')
    wireframe = ax.plot_surface(T,X, np.transpose(m), cmap = cm.plasma) #nipy_spectral, coolwarm , inferno, ocean, jet, gnu, viridis

    if task == 1:
        title = 'Solution to Diffusion Equation Custom $u(0,x) = g(x) = $'
        for i in param:
            title = title + param[i] + ' , '
        ax.set_title(title)
        ax.set_xlabel('T')
        ax.set_ylabel('X')
        ax.set_zlabel('u(t,x)')
    if task == 2:
        title = 'Solution to Advection Equation Custom $u(0,x) = g(x,t)$,' +'$a \mu = $'
        for i in param:
            title = title + param[i] + ' , '
        ax.set_title(title)
        ax.set_xlabel('T')
        ax.set_ylabel('X')
        ax.set_zlabel('u(t,x)')
    if task == 3:
        title = 'Convection Diffusion Custom $u(0,x) = g(x,t)$, '
        for i in param:
            title = title + param[i] + ' , '
        ax.set_title( title )
        ax.set_xlabel('T')
        ax.set_ylabel('X')
        ax.set_zlabel('u(t,x)')
    if task == 4:
        title = 'Viscous Burger´s Equation, '
        for i in param:
            title = title + param[i] + ' , '
        ax.set_title( title )
        ax.set_xlabel('T')
        ax.set_ylabel('X')
        ax.set_zlabel('u(t,x)')
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter('{x:.02f}')
    fig.colorbar(wireframe, shrink = 0.5, aspect = 5)
    plt.show()

def task1(N , M, g, t0, tf, choice):
    # Next line should compute a 2-dimensional vector u(m∆t,n∆x)
    # u(t,x) has u(t,0) = u(t,1) = 0
    # u(0,x) = g(x)
    dt = (tf-t0)/M
    xx = np.linspace(0,1,N+2)
    tt = np.linspace(t0,tf,M+1)
    T,X = np.meshgrid(tt, xx) 
    uold = [g(i/(N), i/(M)) for i in range(1, N+1)]
    m = []
    m.append(np.array([0, *uold, 0]))
    for _ in range(1, np.size(tt)):
        uold = step(choice, Tdx = matrix(N), uold = uold, dt = dt )
        m.append(np.array([0, *uold, 0]))
    m = np.array(m)
    return m, T, X
    
def convdif(a,d,dt,dx): # u' = (d*Tdx - a*Sdx)u
    N = int(1/dx)
    M = int(1/dt)
    g = lambda x,t:  x**3 -(4/3)*x**2 +(1/3)*x 
    uold = [g(i/N, i/M) for i in range(1, N+1)]
    xx = np.linspace(0,1,N)
    tt = np.linspace(0,1,M+1)
    T,X = np.meshgrid(tt, xx) 
    Tdx = (np.diag(1*np.ones(N-1), -1) + np.diag(-2*np.ones(N)) + np.diag(1*np.ones(N-1), 1))/dx**2
    Sdx = (np.diag(-1*np.ones(N-1), -1) + np.diag(0*np.ones(N)) + np.diag(np.ones(N-1), 1))/(2*dx)
    Master = (d*Tdx - a*Sdx)
    Master[0][-1] = d/(dx**2)+a/(2*dx)
    Master[-1][0] = d/(dx**2)-a/(2*dx)
    m = []
    m.append(np.array([*uold]))
    for _ in range(0, np.size(tt)-1):
        uold = step(1, Master, uold, dt)
        m.append(np.array([*uold]))
    m = np.array(m)
    print(np.shape(np.transpose(m)))
    return m, T, X

def viscBurgerSolver(g, d):
    N = 300
    delx = 1/N
    M = 5*N
    delt = 1/M
    uold = [g(i/N, i/M) for i in range(1, N+1)]
    BTedX = BTdx(d,delx, sparse = False )
    Burger = BurgerMatrix(delx, delt, BTedX, sparse = False)
    m = []
    m.append(np.array([*uold]))
    for _ in range(0, M):
        uold = BurgerStep(uold,BTedX, Burger, d, delx, delt, N)
        m.append(np.array([*uold]))
    m = np.array(m)

    xx = np.linspace(0,1,N)
    tt = np.linspace(0,1,M+1)
    T,X = np.meshgrid(tt, xx)
    print(np.shape(m))
    return m, T , X


if __name__ == '__main__':
    c = 1e-6
    N = 200
    tstart = 0
    tend= 1
    t = tend - tstart
    M_CFL =  2*t*(N+1)**2 # CFL-condition dt/dx**2 <= 1/2 \iff t/M * (N+1)^2
    chc = 1
    if chc == 0:
        Mnbr = int(np.ceil(M_CFL))
    if chc == 1:
        Mnbr = 10*N
    # m, T, X = task1(N = 100, M = Mnbr, g = lambda x,t: (x-0.5)**2, t0= tstart, tf=tend, choice = chc )
    # m, T, X = task1(N = N, M = Mnbr-10, g = lambda x,t: c*(x)**3 + c/((t)**3), t0= tstart, tf=tend, choice = chc) #BREAK IT!
    # m, T, X = task1(N = N, M = Mnbr, g = lambda x,t: np.exp(x /t), t0= tstart, tf = tend, choice = chc )
    # m, T, X = task1(N = N, M = Mnbr-16, g = lambda x,t: np.exp(x /t), t0= tstart, tf = tend, choice = chc) #BREAK IT AND GET PANTS!
    # plot(m = m, T = T, X = X, task=1)
    # advectionsSolver(N = N, M = Mnbr, g = lambda x: np.exp(-100*(x-0.5)**2))
    # m, T, X = advectionsSolver(N = 100, M = Mnbr, g = lambda x,t: np.exp(-100*np.power((x-0.5),2)), t0 = 0, tf = 1 )
    # plot(m = m, T = T, X = X)
    #a = -10
    #d = a/(N**2)
    # d = 1000000
    # print('Pe') # d > a * dx/2 = a ^(1/N^2), 
    # print(np.abs(d/a))
    # m, T, X = convdif(a, d, 1/Mnbr, 1/(N+1) )
    # plot(m, T, X, task = 3, a_and_d = '$(a,d) = (1, 1.11e-5)$', g_of_x = '$g(x) = x^3 -x^2 4/3 + x/3 $')
    # L2(m,5, taskParam='$0.9$')



    ##########UPPGIFT 4.1 / 4.1###########################
    # Jag tror mina matriser är fel, har tolkat LW(uold) som en matrismultiplikation U*uold
    # där U härleds till (dt/2dx)*tridiag([dt-1, (2*dx/dt) - (dt/dx) ,dt+1 ])
    # + de hörnens punkter så det blir en circulent matris. Jag har hittat att U:s invers är exakt U^T fast alla nollskilda element byts ut till 1/a_ij. Inte till någon hjälp precis då det är (I - TdxCirc)^-1 som ska lösas i BurgerStep. 

    # Iförsig är den matrisen cirkulant med
    
    #                   Circ( 1/(1 + 2/dx**2), -dx**2, 0, ... , 0, -dx**2 )!!!!!
    
    # Vi testar imorgon godnatt!!

    # Sparse fungerar ej... Men vi behöver omimplementera matriserna och kanske lösa den bättre

    m, T, X = viscBurgerSolver(lambda x,t: np.exp(-100*np.power((x-0.5),2) ),  d = 0.01) #Lägre d fungerar men inte 0.1..
    plot(m, T, X, 4, dcoeff = '$d = 0.01$', g_of_x = '$g(x) = \exp(-100 (x-0.5)^2)$' )