In [2]:

    

import numpy as np
from numpy import pi
from scipy.integrate import ode

mrand=lambda: 2*np.random.rand()-1 # shortcut random number between -1 and 1
def roll(arr,tup):
    '''a shortcut. roll(arr,(S1,S2)) applies roll with S1 steps on 1st axis and then with S2 on the 2nd one'''
    return np.roll(np.roll(arr,tup[0],axis=0),tup[1],axis=1)

def generate_v_field(n=17,max_periods=4):
    '''Generates a random divergence-less velocity field on the periodic domain [0,2*pi]x[0,2*pi].
       n - number of sinusoidal terms
       max_periods - maximal periods each term has over the domain 2pi (large value = high frequency)
    '''
    x,y=np.meshgrid(np.linspace(0,2*pi,N+1)[:-1],
                    np.linspace(0,2*pi,N+1)[:-1])
    u=0*x;
    v=0*x;
    for i in xrange(n):
        A=np.random.randn()
        nx=np.round(max_periods*mrand())
        ny=np.round(max_periods*mrand())
        phi=pi*mrand()
        u+=A*ny*np.sin(nx*x+ny*y+phi)
        v-=A*nx*np.sin(nx*x+ny*y+phi)
    return u,v

def velocity_dot_grad_c(u,v,C):
    '''calculates the advective derivative of C with the velocity field (u,v)/
       That is, u dC/dx+v dC/dy'''
    
    return  (
            (0.25/dx)*(
             (C+roll(C,[0,-1]))*(roll(u,[1,-1])+roll(u,[0,-1]))-
             (C+roll(C,[0, 1]))*(roll(u,[0, 0])+roll(u,[1, 0])))+
            (0.25/dx)*(
             (C+roll(C,[-1, 0]))*(roll(v,[0,0])+roll(v,[0,-1]))-
             (C+roll(C,[ 1, 0]))*(roll(v,[1,-1])+roll(v,[1,0])))
            )

def laplacian(C):
    '''calculates the laplacian of C'''
    return (
            roll(C,[-1,0])+roll(C,[0,-1])+roll(C,[1,0])+roll(C,[0,1])-4*C
            )/(dx**2)

def dvdt_maker(u,v,eta,sources):
    '''returns a function that calculates the time derivative but is usable in the integrator that needs a vector
       input rather than a matrix.'''
    return lambda t,C: Cdot(C.reshape(N,N),u,v,eta,sources).flatten()

def Cdot(C,u,v,eta,sources):
    '''Returns the time derivative of the concentration field C, under diffusion and advection
       C        - concentration field
       (u,v)    - components of velocity field
       eta      - diffusion constant
       sources  - sources (and sinks)'''
    return -velocity_dot_grad_c(u,v,C)+eta*laplacian(C)+sources

def point_sources():
    z=np.zeros([N,N])
    z[np.round(N/4),np.round(N/4)]=1
    z[np.round(3*N/4),np.round(3*N/4)]=-1
    return z

def gauss_sources(width=0.2,randomize_positions=False):
    x,y=np.meshgrid(np.linspace(0,2*pi,N+1)[:-1],
                np.linspace(0,2*pi,N+1)[:-1])
    q=np.exp(-((x-pi)/width)**2-((y-pi)/width)**2)

    if randomize_positions:
        return roll(q,[int(N*mrand()) for i in [1,2]])-roll(q,[int(N*mrand()) for i in [1,2]])
    else:
        return roll(q,[N/4,N/4])-roll(q,[-N/4,-N/4])

