/*
Women in Mathematics Summer School
27th May - 1st June 2013
ICTP - Trieste, Italy
Th. Katsaounis, University of Crete, 28/5/2013
Implementation of the basic finite volume scheme for solving the conservation law
u_t + f(u)_x = 0 x in [a,b]
u(x,0) = u_0(x)
u(a,t) = u(b,t) periodic boundary conditions :
meaning whatever exits the interval from one side it enters from the other side
(CAUTION : make sure the initial condition u0(x) is periodic in the [a,b])
The scheme is the one presented in the lecture
U_i^{n+1} = U_i^n - (Dt/Dx) ( F_{i+1/2) - F_{i-1/2} )
where F_{i+1/2} = F(U_i, U_{i+1})
Feel free to experiment with the code by changing the various parameters.
Try to implement the other numerical flux functions from the notes
The code ouputs the following
1) a series of files named Uxxxx.dat. Each Uxxxx.dat contains a snapshot of the solution at
a certain time step (look inside the Ustats.dat file for more info
2) a statistics file named Ustats.dat with a lot of information about the run.
if a line contains a star in its begining then at that time-instance a Uxxxx.dat was created
*/
#include
#include
#include
#define min(a,b) ((a)<=(b) ? (a) : (b))
#define max(a,b) ((a)>=(b) ? (a) : (b))
// The parameters of the problem
#define a -1.0 // Left endpoint of the interval
#define b 1.0 // Right endpoint of the interval
#define N 100 // Number of subintervals
#define CFL 1.0 // your CFL stability constant
#define T 1.0 // the final time
#define Nout 10 // Output of the solution every Nout-steps
double dt, dx;
// The continuous flux : the Burger's flux
double f(double u){
return u*u/2.0;
}
// The derivative of the continuous flux
double df(double u){
return u;
}
// Numerical flux : Lax-Friedrichs numerical flux
double F(double UR, double UL){
return (f(UR)+f(UL))/2.0 - (0.5*dt/dx)*(UL - UR);
}
// The initial condition u_0 (make sure it's periodic in [a,b]
double u0(double x) {
return cos(M_PI*x);
}
// The primitive of the initial condition to be used
// for computing the initial cell averages
double pu0(double x) {
return sin(M_PI*x)/M_PI;
}
void main(){
double t, dtodx, wspeed=0.0, tv=0.0, mass=0.0;
// XN : holds the mesh
// Un : the solution at the current step
// Up : the solution at the previous step
double X[N], Up[N], Un[N];
int i,nstep, nfid, fid;
char fname[]="Uxxxx.dat", digit[] = "0123456789"; // output filename
const char s9[]=" ";
FILE *fp, *fps;
dt = 1.0e-5; // Initial time-step
t = 0.0; // initial time
dx = (b-a)/N; // uniform mesh
nstep=0; // number of steps
fid = 0; // output file id
fps = fopen("Ustats.dat","w"); // Open the Ustats.dat file
fprintf(fps," Nstep %s t %s Dt %s Wspeed %s TV %s Mass\n", s9,s9,s9,s9,s9,s9);
for (i=0; i<=N; i++) X[i] = a + ((double)i) * dx; // the mesh
for (i=0; iT) {t-=dt; dt=T-t; t=T;} // check if you went over the final time and adjust appropriately
nstep++; // count the steps of the method
dtodx = dt / dx;
wspeed = 0.0; // wspeed = max_i |f'(U_i)| (for CFL condition, look at the notes)
// Now the scheme !!! (one line code)
for (i=0; i