solver_burgers_model

import numpy as np
from numpy.linalg import norm

def newton_system(fnon, jacobian, x0, tolerance, max_iterations, *fnonargs):

    k = 0
    x = x0

    F = eval(fnon)(x,*fnonargs)
    n = F.size

    print(' k    f(xk)')

    while (norm(F,2) > tolerance and k <= max_iterations):
        J = eval(jacobian)(x,n,fnon,F,*fnonargs)

        delta = np.linalg.solve(J,-F)
        x = x + delta

        F = eval(fnon)(x,*fnonargs)

        print('{0:2.0f}  {1:2.2e}'.format(k, norm(F,2)))
        k += 1

    if (k >= max_iterations):
        print('Failed to converged')
    else:
        return x

import numpy as np
import copy

def calcJacobian(x,n,fnon,F0,*fnonargs):
  J = np.zeros((n,n), dtype=np.float64)
  h = 10e-8

  for k in range(0,n):
    xb = copy.copy(x)
    xb[k] = xb[k] + h

    F = eval(fnon)(xb,*fnonargs)

    for i in range(0,n):
        J[i,k] = (F[i] - F0[i]) / h

  return J

def step_function(x):
        if x <= 0.1:
            return 1
        else:
            return 0

# this is just to display the step function
from matplotlib import pyplot as plt
x_xis = np.linspace(0, 1, 400)
y = np.zeros(400)
for i in range(400):
    y[i] = step_function(x_xis[i])

plt.plot(xxis, y)
plt.title('Step Function')
plt.xlabel('Spatial coordinate x')
plt.ylabel('Solution u')
plt.grid(True)
plt.show()



def setupInitialData( m ):
    xL = 0
    xR = 1

    h = (xR - xL) / (m-1)                      # mesh size
    x = np.linspace(xL, xR, m).reshape((m, 1)) # grid points
    v = np.zeros(len(x))
    for i in range(len(x)):
        v[i] = step_function(x[i]) # initial data

    return v

def discreteBurgers(uk, ukp, dt, h, nu, ua, ub):
    # ua is boundary condition
    # ub is boundary condition
    # nu is kinematic viscosity
    m = uk.size
    # f to store values of the function for each point uk in space
    f = np.zeros((m-2, 1))

    # boundary conditions
    uL = ua
    uR = ub

    # left boundary (ERROR OCCURS HERE)
    f[0] = (uk[0] - ukp[1])/dt + uk[0] * (uk[0] - uL)/h - nu * (uk[1] - 2*uk[0] + uL)/h**2

    # Difference equations at each internal node
    for i in range(1, m-3):
        f[i] = (uk[i] - ukp[i+1])/dt + uk[i] * (uk[i] - uk[i-1])/h - nu * (uk[i+1] - 2*uk[i] + uk[i-1])/h**2

    # right boundary
    f[m-3] = (uk[m-3] - ukp[m-2])/dt + uk[m-3] * (uk[m-3] - uk[m-4])/h - nu * (uR - 2*uk[m-3] + uk[m-4])/h**2

    return f


def burgersModel( dtIn, nTimeSteps, mIn ):

  # Input:  dtIn       - time step size
  #         nTimeSteps - number of time steps
  #         mIn        - grid dimension
  # Output: sols       - array of solution vectors u^k

  dt = dtIn   # time step size
  m = mIn     # dimension of spatial mesh for [X1,X2]
  n = m-2     # number of nonlinear equations (excluding boundary condition)
  h = 1/(m-1) # grid size

  # Initial and boundary conditions
  uOld = initialData( m )
  u0 = np.zeros( (n,1) )

  # store solutions in time
  sols = []
  sols.append(uOld)

  for k in range( 0, nTimeSteps ):
    print('Time step: {0:1.0f}'.format(k))

    # Use previous time step as the initial guess for Newton
    u0 = uOld[1:n+1] # take the values in the middle, ignoring left and right boundary values

    # Call the Newton solver
    #uk, ukp, dt, h, nu, ua, ub
    nu = 0.01
    ua = 1 #u(0,t) = 1
    ub = 0 #u(1,t) = 0

    u = newton_system("discreteBurgers", "calcJacobian", u0, 1e-6, 10, uOld, dt, h, nu, ua, ub)

    # Update the previous time step
    uOld[1:n+1] = u[0:n] # insert new values as old now into the middle, not updating left and right boundaries

    # Store vector of solutions
    sols.append(copy.copy(uOld))

  return sols;

mPts = 101
tPts = 10
sols = burgersModel( 0.01, tPts, mPts )