def calculate_steady_plume(eta,
                           number_of_points=200,
                           tolerance=1.0e-8,
                           max_iterations=100,
                           verbose=False,
                           width=0.2,
                           randomize_positions=False,
                           timestep=1,
                           seed='auto'):
    '''
        Calculates a steady plume in a divergenceless field. 
        Returns (c,u,v,sources,flag, errors):
            c       - the steady concentration field
            (u, v)  - the velocity field
            sources - the sources and sinks in the equation
            flag    - True if the solution is steady to within the tolerance
            errors  - the time evolution of the error (=norm of time derivative). 
                      when error<tolerance the solver stops and exits
        Arguments:
            eta                 - diffusion constant
            number_of_points    - in each dimension (result is an NxN matrix)
            width               - width of the gaussian sources
            randomize_positions - whether to choose the positions of the source and sink randomly
                                  if False, they are at the middle of the top-right and bottom-left quartiles
            verbose             - whether to print progress in each iteration
            timestep            - check if solution has converged every timestep (in simulation time)
            seed                - random seed (must be int, otherwise no explicit seed specified to numpy)
    '''
    global dx, N
    if type(seed)==type(3):
        np.random.seed(seed)
        
    N=number_of_points
    dx=2*pi/N
    u,v=generate_v_field()
    sources=gauss_sources(randomize_positions=randomize_positions,width=width)
    fdot=dvdt_maker(u,v,eta,sources)
    
    solver=ode(fdot).set_integrator('dopri5')
    solver.set_initial_value(0*u.flatten())
    i=0;
    #prev=np.zeros(N*N);
    error=1
    errs=[]
    while solver.successful() and i<max_iterations and error > tolerance:
        solver.integrate(solver.t+timestep)
        error=(dx**2)*np.sum(np.abs(fdot(0,solver.y)))
        #prev=solver.y
        errs.append(error)
        i+=1
        if verbose:
            print 'step %d, error=%e' %(i, error)
    return (solver.y.reshape(N,N), u, v,  sources, error<tolerance ,errs)


    

In [3]:

    

res=[calculate_steady_plume(0.1,number_of_points=q,verbose=True,tolerance=1e-6,seed=413, timestep=0.1,max_iterations=1e100) for q in [50, 100,200]]


    

    




