<jupyter_start><jupyter_text>V-Shaped Multigrid approach for a 1D diffusion problem Abhijeet Sutar<jupyter_code>import numpy as np
from matplotlib import pyplot as plt
from numba import jit
from IPython.display import clear_output
@jit(nopython = True)
def ABgen(nx, dx, k, Ar, q, TA, TB):
    A = np.zeros((nx,nx))
    B = np.zeros(nx)
    
    ae = k * Ar / dx
    aw = 0 - ae
    ae2 = 2 * ae
    Sp = 0 - ae2
    Su = q * Ar * dx
    ##for A
    for i in range(nx):
        if (i == 0):
            A[i,i] = 3 * ae
            A[i,i+1] = aw
        elif (i>0) and (i < (nx-1)):
            A[i,i-1] = aw
            A[i,i] = ae2
            A[i,i+1] = aw
        elif (i == (nx-1)) :
            A[i,i-1] = aw
            A[i,i] = 3 * ae
        
    ##for B
    for i in range(nx):
        if (i == 0):
            B[i] = Su + ae2 * TA
        elif (i > 0) and (i < (nx-1)) :
            B[i] = Su
        elif(i == (nx-1)):
            B[i] = Su + ae2 * TB
            
    return A, B
@jit(nopython = True)
def TDMA(nx, dx, k, Ar, q, TA, TB):
    A , B = ABgen(nx, dx, k, Ar, q, TA, TB)
    n = len(B)
    x = np.zeros(n)
    
    for i in range(1,n):
        q = A[i, i-1] / A[i-1,i-1]
        #A[i,i-1] = A[i, i-1] / A[i-1,i-1]
        A[i,i] = A[i,i] - (A[i,i-1] * q)
                   
        B[i] = B[i] - q * B[i-1]
        
    q = B[-1] / A[-1,-1]
    x[-1] = q
    for i in range(n-2, -1, -1):
        q = ( B[i] - A[i,i+1] * q ) / A[i,i]
        x[i] = q
    return x
@jit(nopython = True)
def GaussSeidel(A, B, y):
    ## y is initial guess
    matsize = len(y)
    for i in range(matsize):
        y[i] = 0
        for j in range(matsize):
            if (j!=i) and (A[i,j]!=0):
                y[i] += -A[i,j] / A[i,i] * y[j]
        
        y[i] += B[i] / A[i,i]
        
    return y
    
            
def StdGS1grid(nx, dx, k, Ar, q, TA, TB):
    Ah , Bh = ABgen(nx, dx, k, Ar, q, TA, TB)
    erms = 1
    yh = np.zeros_like(Bh)
    yh[:] = min(TA,TB)
    stepcount = 0
    while(erms > 1e-6):
        yh = GaussSeidel(Ah,Bh,yh)
        rh = residual(Ah,Bh,yh)
        erms = rmsf(rh)
        stepcount +=1
        clear_output(wait=True)
        print(stepcount, erms, flush=True)
    print(stepcount)
    return yh, erms
@jit(nopython = True)
def rmsf(R):
    r = np.sqrt(np.mean(np.square(R)))
    return r
def residual(A, B, y):
    r = np.zeros_like(y)
    r = B - A @ y
    return r
@jit(nopython= True)
def Rextrapolate2(y):
    ## y needs even number of elements
    y2h = np.zeros(int(len(y)/2))
    for i in range(len(y2h)):
        y2h[i] = (y[2*i] + y[2*i+1]) / 2
    
    return y2h
@jit(nopython = True)
def Interpolateback2(e):
    ## e is a vector
    xaxis = len(e)
    yaxis = 2 * len(e) 
    I = np.zeros((yaxis, xaxis))
    j_arr = np.arange(len(e)) + 1
    for j in j_arr:
        if (j != 1):
            I[2*j-3,j-1] = 0.25
        I[2*j-2,j-1] = 0.75
        I[2*j-1,j-1] = 0.75
        if (j != j_arr[-1]):
            I[2*j  ,j-1] = 0.25
        
    return I
