/***********************************************************************/ /* LICENSED MATERIALS - PROPERTY OF IBM */ /* "RESTRICTED MATERIALS OF IBM" */ /* */ /* 5765-422 */ /* 5765-F84 */ /* 5765-G18 */ /* (C) COPYRIGHT IBM CORP. 1995, 2003. ALL RIGHTS RESERVED. */ /* */ /* U.S. GOVERNMENT USERS RESTRICTED RIGHTS - USE, DUPLICATION */ /* OR DISCLOSURE RESTRICTED BY GSA ADP SCHEDULE CONTRACT WITH */ /* IBM CORP. */ /***********************************************************************/ #include "diffus.h" #include #include #include /* first list all external variables */ double dif_ly_ratio; int dif_nx, dif_ny, dif_npts, dif_ntemps; double dif_delta_t; double *dif_x=NULL, *dif_y=NULL; /***************************************************************************** dif_ly_ratio is the ratio of the x and y lengths of the beam dif_nx is the number of sine expansion coefficients to use in the x direction dif_ny is the number of sine expansion coefficients to use in the y direction dif_delta_t is the size of the time step to be display on output dif_ntemps is the total number of temperatures to display per point dif_npts is the total number of points to output dif_x is the x coordinates of the points dif_y is the y coordinates of the points *****************************************************************************/ /* now include all the private variables */ static int init_f=1, diff_f=1; /***************************************************************************** init_f chooses the functional form of initial distribution of temperature diff_f chooses the functional form for spatially dependent head diffusion coefficient *****************************************************************************/ void init_diffusion( void) { int numx=5, numy=5, nx=7, ny=7, numt=20; double ly_ratio=1.0, delta_t=0.1; double delx, dely; int i, j, ij; FILE *in; char *path="diffus.dat"; char *pline; char line[256]; /*==========================================================================*/ /* Start of executable code */ /* if the diffus.dat file exists and is readable, override the defaults */ in = fopen(path,"r"); if ( in != NULL ) { pline = fgets(line,256,in); if (pline) ly_ratio = atof(line); /* read in ly_ration */ pline = fgets(line,256,in); if (pline) delta_t = atof(line); /* read in delta_t */ pline = fgets(line,256,in); if (pline) nx = atoi(line); /* read in polynomials in x direction */ pline = fgets(line,256,in); if (pline) ny = atoi(line); /* read in polynomials in y direction */ pline = fgets(line,256,in); if (pline) numx = atoi(line); /* read in output points in x direction */ pline = fgets(line,256,in); if (pline) numy = atoi(line); /* read in output points in y direction */ pline = fgets(line,256,in); if (pline) numt = atoi(line); /* read in output number of time steps */ pline = fgets(line,256,in); if (pline) init_f = atoi(line); /* choose initial temperature distribution */ pline = fgets(line,256,in); if (pline) diff_f = atoi(line); /* choose diffusion coefficient option */ } dif_ly_ratio = ly_ratio; dif_npts = numx*numy; dif_delta_t = delta_t; dif_ntemps = numt; dif_nx = nx; dif_ny = ny; /* get room for the x and y coordinates of our solution */ dif_x = (double *) malloc( numx*numy*sizeof(double)); dif_y = (double *) malloc( numx*numy*sizeof(double)); /* assign a simple linear array of points */ delx = PI/ ( numx + 1.0); dely = PI/ ( numy + 1.0); for( i=0; i < numx; i++) for( j=0; j < numy; j++) { ij = A_INDEX(i,j,numx); dif_x[ij] = delx * (double) (i+1); dif_y[ij] = dely * (double) (j+1); } return; } double init_temp( double x, double y) { /************************************************************************* the problem has been scaled to go from 0 to pi in both the x and y directions. To more easily calculate the expansion coefficients we define the function to be odd about pi and use the rang 0 < x < 2*pi x is the input x coordinate y is the input y coordinate init_temp is the initial temperature at that point **********************************************************************/ /* local variables */ int isign; double sign, val; isign = 1; /* force the function to be odd in x about pi */ if( x > PI ) { isign = -isign; x = TWOPI - x; } /* force the function to be odd in y about pi */ if( y > PI ) { isign = -isign; y = TWOPI - y; } sign = (double) isign; /* choose very simple temperature profile cases */ switch (init_f) { case 1: val = sign*(x*(PI-x))*y*(PI-y); break; case 2: val = sign*(x*(PI-x))*y*(PI-y)*y; case 3: val = sign*(x*(PI-x))*y*(PI-y)*x; case 4: val = sign*(x*(PI-x))*y*(PI-y)*x*y; case 5: val = sign*(x*(PI-x))*y*(PI-y); case 6: val = sign*(x*(PI-x))*(x*(PI-x)) *y*(PI-y); case 7: val = sign*(x*(PI-x))*(y*(PI-y))*(y*(PI-y)); default: val = sign*sin(x)*sin(y); } return val; } double diff_coef(double x, double y) { /************************************************************************* the problem has been scaled to go from 0 to pi in both the x and y directions. To simplify the matrix element calculations, we define the function to be even about pi x is the input x coordinate y is the input y coordinate diff_coef is value of the diffusion coefficient at (x,y) *************************************************************************/ double val; /* Start of executable code */ if( x > PI ) x = TWOPI - x; if( y > PI ) y = TWOPI - y; /* choose very simple diffusion coefficient cases */ switch (diff_f) { case 1: val = .50 + (x + y) / (2 * TWOPI); break; case 2: val = ((1.0 + x)*( PI -x + 1.0) *(y + PI))/ 3*PI; break; case 3: val = (y + PI) * PI/( ( PI + x) * (2* PI - x)); break; default: val = 1.0; } return val; }