[yosuke1984]

Hi all

I'm testing discritisation schemes for advection term, but struggling. Some how, second upwind and third upwind method diverges. I am using identical grid arrangement and courant conditions for all, and others ( quick, first upwind, kk-scheme) worked.

My reference says, higher order scheme gives smaller numerical diffusion but from what I can see from the result, it's opposite. The discritised equations were worked out from the expression on wikipedia and confirmed the reference as well.

Can any body point me to the right direction.

Here is my code and the result.


http://www.geocities.jp/yosuke_31415...rst_upwind.gif
fig 1: first upwind


http://www.geocities.jp/yosuke_3141592/pict/second_upwind.gif
fig 2: second upwind

Best Regards

y
---------------------------------------------------------------------------------------------------------------
#include<stdio.h>
#include<math.h>
#define NODES 128
#define V 2.0;

void init( const double *dt, const double *dx, double *t, double *u, double *f );
void advec_diff( FILE *wfp, const double *dt, const double *dx, const double *nu, double *u, double *f );
double second_upwind( int *i, double *f, double *u, const double *dx );

int main(){
char filename[128];
double t=0;
const double t_end = 0.5;
const double dt=0.01;
const double dx=0.1;
const double nu=0.0;
double u[NODES] ;
double f[NODES] ;

init( &dt, &dx, &t, u, f);

printf("C=%lf\n", u[0]*dt/dx);

FILE *wfp;
for( t=0; t<t_end; t+=dt){
sprintf( filename, "result-%lf.data", t);
if( (wfp=fopen( filename, "w")) == NULL ){
printf("File write ERROR\n");
}
advec_diff( wfp, &dt, &dx, &nu, u, f );
}
fclose(wfp);
return 0;}

void init( const double *dt, const double *dx, double *t, double *u, double *f ){
int i;
*t = 0.0;
for( i=0; i<NODES; i++){
u[i] = V;
f[i] = 0.0;
}
f[NODES/2-1] = 1.0;
f[NODES/2-2] = 1.0;
f[NODES/2-3] = 1.0;
f[NODES/2-4] = 1.0;
f[NODES/2-5] = 1.0;
}

void advec_diff( FILE *wfp, const double *dt, const double *dx, const double *nu, double *u, double *f ){
int i;
double dmy_f[NODES];
for( i=2; i<NODES-2; i++){
dmy_f[i] = f[i] + *dt * ( -second_upwind( &i, f, u, &*dx ) + *nu*( f[i-1] - 2.0*f[i] + f[i+1]) /( *dx**dx ) );
fprintf( wfp, " %lf\n", dmy_f[i]);
}
for( i=2; i<NODES-2; i++){
f[i] = dmy_f[i];
}
}

double second_upwind( int *i, double *f, double *u, const double *dx ){
double tmp;
tmp = u[*i]*( f[*i-2]-4.0*f[*i-1]+4.0*f[*i+1]-f[*i+2] )/(4.0**dx);
- fabs( u[*i] )*( -f[*i-2]+4.0*f[*i-1]-6.0*f[*i]+4.0*f[*i+1]-f[*i+2] )/ (4.0**dx);
return tmp;
}