step 1, error=2.755898e-01
step 2, error=2.823388e-01
step 3, error=2.651124e-01
step 4, error=2.554579e-01
step 5, error=2.402448e-01
step 6, error=2.044713e-01
step 7, error=1.699573e-01
step 8, error=1.407794e-01
step 9, error=1.157042e-01
step 10, error=9.600307e-02
step 11, error=8.362041e-02
step 12, error=7.367565e-02
step 13, error=6.475715e-02
step 14, error=5.722738e-02
step 15, error=5.155340e-02
step 16, error=4.704025e-02
step 17, error=4.247996e-02
step 18, error=3.809899e-02
step 19, error=3.424559e-02
step 20, error=3.105847e-02
step 21, error=2.826977e-02
step 22, error=2.575632e-02
step 23, error=2.357520e-02
step 24, error=2.176462e-02
step 25, error=2.018942e-02
step 26, error=1.876101e-02
step 27, error=1.748843e-02
step 28, error=1.636939e-02
step 29, error=1.538129e-02
step 30, error=1.447624e-02
step 31, error=1.361676e-02
step 32, error=1.278678e-02
step 33, error=1.198062e-02
step 34, error=1.121274e-02
step 35, error=1.049053e-02
step 36, error=9.818859e-03
step 37, error=9.200208e-03
step 38, error=8.628983e-03
step 39, error=8.097953e-03
step 40, error=7.599403e-03
step 41, error=7.129522e-03
step 42, error=6.686048e-03
step 43, error=6.267890e-03
step 44, error=5.874546e-03
step 45, error=5.506933e-03
step 46, error=5.164000e-03
step 47, error=4.843318e-03
step 48, error=4.542520e-03
step 49, error=4.259410e-03
step 50, error=3.992438e-03
step 51, error=3.740875e-03
step 52, error=3.503741e-03
step 53, error=3.280925e-03
step 54, error=3.071735e-03
step 55, error=2.875589e-03
step 56, error=2.691808e-03
step 57, error=2.519243e-03
step 58, error=2.357475e-03
step 59, error=2.205225e-03
step 60, error=2.062700e-03
step 61, error=1.928631e-03
step 62, error=1.803350e-03
step 63, error=1.685517e-03
step 64, error=1.575597e-03
step 65, error=1.472237e-03
step 66, error=1.375833e-03
step 67, error=1.285319e-03
step 68, error=1.200637e-03
step 69, error=1.121402e-03
step 70, error=1.047136e-03
step 71, error=9.777217e-04
step 72, error=9.127184e-04
step 73, error=8.518654e-04
step 74, error=7.950498e-04
step 75, error=7.417732e-04
step 76, error=6.921154e-04
step 77, error=6.455887e-04
step 78, error=6.022256e-04
step 79, error=5.615048e-04
step 80, error=5.237305e-04
step 81, error=4.881728e-04
step 82, error=4.551906e-04
step 83, error=4.242117e-04
step 84, error=3.954156e-04
step 85, error=3.684522e-04
step 86, error=3.433191e-04
step 87, error=3.199028e-04
step 88, error=2.980014e-04
step 89, error=2.776336e-04
step 90, error=2.585684e-04
step 91, error=2.408606e-04
step 92, error=2.242810e-04
step 93, error=2.088366e-04
step 94, error=1.944271e-04
step 95, error=1.810569e-04
step 96, error=1.686172e-04
step 97, error=1.570506e-04
step 98, error=1.463796e-04
step 99, error=1.364395e-04
step 100, error=1.272732e-04
step 101, error=1.186671e-04
step 102, error=1.107070e-04
step 103, error=1.032202e-04
step 104, error=9.633032e-05
step 105, error=8.995462e-05
step 106, error=8.406771e-05
step 107, error=7.865656e-05
step 108, error=7.362936e-05
step 109, error=6.893687e-05
step 110, error=6.452524e-05
step 111, error=6.034179e-05
step 112, error=5.644146e-05
step 113, error=5.272488e-05
step 114, error=4.927381e-05
step 115, error=4.609112e-05
step 116, error=4.313326e-05
step 117, error=4.042074e-05
step 118, error=3.785931e-05
step 119, error=3.545539e-05
step 120, error=3.318545e-05
step 121, error=3.101311e-05
step 122, error=2.898479e-05
step 123, error=2.703784e-05
step 124, error=2.523600e-05
step 125, error=2.358823e-05
step 126, error=2.214966e-05
step 127, error=2.078745e-05
step 128, error=1.952367e-05
step 129, error=1.829481e-05
step 130, error=1.714593e-05
step 131, error=1.601845e-05
step 132, error=1.495148e-05
step 133, error=1.392368e-05
step 134, error=1.296256e-05
step 135, error=1.216909e-05
step 136, error=1.149706e-05
step 137, error=1.086258e-05
step 138, error=1.023665e-05
step 139, error=9.643044e-06
step 140, error=9.046805e-06
step 141, error=8.465749e-06
step 142, error=7.887432e-06
step 143, error=7.320393e-06
step 144, error=6.808872e-06
step 145, error=6.481904e-06
step 146, error=6.197631e-06
step 147, error=5.917874e-06
step 148, error=5.631982e-06
step 149, error=5.338171e-06
step 150, error=5.030379e-06
step 151, error=4.712415e-06
step 152, error=4.379626e-06
step 153, error=4.041860e-06
step 154, error=3.789268e-06
step 155, error=3.700487e-06
step 156, error=3.620071e-06
step 157, error=3.514897e-06
step 158, error=3.395676e-06
step 159, error=3.251320e-06
step 160, error=3.086250e-06
step 161, error=2.899727e-06
step 162, error=2.690048e-06
step 163, error=2.469572e-06
step 164, error=2.394652e-06
step 165, error=2.425836e-06
step 166, error=2.435251e-06
step 167, error=2.424705e-06
step 168, error=2.381819e-06
step 169, error=2.316485e-06
step 170, error=2.217378e-06
step 171, error=2.093625e-06
step 172, error=1.941946e-06
step 173, error=1.771730e-06
step 174, error=1.821987e-06
step 175, error=1.895211e-06
step 176, error=1.942626e-06
step 177, error=1.958964e-06
step 178, error=1.944503e-06
step 179, error=1.895887e-06
step 180, error=1.815546e-06
step 181, error=1.703971e-06
step 182, error=1.563632e-06
step 183, error=1.441134e-06
step 184, error=1.548175e-06
step 185, error=1.630517e-06
step 186, error=1.684310e-06
step 187, error=1.705315e-06
step 188, error=1.693423e-06
step 189, error=1.647364e-06
step 190, error=1.569350e-06
step 191, error=1.459322e-06
step 192, error=1.321439e-06
step 193, error=1.280904e-06
step 194, error=1.392538e-06
step 195, error=1.477527e-06
step 196, error=1.530817e-06
step 197, error=1.551675e-06
step 198, error=1.538498e-06
step 199, error=1.491132e-06
step 200, error=1.411026e-06
step 201, error=1.299193e-06
step 202, error=1.159746e-06
step 203, error=1.201745e-06
step 204, error=1.313629e-06
step 205, error=1.395435e-06
step 206, error=1.445661e-06
step 207, error=1.461475e-06
step 208, error=1.443561e-06
step 209, error=1.391118e-06
step 210, error=1.306290e-06
step 211, error=1.190628e-06
step 212, error=1.048019e-06
step 213, error=1.168693e-06
step 214, error=1.276078e-06
step 215, error=1.353148e-06
step 216, error=1.397267e-06
step 217, error=1.407494e-06
step 218, error=1.383005e-06
step 219, error=1.324485e-06
step 220, error=1.233230e-06
step 221, error=1.111556e-06
step 222, error=1.033955e-06
step 223, error=1.164242e-06
step 224, error=1.265683e-06
step 225, error=1.335604e-06
step 226, error=1.371969e-06
step 227, error=1.373832e-06
step 228, error=1.341206e-06
step 229, error=1.274835e-06
step 230, error=1.176490e-06
step 231, error=1.048674e-06
step 232, error=1.048125e-06
step 233, error=1.171922e-06
step 234, error=1.266370e-06
step 235, error=1.328656e-06
step 236, error=1.357230e-06
step 237, error=1.351119e-06
step 238, error=1.310451e-06
step 239, error=1.236158e-06
step 240, error=1.130119e-06
step 241, error=9.951516e-07
step 1, error=2.513782e-01
step 2, error=2.513618e-01
step 3, error=2.513168e-01
step 4, error=2.500511e-01
step 5, error=2.348582e-01
step 6, error=2.041876e-01
step 7, error=1.701843e-01
step 8, error=1.418545e-01
step 9, error=1.191488e-01
step 10, error=1.018163e-01
step 11, error=8.931329e-02
step 12, error=7.901231e-02
step 13, error=7.010522e-02
step 14, error=6.293125e-02
step 15, error=5.718772e-02
step 16, error=5.204897e-02
step 17, error=4.711721e-02
step 18, error=4.268209e-02
step 19, error=3.894967e-02
step 20, error=3.584437e-02
step 21, error=3.301888e-02
step 22, error=3.042110e-02
step 23, error=2.808850e-02
step 24, error=2.603133e-02
step 25, error=2.422775e-02
step 26, error=2.261281e-02
step 27, error=2.114310e-02
step 28, error=1.979066e-02
step 29, error=1.853722e-02
step 30, error=1.737056e-02
step 31, error=1.627844e-02
step 32, error=1.525376e-02
step 33, error=1.429307e-02
step 34, error=1.339247e-02
step 35, error=1.255342e-02
step 36, error=1.177534e-02
step 37, error=1.105518e-02
step 38, error=1.038691e-02
step 39, error=9.764896e-03
step 40, error=9.183791e-03
step 41, error=8.640472e-03
step 42, error=8.132825e-03
step 43, error=7.657751e-03
step 44, error=7.210709e-03
step 45, error=6.790570e-03
step 46, error=6.395670e-03
step 47, error=6.024262e-03
step 48, error=5.674620e-03
step 49, error=5.345190e-03
step 50, error=5.034762e-03
step 51, error=4.742384e-03
step 52, error=4.467018e-03
step 53, error=4.207792e-03
step 54, error=3.963910e-03
step 55, error=3.734494e-03
step 56, error=3.518704e-03
step 57, error=3.315628e-03
step 58, error=3.124470e-03
step 59, error=2.944460e-03
step 60, error=2.774909e-03
step 61, error=2.615276e-03
step 62, error=2.464985e-03
step 63, error=2.323551e-03
step 64, error=2.190444e-03
step 65, error=2.065189e-03
step 66, error=1.947291e-03
step 67, error=1.836317e-03
step 68, error=1.731873e-03
step 69, error=1.633566e-03
step 70, error=1.541098e-03
step 71, error=1.454202e-03
step 72, error=1.372668e-03
step 73, error=1.296263e-03
step 74, error=1.224589e-03
step 75, error=1.157317e-03
step 76, error=1.094125e-03
step 77, error=1.034510e-03
step 78, error=9.782370e-04
step 79, error=9.251215e-04
step 80, error=8.749942e-04
step 81, error=8.276757e-04
step 82, error=7.830034e-04
step 83, error=7.408078e-04
step 84, error=7.009781e-04
step 85, error=6.633563e-04
step 86, error=6.278115e-04
step 87, error=5.942251e-04
step 88, error=5.624995e-04
step 89, error=5.325307e-04
step 90, error=5.042251e-04
step 91, error=4.774862e-04
step 92, error=4.522320e-04
step 93, error=4.283994e-04
step 94, error=4.059103e-04
step 95, error=3.846829e-04
step 96, error=3.646573e-04
step 97, error=3.457656e-04
step 98, error=3.279422e-04
step 99, error=3.111163e-04
step 100, error=2.952233e-04
step 101, error=2.802019e-04
step 102, error=2.659965e-04
step 103, error=2.525511e-04
step 104, error=2.398057e-04
step 105, error=2.277237e-04
step 106, error=2.162748e-04
step 107, error=2.054121e-04
step 108, error=1.951019e-04
step 109, error=1.853176e-04
step 110, error=1.760324e-04
step 111, error=1.672199e-04
step 112, error=1.588557e-04
step 113, error=1.509179e-04
step 114, error=1.433849e-04
step 115, error=1.362332e-04
step 116, error=1.294411e-04
step 117, error=1.229913e-04
step 118, error=1.168659e-04
step 119, error=1.110494e-04
step 120, error=1.055259e-04
step 121, error=1.002803e-04
step 122, error=9.529941e-05
step 123, error=9.056804e-05
step 124, error=8.607421e-05
step 125, error=8.180621e-05
step 126, error=7.775317e-05
step 127, error=7.390429e-05
step 128, error=7.024750e-05
step 129, error=6.677374e-05
step 130, error=6.347573e-05
step 131, error=6.034564e-05
step 132, error=5.737497e-05
step 133, error=5.455687e-05
step 134, error=5.188576e-05
step 135, error=4.935447e-05
step 136, error=4.695476e-05
step 137, error=4.467777e-05
step 138, error=4.251136e-05
step 139, error=4.044991e-05
step 140, error=3.848830e-05
step 141, error=3.662166e-05
step 142, error=3.484558e-05
step 143, error=3.315549e-05
step 144, error=3.155858e-05
step 145, error=3.001761e-05
step 146, error=2.856188e-05
step 147, error=2.718597e-05
step 148, error=2.585953e-05
step 149, error=2.461863e-05
step 150, error=2.341641e-05
step 151, error=2.228507e-05
step 152, error=2.120999e-05
step 153, error=2.018916e-05
step 154, error=1.923077e-05
step 155, error=1.829959e-05
step 156, error=1.747595e-05
step 157, error=1.658336e-05
step 158, error=1.578846e-05
step 159, error=1.503480e-05
step 160, error=1.431238e-05
step 161, error=1.365593e-05
step 162, error=1.297018e-05
step 163, error=1.242529e-05
step 164, error=1.175087e-05
step 165, error=1.120185e-05
step 166, error=1.067065e-05
step 167, error=1.014361e-05
step 168, error=9.677255e-06
step 169, error=9.196094e-06
step 170, error=8.824488e-06
step 171, error=8.355267e-06
step 172, error=7.978906e-06
step 173, error=7.568039e-06
step 174, error=7.280419e-06
step 175, error=6.926541e-06
step 176, error=6.538012e-06
step 177, error=6.275363e-06
step 178, error=5.932687e-06
step 179, error=5.670806e-06
step 180, error=5.382111e-06
step 181, error=5.206888e-06
step 182, error=4.954194e-06
step 183, error=4.645765e-06
step 184, error=4.472294e-06
step 185, error=4.233829e-06
step 186, error=4.182307e-06
step 187, error=3.905006e-06
step 188, error=3.643930e-06
step 189, error=3.499254e-06
step 190, error=3.312355e-06
step 191, error=3.157194e-06
step 192, error=3.037160e-06
step 193, error=2.885574e-06
step 194, error=2.774791e-06
step 195, error=2.777700e-06
step 196, error=2.494910e-06
step 197, error=2.403374e-06
step 198, error=2.302083e-06
step 199, error=2.158674e-06
step 200, error=2.100868e-06
step 201, error=1.965352e-06
step 202, error=1.960265e-06
step 203, error=1.779940e-06
step 204, error=1.860854e-06
step 205, error=1.673185e-06
step 206, error=1.667337e-06
step 207, error=1.469604e-06
step 208, error=1.539991e-06
step 209, error=1.301020e-06
step 210, error=1.279765e-06
step 211, error=1.217249e-06
step 212, error=1.198343e-06
step 213, error=1.127737e-06
step 214, error=1.050783e-06
step 215, error=1.033929e-06
step 216, error=9.903024e-07
step 1, error=2.513274e-01
step 2, error=2.513274e-01
step 3, error=2.513177e-01
step 4, error=2.497634e-01
step 5, error=2.332904e-01
step 6, error=2.018945e-01
step 7, error=1.679283e-01
step 8, error=1.401612e-01
step 9, error=1.180239e-01
step 10, error=1.012942e-01
step 11, error=8.907789e-02
step 12, error=7.884913e-02
step 13, error=6.999625e-02
step 14, error=6.298784e-02
step 15, error=5.735468e-02
step 16, error=5.216143e-02
step 17, error=4.724304e-02
step 18, error=4.286933e-02
step 19, error=3.925860e-02
step 20, error=3.617236e-02
step 21, error=3.335187e-02
step 22, error=3.076858e-02
step 23, error=2.845707e-02
step 24, error=2.642708e-02
step 25, error=2.464070e-02
step 26, error=2.303097e-02
step 27, error=2.155358e-02
step 28, error=2.018785e-02
step 29, error=1.892198e-02
step 30, error=1.774438e-02
step 31, error=1.664446e-02
step 32, error=1.560914e-02
step 33, error=1.463480e-02
step 34, error=1.372222e-02
step 35, error=1.287319e-02
step 36, error=1.208651e-02
step 37, error=1.135802e-02
step 38, error=1.068182e-02
step 39, error=1.005206e-02
step 40, error=9.464097e-03
step 41, error=8.914264e-03
step 42, error=8.396107e-03
step 43, error=7.907732e-03
step 44, error=7.448215e-03
step 45, error=7.016272e-03
step 46, error=6.610227e-03
step 47, error=6.228233e-03
step 48, error=5.868524e-03
step 49, error=5.529583e-03
step 50, error=5.210102e-03
step 51, error=4.909098e-03
step 52, error=4.625658e-03
step 53, error=4.358859e-03
step 54, error=4.107763e-03
step 55, error=3.871472e-03
step 56, error=3.649071e-03
step 57, error=3.439678e-03
step 58, error=3.242479e-03
step 59, error=3.056738e-03
step 60, error=2.881789e-03
step 61, error=2.717022e-03
step 62, error=2.561882e-03
step 63, error=2.415847e-03
step 64, error=2.278359e-03
step 65, error=2.148908e-03
step 66, error=2.027037e-03
step 67, error=1.912304e-03
step 68, error=1.804300e-03
step 69, error=1.702699e-03
step 70, error=1.607205e-03
step 71, error=1.517547e-03
step 72, error=1.433453e-03
step 73, error=1.354512e-03
step 74, error=1.280376e-03
step 75, error=1.210710e-03
step 76, error=1.144986e-03
step 77, error=1.082941e-03
step 78, error=1.024362e-03
step 79, error=9.690435e-04
step 80, error=9.168017e-04
step 81, error=8.674727e-04
step 82, error=8.208865e-04
step 83, error=7.768815e-04
step 84, error=7.353115e-04
step 85, error=6.960477e-04
step 86, error=6.589581e-04
step 87, error=6.239284e-04
step 88, error=5.908498e-04
step 89, error=5.596198e-04
step 90, error=5.301364e-04
step 91, error=5.022978e-04
step 92, error=4.760146e-04
step 93, error=4.511973e-04
step 94, error=4.277579e-04
step 95, error=4.056061e-04
step 96, error=3.846595e-04
step 97, error=3.648450e-04
step 98, error=3.461022e-04
step 99, error=3.283706e-04
step 100, error=3.115948e-04
step 101, error=2.957225e-04
step 102, error=2.807066e-04
step 103, error=2.664992e-04
step 104, error=2.530455e-04
step 105, error=2.402983e-04
step 106, error=2.282135e-04
step 107, error=2.167518e-04
step 108, error=2.058780e-04
step 109, error=1.955596e-04
step 110, error=1.857678e-04
step 111, error=1.764735e-04
step 112, error=1.676508e-04
step 113, error=1.592762e-04
step 114, error=1.513250e-04
step 115, error=1.437759e-04
step 116, error=1.366077e-04
step 117, error=1.298016e-04
step 118, error=1.233383e-04
step 119, error=1.172004e-04
step 120, error=1.113713e-04
step 121, error=1.058341e-04
step 122, error=1.005761e-04
step 123, error=9.558246e-05
step 124, error=9.083943e-05
step 125, error=8.633498e-05
step 126, error=8.205681e-05
step 127, error=7.799373e-05
step 128, error=7.413434e-05
step 129, error=7.046968e-05
step 130, error=6.698984e-05
step 131, error=6.368579e-05
step 132, error=6.054973e-05
step 133, error=5.757476e-05
step 134, error=5.475431e-05
step 135, error=5.208167e-05
step 136, error=4.954831e-05
step 137, error=4.714527e-05
step 138, error=4.485944e-05
step 139, error=4.268450e-05
step 140, error=4.061508e-05
step 141, error=3.864552e-05
step 142, error=3.677125e-05
step 143, error=3.498798e-05
step 144, error=3.329118e-05
step 145, error=3.167692e-05
step 146, error=3.014107e-05
step 147, error=2.867987e-05
step 148, error=2.729027e-05
step 149, error=2.596876e-05
step 150, error=2.471233e-05
step 151, error=2.351811e-05
step 152, error=2.238319e-05
step 153, error=2.130502e-05
step 154, error=2.027970e-05
step 155, error=1.930510e-05
step 156, error=1.837832e-05
step 157, error=1.749673e-05
step 158, error=1.665782e-05
step 159, error=1.585980e-05
step 160, error=1.510057e-05
step 161, error=1.437844e-05
step 162, error=1.369164e-05
step 163, error=1.303826e-05
step 164, error=1.241667e-05
step 165, error=1.182541e-05
step 166, error=1.126264e-05
step 167, error=1.072667e-05
step 168, error=1.021630e-05
step 169, error=9.730209e-06
step 170, error=9.267303e-06
step 171, error=8.826401e-06
step 172, error=8.406389e-06
step 173, error=8.006282e-06
step 174, error=7.625126e-06
step 175, error=7.262035e-06
step 176, error=6.916188e-06
step 177, error=6.586744e-06
step 178, error=6.272940e-06
step 179, error=5.974079e-06
step 180, error=5.689448e-06
step 181, error=5.418388e-06
step 182, error=5.160272e-06
step 183, error=4.914505e-06
step 184, error=4.680520e-06
step 185, error=4.457752e-06
step 186, error=4.245665e-06
step 187, error=4.043747e-06
step 188, error=3.851500e-06
step 189, error=3.668482e-06
step 190, error=3.494410e-06
step 191, error=3.328721e-06
step 192, error=3.171009e-06
step 193, error=3.020911e-06
step 194, error=2.878062e-06
step 195, error=2.742119e-06
step 196, error=2.612750e-06
step 197, error=2.489652e-06
step 198, error=2.372540e-06
step 199, error=2.261136e-06
step 200, error=2.155175e-06
step 201, error=2.054403e-06
step 202, error=1.958596e-06
step 203, error=1.867518e-06
step 204, error=1.780967e-06
step 205, error=1.698832e-06
step 206, error=1.621056e-06
step 207, error=1.547106e-06
step 208, error=1.477061e-06
step 209, error=1.410643e-06
step 210, error=1.347745e-06
step 211, error=1.288126e-06
step 212, error=1.231595e-06
step 213, error=1.178030e-06
step 214, error=1.127330e-06
step 215, error=1.079416e-06
step 216, error=1.034180e-06
step 217, error=9.914823e-07

