/****************** 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 1.0 // Diffusion Coefficient

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

void Set_Up_Boundary_Conditions(double * Theta_n);
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
{
if( dt < 0.5*dx*dx )
{
   // 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 x,y;

// Set up Initial Condition
for(j=1;j<(NY-1);j++)
{
for(i=1;i<(NX-1);i++)
{
   l = j*NX + i;

   // A Circle at .5,.5 with r = .25
   x = (i-1)*deltax + deltax/2.0;
   y = (j-1)*deltay + deltay/2.0;

   if ( ( (x-.5)*(x-.5) + (y-.5)*(y-.5) - .25*.25 ) <= 0.0 ) // Inside
   Theta_n[l] = 1.0;
   else
   Theta_n[l] = 0.0;

}

STOKES DRAG (Drag force acting on a sphere suspended in a fluid)...