/****************** 2D Advection Diffusion Equation with Constant Coeficients****************************/

/*********************READ ME*************************************/
/*
* To run the code:
*
* 1.Open the terminal in the folder containing Code.c
* 2.Type in : gcc Code.c -lm
* 3.Type in : ./a.out
*
* Now .VTK files will be created in this folder
*
* To plot the Results use Visit
* And Plot Theta Variable
*
* ****************************************************************/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<stdbool.h>

#define PI 3.14156

// Problem Parameters
#define NX (100+2)    // No of Cells in X direction including ghost cells
#define NY NX       // No of Cells in Y direction including ghost cells
#define N NX*NY    // Total Number of Cells including ghost cells
#define L 1.0E0    // Length of Square domain
#define K 0.0001E0 // Diffusion Coefficient
#define u 1.0E0    // Since u = 1.0 we need not explicitly use it
#define v 1.0E0    // Since v = 1.0 we need not explicitly use it

#define dx L/(NX - 2)
#define dy L/(NY - 2)
#define dt 1.0E-5

void Set_Up_Boundary_Conditions(double * Theta_n);
void Solve_Matrix_System_BiCGSTAB(double * Theta_np1, double * b);
void Get_Ax_vector(double * x,double * Ax_vector);
void Write_Output_VTK(int SiNo,double * Theta_n);

int main()
{
double T_f = .1;       // Final time

double deltax = dx;
double deltay = dy;
double deltat = dt;

//Checking Stability
//For the reason Read Finite Volume Methods for Hyperbolic Problems Leveque Page 449
{
if( ( (deltat/deltax) + (deltat/deltay) ) < 1.0E0 )
{
   // The current dt is stable
}
else
{
   printf("The current dt used is too large!!! Exit\n");
   return 0;
}
}

int END = (int)(T_f/dt) + 1;
double time = 0.0E0;

int i,j,k,n;                                  // Loop counters

int l;

// Initialisation of Theta varialbles
double * Theta_n = malloc(sizeof(double)*N); // Theta at time n
double * Theta_np1 = malloc(sizeof(double)*N); // Theta at time n + 1

double * b       = calloc(N,sizeof(double)); // Setting b vector

// Set up Initial Condition
for(j=1;j<(NY-1);j++)
{
for(i=1;i<(NX-1);i++)
{
   l = j*NX + i;
   // Given Initial Condition is sin^2(4*pi*x*y)
   Theta_n[l] = sin(4.0E0*PI*( (i-1)*dx + dx/2.0E0 )*( (j-1)*dy+ dy/2.0E0 ) )*sin(4.0E0*PI*( (i-1)*dx + dx/2.0E0 )*( (j-1)*dy + dy/2.0E0 ) );
}
}

// Boundary Conditions are taken care of
Set_Up_Boundary_Conditions(Theta_n);

Write_Output_VTK(0,Theta_n);

// Starting the main loop
for(n=1;n<END;n++)
{

// Boundary Conditions are taken care of
Set_Up_Boundary_Conditions(Theta_n);

// Solving the advection equation using Donor Cell Upwind
for(i=1;i<(NX-1);i++)
{
   for(j=1;j<(NY-1);j++)
   {
   l = j*NX+i;
   Theta_n[l] = Theta_n[l] + (dt/dx)*(Theta_n[l-1] - Theta_n[l]) + (dt/dy)*(Theta_n[l-NX] - Theta_n[l]);
   }
}

// Solving the Diffusion part using Crank-Nicholson Scheme

// Boundary Conditions are taken care of
Set_Up_Boundary_Conditions(Theta_n);

// 1.Set up b vector
for(j=1;j<(NY-1);j++)
{
   for(i=1;i<(NX-1);i++)
   {

1. Use the Taylor expansions for f(x0 − 2h), f(x0 − h), f(x0 + h),...