In [11]:

    

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
x,y=np.meshgrid(np.linspace(0,2*pi,N+1)[:-1],
                    np.linspace(0,2*pi,N+1)[:-1])
fig , axs= plt.subplots(nrows=2,ncols=2,figsize=(14,14))
axs[0,0].streamplot(N*x/(2*pi),N*y/(2*pi),res[-1][1],res[-1][2],color='b',density=4)
axs[0,0].set_title('Velocity')
axs[0,1].imshow(res[-1][0],cmap=plt.get_cmap('jet'))
axs[0,1].set_title('Concentration')
axs[1,0].imshow(res[-1][3],cmap=plt.get_cmap('jet'))
axs[1,0].set_title('Sources')
axs[1,1].imshow(res[-1][0],cmap=plt.get_cmap('jet'))
axs[1,1].streamplot(N*x/(2*pi),N*y/(2*pi),res[-1][1],res[-1][2],color='w',density=2)
axs[1,1].set_title('Velocity + conc.')
for ax in axs.flatten():
    ax.set_xlim([0,N])
    ax.set_ylim([0,N])
    ax.set_xticks([])
    ax.set_yticks([])
plt.show()


    

In [15]:

    

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
fig, ax= plt.subplots(figsize=(10,10))
ax.semilogy(res[-1][-1],'g--',label='n=100')
ax.semilogy(res[0][-1],'r',label='n=50')
ax.semilogy(res[-2][-1],'b',label='n=200')
plt.legend()
ax.set_xlabel('iterations')
ax.set_ylabel('error')
plt.show()