/***********************************************************************/ /* LICENSED MATERIALS - PROPERTY OF IBM */ /* "RESTRICTED MATERIALS OF IBM" */ /* */ /* 5765-422 */ /* (C) COPYRIGHT IBM CORP. 1995. ALL RIGHTS RESERVED. */ /* */ /* U.S. GOVERNMENT USERS RESTRICTED RIGHTS - USE, DUPLICATION */ /* OR DISCLOSURE RESTRICTED BY GSA ADP SCHEDULE CONTRACT WITH */ /* IBM CORP. */ /***********************************************************************/ /******************************************************************** Purpose, to find the cooling rate for a specified set of points in an anisotropic retangular beam, immersed in a constant heat bath with a temperature of 0. *********************************************************************/ #include "diffus.h" int main( int argc, char **argv) { int n, ix, iy, k, i, j, it, ib, ig; int num_errors=0, *iwork,lwork,ilwork; /*********************************************************************** a is the sine expansion coefficients of the initial temperature profile b is the eigenvectors of the diffusion operator in the sine function basis ab is the initial temperature profile expanded in the eigenvectors of the diffusion operator f is the matrix elements of the diffusion operation in sine function basis lambda is the eigenvalues of the diffusion operator df is the Fourier transform of the diffusion coefficient function **************************************************************************/ double *lambda, **xg, *gap, *work; double **f, **b, **a; int f_d[DESC_DIM], b_d[DESC_DIM], a_d[DESC_DIM]; G_index *f_i, *b_i, *a_i; double **x, **ab, **abt; int x_d[DESC_DIM], ab_d[DESC_DIM], abt_d[DESC_DIM]; G_index *x_i, *ab_i, *abt_i; double **xsines, **ysines, **temp; int *ifail, *iclustr; int biga_index, num_eigvalues, num_vectors, info; /***************************************************************************** read in the problem size, initialize the problem dimensions, choose functional form for the spatially dependent heat diffusion coefficients, choose funtional form of initial temperature distribution and choose the number of points, in both space and time, of the solution to print out. *****************************************************************************/ init_diffusion(); /******************************************************************* read in the how many sine functions to include in both the x and y directions Since we need to get the Fourier coefficients of the sums and differences of the indices, we need to include twice as many Fourier coefficients as number of sine expansion coefficients. Also, since the top half of the Fourier transfor is bogus, an artifact of artificially extending the diffusion coefficient function and the initial temperature distribution, we need to double the number of Fourier coefficients again. Finally, the extra sum of one comes from the fact that the sine function expansion starts a 1 and the Fourier expansion starts at 0. ************************************************************************/ /* n is the order of the diffusion operator equation */ n = dif_nx * dif_ny; /* initialize BLACS and calculate default block sizes */ initialize_scale(n, &biga_index); /* allocate room for the eigenvalues of diffusion operator */ lambda = ( double *) malloc( sizeof(double) * n ); /* allocate room for sines of x and y coordinates */ allocate_rarray( &xsines, dif_npts, dif_nx); allocate_rarray( &ysines, dif_npts, dif_ny); /* allocate room for temperature history at selected points */ allocate_rarray( &temp, dif_npts, dif_ntemps); /* allocate room for global temperature profile expansion vector at time t */ allocate_rarray( &xg, dif_nx, dif_ny); /* allocate arrays for eigenvalue decomposition */ ifail = ( int * ) malloc( n * sizeof(int)); iclustr = ( int * ) malloc( n * sizeof(int)); gap = ( double * ) malloc( sc_nnodes * sizeof(double)); /* check for allocation errors */ if( lambda == NULL || xsines == NULL || ysines == NULL || temp == NULL || xg == NULL || ifail == NULL || iclustr == NULL || gap == NULL) { num_errors++; } igamx2d(sc_icontext,"A"," ",1,1,&num_errors,1, -1,-1,-1,-1,-1); if( num_errors > 0 ) if( sc_iam == 0 ) { printf("Error in allocating arrays in main\n"); blacs_abort(sc_icontext, 1); } /* get matrices */ /* diffusion operator matrix */ initialize_rarray(&f, f_d, &f_i, n, n, biga_index); /* eigenvectors of diffusion operator matrix */ initialize_rarray(&b, b_d, &b_i, n, n, biga_index); /* initial temperature profile, in row vector */ initialize_rarray(&a, a_d, &a_i, 1, n, biga_index); /* initial temperature profile, in eigenfunction basis, in row vector */ initialize_rarray(&ab, ab_d, &ab_i, 1, n, biga_index); /* temperature profile, at time t, in eigenfunction basis, in row vector */ initialize_rarray(&abt, abt_d, &abt_i, 1, n, biga_index); /* temperarure profile in at time t in sine expansion basis, in row vector */ initialize_rarray(&x, x_d, &x_i, 1, n, biga_index); /* represent initial temperature in sine function expansion */ expand_temp_profile(a,a_i,a_d,1); /* rlocal_to_rglobal(a, a_d, *xg, 1 ); if( sc_iam == 0 ) { for(i=0;inum_col_blks; ib++) for( ig=ab_i->g_start_col_blk[ib]; ig <= ab_i->g_end_col_blk[ib]; ig++) { i++; abt[i][0] = exp( -lambda[ig] * (it+1) * dif_delta_t) * ab[i][0]; } /* abt now has the expansion of the temperature profile terms of the */ /* eigenvectors of the diffusion operator */ /* multiply the eigenvector matrix with abt to give the sine function */ /* expansion of the temperature profile at time t, x */ pdgemv("N", n, n, 1.e0, *b, 1, 1, b_d, *abt, 1, 1, abt_d, 1, 0.e0, *x, 1, 1, x_d, 1); /* gather all of the local pieces of the array x to the array xg */ rlocal_to_rglobal(x, x_d, *xg, 1 ); for ( k=0; k< dif_npts; k++) temp[it][k] = 0.e0; for ( iy=0; iy< dif_ny; iy++) for ( ix=0; ix< dif_ny; ix++) for ( k=0; k< dif_npts; k++) temp[it][k] = temp[it][k] + xsines[ix][k] * ysines[iy][k] * xg[iy][ix]; } /* end of time loop */ if ( sc_iam == 0) { printf(" point # X Y\n"); for ( i=0; i< dif_npts; i++) printf(" %6d %11.4f %11.4f\n",i, dif_x[i], dif_y[i]); printf("\n"); for( k=0; k < dif_npts; k += 6) { printf("\n"); printf(" Points\n"); printf(" TIME "); for( i=k; i< min(k+6, dif_npts); i++) printf(" #%4d ",i); printf("\n"); for( i=0; i< dif_ntemps; i++) { printf("%11.5f", (double) (i+1) * dif_delta_t); for(j=k; j< min(k+6, dif_npts); j++) printf("%11.5f",temp[i][j]); printf("\n"); } } } blacs_barrier(sc_icontext,"A"); sleep(5); blacs_exit(0); free(lambda); free( ifail); free( iclustr); free( gap); return 0; }