/* expdiff.c explicit diffusion simulation J. Watlington, 10/16/95 */ #include #include #include #define DEBUG #define DEBUG_PERIOD 20.0 #define NUM_POINTS 500 #define NUM_BUFFERS 2 #define MIN_TOTAL_CHANGE 0.001 #define MAX_TIME 200.0 #define DIFF_COEFF 1.0 #define DELTA_X 1.0 #define DEFAULT_DELTA_T 1.0 #define START_IMPULSE_AMPL 1.0 #define START_IMPULSE_LOC ( NUM_POINTS / 2 ) typedef long double sim_type; void display_lattice( sim_type *lattice ); void main( int argc, char **argv ) { int x; int buf; int cur_buf = 1; int prev_buf = 0; long double t = 0.0; long double delta_t = DEFAULT_DELTA_T; long double total_change = 10.0 * MIN_TOTAL_CHANGE; sim_type present_value; sim_type new_value; sim_type state[ NUM_BUFFERS ][ NUM_POINTS + 2 ]; /* Parse arguments, if provided */ if( argc > 2 ) { printf( "usage: %s []\n", argv[0] ); exit( -1 ); } else if( argc > 1 ) { sscanf( argv[1], "%Lf", &delta_t ); } printf( "running till total change less than %f\n", MIN_TOTAL_CHANGE ); printf( "using a %f second time step\n", delta_t ); /* Initialize the simulation state to initial conditions */ for( buf = 0; buf < NUM_BUFFERS; buf++ ) for( x = 0; x < (NUM_POINTS + 2); x++ ) state[ buf ][ x ] = 0.0; state[ prev_buf ][ START_IMPULSE_LOC ] = START_IMPULSE_AMPL; #ifdef DEBUG printf( "\nAt time t = %f, here are the results :\n", t ); display_lattice( state[ prev_buf ] ); #endif /* And loop until the simulation stops changing... */ while( (total_change > MIN_TOTAL_CHANGE) && (t < MAX_TIME) ) { /* Perform one iteration through the lattice, calculating the new values. */ for( total_change = 0.0, x = 1; x <= NUM_POINTS; x++ ) { present_value = state[ prev_buf ][ x ]; new_value = present_value + (((DIFF_COEFF * delta_t)/(DELTA_X * DELTA_X)) * ( state[ prev_buf ][ x - 1 ] - (2.0 * present_value) + state[ prev_buf ][ x + 1 ] )); state[ cur_buf ][ x ] = new_value; total_change += fabs( new_value - present_value ); } /* The edge values at x = 0 and x = NUM_POINTS+1 are mirrors of the data value next to them, simulating a perfectly insulated boundary. */ state[ cur_buf ][ 0 ] = state[ cur_buf ][ 1 ]; state[ cur_buf ][ NUM_POINTS + 1 ] = state[ cur_buf ][ NUM_POINTS ]; /* Flip buffers */ if( ++cur_buf >= NUM_BUFFERS ) cur_buf = 0; if( ++prev_buf >= NUM_BUFFERS ) prev_buf = 0; t += delta_t; #ifdef DEBUG if( fmod( t, DEBUG_PERIOD ) == 0.0 ) { printf( "\nAt time t = %f, here are the results :\n", t ); display_lattice( state[ prev_buf ] ); } #endif } if( t >= MAX_TIME ) printf( "The solution never stabilized !!!\n" ); printf( "\nAt time t = %f, here are the results :\n", t ); display_lattice( state[ prev_buf ] ); } #define NUM_COLS 10 #define NUM_ROWS ( NUM_POINTS / NUM_COLS ) void display_lattice( sim_type *lattice ) { int y; int x; int index; for( index = 1, y = 0; y < NUM_ROWS; y++ ) { printf( "%d:", (NUM_COLS * y) ); for( x = 0; x < NUM_COLS; x++, index++ ) { printf( " %1.4f", lattice[ index ] ); } printf( "\n" ); } }