#
# This simple program approximates first- and second-order spatial derivatives using
# finite difference schemes.
# The function considered is exp(x) and we compute derivatives at x=1.0
# Order of convergence of the following schemes are studied:
# 1. Forward differece for first derivative
# 2. Backward differece for first derivative
# 3. 2nd order central differece for first derivative
# 4. 2nd order central differece for second derivative
#
# This is a part of teaching materials provided for the course
# "Computational Heat & Fluid Flow (ME 605)" taught at IIT Goa
# during the winter semester of 2018
#
# Written by: Y Sudhakar, IIT Goa
# email: sudhakar@iitgoa.ac.in
#
# Tested with python 3.6.3 on 15 August 2018
#
# Assignment for the students
# 1. second order one-sided difference for first derivative
# 2. Fourth order central difference approximation for first derivative
# 3. Fourth order central difference approximation for second derivative
#

import numpy as np
import math as m
import matplotlib.pyplot as plt

num_refine = 10;  # number of refinement level considered to plot order of convergence

h      = np.ones(num_refine)      # mesh spacing for each refinement
bd     = np.ones(num_refine)      # approximation by backward differencing
bderr  = np.ones(num_refine)      # error introduced in backward differencing
fd     = np.ones(num_refine)      # approximation by forward differencing
fderr  = np.ones(num_refine)      # error introduced in forward differencing
cd     = np.ones(num_refine)      # approximation by central differencing
cderr  = np.ones(num_refine)      # error introduced in central differencing
cd2    = np.ones(num_refine)      # approximation of second derivative by central differencing
cd2err = np.ones(num_refine)      # error introduced in central differencing for second derivative

# exact value at x=1.0
exac = m.exp(1.0)

# compute error introduced in the approximation by all the schemes
for i in range(len(h)):
    h[i]=0.1/(2**i)
    xp=1.0            # location and solution value at the considered node
    fp=m.exp(xp)
    xw=1.0-h[i]       # left side node (denoted with 'w' meaning 'west' as in finite volume context)
    fw=m.exp(xw)
    xe=1.0+h[i]       # right side node (denoted with 'e' meaning 'east' as in finite volume context)
    fe=m.exp(xe)
    bd[i]=(fp-fw)/h[i]         # backward difference approximation and error
    bderr[i]=abs(exac-bd[i])
    fd[i]=(fe-fp)/h[i]         # forward difference approximation and error
    fderr[i]=abs(exac-fd[i])
    cd[i]=0.5*(fe-fw)/h[i]     # central difference approximation and error
    cderr[i]=abs(exac-cd[i])
    cd2[i]=(fw+fe-2.0*fp)/h[i]/h[i]  # second order derivatives using central differencing
    cd2err[i]=abs(exac-cd2[i])

# print(exac)
# print(fd, bd, cd)

# Lines showing first and second order convergence
fir    = np.ones(num_refine)
sec    = np.ones(num_refine)
fou    = np.ones(num_refine)
fir[0] = 0.05
sec[0] = 0.05 #0.01
fou[0] = 0.05
i = 1
while i < len(fir):
    fir[i] = fir[i-1]*(h[i]/h[i-1])
    sec[i] = sec[i-1]*((h[i]/h[i-1])**2)
    fou[i] = fou[i-1]*((h[i]/h[i-1])**4)
    i = i+1

# plot all quantities
plt.plot(h,bderr,label='Backward')
plt.plot(h,fderr,'o',label='Forward')
plt.plot(h,cderr,label='Central')
plt.plot(h,cd2err,label='Central 2nd der.')
plt.plot(h,fir,'--',label='1st order')
plt.plot(h,sec,'--',label='2nd order')
plt.plot(h,fou,'o-',label='4th order')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$\Delta x$')
plt.ylabel('Truncation Error')
plt.legend()
plt.show()
#plt.savefig('truncation-error_FD.pdf')