def V2multgrid(nx, dx, k, Ar, q, TA, TB):
    if ((nx % 4) != 0):
        print("use nx value divisible by 4")
        return 
    Ah , Bh = ABgen(nx, dx, k, Ar, q, TA, TB)
    erms = 1
    yh = np.zeros_like(Bh)
    yh[:] = min(TA,TB)
    stepcount = 0
    while(erms > 1e-6):
        ## fine grid
        for i in range(5):
            yh = GaussSeidel(Ah,Bh,yh)
            stepcount +=1
        rh = residual(Ah,Bh,yh)
        r2h = Rextrapolate2(rh)
        
        ## half resolution
        A2h, B2h = ABgen(int(nx/2), 2*dx, k, Ar, q, TA, TB)
        e2h = np.zeros_like(r2h)
        for i in range(10):
            e2h = GaussSeidel(A2h,r2h,e2h)
            stepcount +=1
        r2h = residual(A2h,r2h,e2h)
        r4h = Rextrapolate2(r2h)
        
        ## quarter resolution
        A4h, B4h = ABgen(int(nx/4), 4*dx, k, Ar, q, TA, TB)
        e4h = np.zeros_like(r4h)
        for i in range(10):
            e4h = GaussSeidel(A4h,r4h,e4h)
            stepcount +=1
        r4h = residual(A4h,r4h,e4h)
        
        ## Interpolate
        ep2h = Interpolateback2(e4h)@e4h
        ec2h = e2h + ep2h
        
        ##smoothening
        for i in range(2):
            ec2h = GaussSeidel(A2h,r2h,ec2h)
            stepcount +=1
        
        ##final Interpolation
        eph = Interpolateback2(ec2h)@ec2h
        yc = yh + eph
        
        rfh = residual(Ah,Bh,yc)
        erms = rmsf(rfh)
        
        clear_output(wait=True)
        print(stepcount, erms, flush=True)
    
    print(stepcount)
    return yc, erms 
# RR = 01(B17ME001)
%timeit ymultG, ermsmultG = V2multgrid(20, 0.05, 5, 0.01, 20e3, 101, 501)
%timeit y1G, erms1G = StdGS1grid(20, 0.05, 5, 0.01, 20e3, 101, 501)
%timeit x = TDMA(20, 0.05, 5, 0.01, 20e3, 101, 501)
ymultG, ermsmultG = V2multgrid(20, 0.05, 5, 0.01, 20e3, 101, 501)
y1G, erms1G = StdGS1grid(20, 0.05, 5, 0.01, 20e3, 101, 501)
x = TDMA(20, 0.05, 5, 0.01, 20e3, 101, 501)
format_row = "{:<25}" * (4)
print(format_row.format("TDMA solution", "Multi Grid", "Error Multi grid", "Error Percentage MG"))
for i in range(len(x)):
    print(format_row.format(x[i], ymultG[i], x[i] - ymultG[i], (x[i] - ymultG[i])/x[i]*100))
format_row = "{:<25}" * (4)
print(format_row.format("TDMA solution", "Single Grid GS", "Error Single grid", "Error Percentage SG"))
for i in range(len(x)):
    print(format_row.format(x[i], y1G[i], x[i] - y1G[i], (x[i] - y1G[i])/x[i]*100))
fig1 = plt.figure(figsize = (8,11), dpi = 100)
plt.plot(x, marker = "x", label = "TDMA", linewidth = 0.5)
plt.plot(ymultG, marker = "*", label="V-cycle Multi Grid", linewidth = 0.5)
plt.grid()
plt.legend()
plt.xlabel("nodes")
plt.ylabel("Temperature in deg C")

y = np.zeros(5)
A = np.array([[375,-125,0,0,0], [-125,250,-125,0,0],[0,-125,250,-125,0],[0,0,-125,250,-125],[0,0,0,-125,375]])
B = np.array([29000,4000,4000,4000,54000])
A, B
n = len(B)
x = np.zeros(n)

for i in range(1,n):
    q = A[i, i-1] / A[i-1,i-1]
    #A[i,i-1] = A[i, i-1] / A[i-1,i-1]
    A[i,i] = A[i,i] - (A[i,i-1] * q)

    B[i] = B[i] - q * B[i-1]

q = B[-1] / A[-1,-1]
x[-1] = q
for i in range(n-2, -1, -1):
    q = ( B[i] - A[i,i+1] * q ) / A[i,i]
    x[i] = q
    
x<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18522103904485318, 'fraction_numerical': 0.06340755082284608, 'mean_word_length': 1.8832558139534883, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '544928', 'license': 'MIT', 'n_tokens_mistral': 3076, 'n_tokens_neox': 2843, 'n_words': 660, 'path': 'V_Multigrid_1D_Diffusion.ipynb', 'repo_name': 'Abhi-Sutar/BasicCFDpy', 'size': '6198'}