/*******************************************************************************
2D advection example program which advects a Gaussian u(x,y) at a fixed velocity
Outputs: initial.dat – inital values of u(x,y)
final.dat – final values of u(x,y)
The output files have three columns: x, y, u
Compile with: gcc -o advection2D -std=c99 advection2D.c -lm
Notes: The time step is calculated using the CFL condition
********************************************************************************/
/*********************************************************************
Include header files
**********************************************************************/
#include
#include
/*********************************************************************
Main function
**********************************************************************/
int main(){
/* Grid properties */
const int NX=1000; // Number of x points
const int NY=1000; // Number of y points
const float xmin=0.0; // Minimum x value
const float xmax=1.0; // Maximum x value
const float ymin=0.0; // Minimum y value
const float ymax=1.0; // Maximum y value
/* Parameters for the Gaussian initial conditions */
const float x0=0.1; // Centre(x)
const float y0=0.1; // Centre(y)
const float sigmax=0.03; // Width(x)
const float sigmay=0.03; // Width(y)
const float sigmax2 = sigmax * sigmax; // Width(x) squared
const float sigmay2 = sigmay * sigmay; // Width(y) squared
/* Boundary conditions */
const float bval_left=0.0; // Left boudnary value
const float bval_right=0.0; // Right boundary value
const float bval_lower=0.0; // Lower boundary
const float bval_upper=0.0; // Upper bounary
/* Time stepping parameters */
const float CFL=0.9; // CFL number
const int nsteps=1500; // Number of time steps
/* Velocity */
const float velx=0.01; // Velocity in x direction
const float vely=0.01; // Velocity in y direction
/* Arrays to store variables. These have NX+2 elements
to allow boundary values to be stored at both ends */
float x[NX+2]; // x-axis values
float y[NX+2]; // y-axis values
float u[NX+2][NY+2]; // Array of u values
float dudt[NX+2][NY+2]; // Rate of change of u
float x2; // x squared (used to calculate iniital conditions)
float y2; // y squared (used to calculate iniital conditions)
/* Calculate distance between points */
float dx = (xmax-xmin) / ( (float) NX);
float dy = (ymax-ymin) / ( (float) NY);
/* Calculate time step using the CFL condition */
/* The fabs function gives the absolute value in case the velocity is -ve */
float dt = CFL / ( (fabs(velx) / dx) + (fabs(vely) / dy) );
/*** Report information about the calculation ***/
printf(“Grid spacing dx = %g\n”, dx);
printf(“Grid spacing dy = %g\n”, dy);
printf(“CFL number = %g\n”, CFL);
printf(“Time step = %g\n”, dt);
printf(“No. of time steps = %d\n”, nsteps);
printf(“End time = %g\n”, dt*(float) nsteps);
printf(“Distance advected x = %g\n”, velx*dt*(float) nsteps);
printf(“Distance advected y = %g\n”, vely*dt*(float) nsteps);
/*** Place x points in the middle of the cell ***/
/* LOOP 1 */
for (int i=0; i