kmazouzi / FormationMeso

Download as
Browse Code ยป

Commit 1e6ef8e72e286128c10aeb45273e3a2705c92230

Authored by kmazouzi
1 parent 6a98a5afa7
Exists in master

Steady state heat

Showing 2 changed files with 436 additions and 0 deletions Side-by-side Diff

lab3/omp_heat2D.c View file @ 1e6ef8e
  1 +/****************************************************************************
  2 + * DESCRIPTION:
  3 + * Serial HEAT2D Example - C Version
  4 + * This example is based on a simplified
  5 + * two-dimensional heat equation domain decomposition. The initial
  6 + * temperature is computed to be high in the middle of the domain and
  7 + * zero at the boundaries. The boundaries are held at zero throughout
  8 + * the simulation. During the time-stepping, an array containing two
  9 + * domains is used; these domains alternate between old data and new data.
  10 + *
  11 + * The physical region, and the boundary conditions, are suggested
  12 + by this diagram;
  13 +
  14 + u = 0
  15 + +------------------+
  16 + | |
  17 + u = 100 | | u = 100
  18 + | |
  19 + +------------------+
  20 + u = 100
  21 +
  22 +Interrior point :
  23 + u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
  24 +
  25 + ****************************************************************************/
  26 +#include <stdio.h>
  27 +#include <stdlib.h>
  28 +#include <math.h>
  29 +#include <omp.h>
  30 +
  31 +#define N 500
  32 +#define M 500
  33 +
  34 +#define ITER_PRINT 100
  35 +#define PRINT_DATA 1
  36 +
  37 +#define EPSILON 1e-1
  38 +
  39 +
  40 +void update(int nx,int ny, float *u, float *unew, float * diff);
  41 +void inidat(int nx, int ny, float *u, float *unew);
  42 +void prtdat(int nx, int ny, float *u,const char *fnam);
  43 +
  44 +
  45 +
  46 +
  47 +int main(int argc, char *argv[])
  48 +{
  49 +
  50 + float diff=1.0;
  51 +
  52 + float *u = (float *)malloc(N * M * sizeof(float));
  53 + float *unew = (float *)malloc(N * M * sizeof(float));
  54 +
  55 + if(u==0 || unew ==0)
  56 + {
  57 + perror("Can't allocated data\n");
  58 + return -1;
  59 + }
  60 +
  61 + printf ( "\n" );
  62 + printf ( "HEATED_PLATE\n" );
  63 + printf ( " Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() );
  64 + printf ( " A program to solve for the steady state temperature distribution\n" );
  65 + printf ( " over a rectangular plate.\n" );
  66 + printf ( " Spatial grid of %d by %d points.\n\n", M, N );
  67 +
  68 +
  69 + /* Initialize grid and create input file */
  70 + printf("Initializing grid\n");
  71 +
  72 + inidat(N, M,u,unew);
  73 +
  74 + prtdat(N, M,u, "initial.dat");
  75 +
  76 +
  77 + printf("Start computing\n");
  78 +
  79 + int iter=0;
  80 +
  81 + /*
  82 + * iterate until the new solution unew differs from the old solution u
  83 + * by no more than EPSILON.
  84 + * */
  85 +
  86 + while(diff> EPSILON) {
  87 +
  88 + update(N, M, u, unew,&diff);
  89 +
  90 + if(iter%ITER_PRINT==0)
  91 + printf("Iteration %d, diff = %f\n ", iter,diff);
  92 +
  93 + iter++;
  94 + }
  95 +
  96 + prtdat(N, M, u, "final.dat");
  97 +
  98 + free(u);
  99 + free(unew);
  100 +}
  101 +
  102 +
  103 +
  104 +/****************************************************************************
  105 + * subroutine update
  106 + ****************************************************************************/
  107 +void update(int nx,int ny, float *u, float *unew, float * diff)
  108 +{
  109 + int ix, iy;
  110 + *diff=0.0;
  111 +
  112 +#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy)
  113 + for (ix = 1; ix < nx-1; ix++) {
  114 + for (iy = 1; iy < ny-1; iy++) {
  115 + unew[ix*ny+iy] =
  116 + (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] +
  117 + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0;
  118 +
  119 + }
  120 + }
  121 +
  122 +//compute reduction
  123 +
  124 +
  125 + float mydiff;
  126 +
  127 +#pragma omp parallel shared(nx,ny,u,unew, diff) private (ix,iy,mydiff)
  128 + {
  129 +mydiff=0.0;
  130 +#pragma omp for
  131 + for (ix = 1; ix < nx-1; ix++) {
  132 + for (iy = 1; iy < ny-1; iy++) {
  133 + if (mydiff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
  134 + {
  135 + mydiff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
  136 + }
  137 + }
  138 + }
  139 +
  140 +
  141 +# pragma omp critical
  142 + {
  143 + if (*diff < mydiff )
  144 + {
  145 + *diff = mydiff;
  146 + }
  147 + }
  148 +
  149 +
  150 +#pragma omp for
  151 + for (ix = 1; ix < nx-1; ix++) {
  152 + for (iy = 1; iy < ny-1; iy++) {
  153 + u[ix*ny+iy] = unew[ix*ny+iy];
  154 + }
  155 + }
  156 +
  157 +
  158 +
  159 +
  160 + }
  161 +}
  162 +
  163 +/*****************************************************************************
  164 + * Initialize Data
  165 + *****************************************************************************/
  166 +void inidat(int nx, int ny, float *u, float *unew)
  167 +{
  168 + int ix, iy;
  169 +
  170 + /*
  171 + *Set boundary data and interrior values
  172 + * */
  173 + for (ix = 0; ix < nx; ix++)
  174 + for (iy = 0; iy < ny; iy++) {
  175 +
  176 + if(ix==0)
  177 + {
  178 + u[ix*ny+iy]=0.0;
  179 + }
  180 + else
  181 + if(iy==0 && ix!=0)
  182 + {
  183 + u[ix*ny+iy]=100.0;
  184 + }else
  185 +
  186 + if(ix==nx-1)
  187 + {
  188 + u[ix*ny+iy]=100.0;
  189 + }else
  190 +
  191 + if(iy==ny-1 && ix!=0)
  192 + {
  193 + u[ix*ny+iy]=100.0;
  194 + }else
  195 +
  196 + u[ix*ny+iy]=( float ) ( 2 * nx + 2 * ny - 4 );
  197 + }
  198 +}
  199 +
  200 +/**************************************************************************
  201 + * Print Data to files
  202 + **************************************************************************/
  203 +void prtdat(int nx, int ny, float *u,const char *fnam)
  204 +{
  205 +
  206 + int ix, iy;
  207 + FILE *fp;
  208 +
  209 + if(ITER_PRINT==0)return;
  210 +
  211 + fp = fopen(fnam, "w");
  212 +
  213 + for (ix = 0 ; ix < nx; ix++) {
  214 + for (iy =0; iy < ny; iy++) {
  215 +
  216 + fprintf(fp, "%8.3f", u[ix*ny+iy]);
  217 +
  218 + if(iy!=ny-1)
  219 + {
  220 + fprintf(fp, " ");
  221 + }else
  222 + {
  223 + fprintf(fp, "\n");
  224 + }
  225 + }
  226 + }
  227 +
  228 + fclose(fp);
  229 +}
lab3/ser_heat2D.c View file @ 1e6ef8e
  1 +/****************************************************************************
  2 + * DESCRIPTION:
  3 + * Serial HEAT2D Example - C Version
  4 + * This example is based on a simplified
  5 + * two-dimensional heat equation domain decomposition. The initial
  6 + * temperature is computed to be high in the middle of the domain and
  7 + * zero at the boundaries. The boundaries are held at zero throughout
  8 + * the simulation. During the time-stepping, an array containing two
  9 + * domains is used; these domains alternate between old data and new data.
  10 + *
  11 + * The physical region, and the boundary conditions, are suggested
  12 + by this diagram;
  13 +
  14 + u = 0
  15 + +------------------+
  16 + | |
  17 + u = 100 | | u = 100
  18 + | |
  19 + +------------------+
  20 + u = 100
  21 +
  22 +Interrior point :
  23 + u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
  24 +
  25 + ****************************************************************************/
  26 +#include <stdio.h>
  27 +#include <stdlib.h>
  28 +#include <math.h>
  29 +
  30 +#define N 200
  31 +#define M 200
  32 +
  33 +#define ITER_PRINT 100
  34 +#define PRINT_DATA 1
  35 +
  36 +#define EPSILON 1e-1
  37 +
  38 +
  39 +void update(int nx,int ny, float *u, float *unew, float * diff);
  40 +void inidat(int nx, int ny, float *u, float *unew);
  41 +void prtdat(int nx, int ny, float *u,const char *fnam);
  42 +
  43 +
  44 +
  45 +
  46 +int main(int argc, char *argv[])
  47 +{
  48 +
  49 + float diff=1.0;
  50 +
  51 + float *u = (float *)malloc(N * M * sizeof(float));
  52 + float *unew = (float *)malloc(N * M * sizeof(float));
  53 +
  54 + if(u==0 || unew ==0)
  55 + {
  56 + perror("Can't allocated data\n");
  57 + return -1;
  58 + }
  59 +
  60 + printf ( "\n" );
  61 + printf ( "HEATED_PLATE\n" );
  62 + printf ( " Serial version\n" );
  63 + printf ( " A program to solve for the steady state temperature distribution\n" );
  64 + printf ( " over a rectangular plate.\n" );
  65 + printf ( "\n" );
  66 + printf ( " Spatial grid of %d by %d points.\n\n", M, N );
  67 +
  68 +
  69 + /* Initialize grid and create input file */
  70 + printf("Initializing grid\n");
  71 +
  72 + inidat(N, M,u,unew);
  73 +
  74 + prtdat(N, M,u, "initial.dat");
  75 +
  76 +
  77 + printf("Start computing\n");
  78 +
  79 + int iter=0;
  80 +
  81 + /*
  82 + * iterate until the new solution unew differs from the old solution u
  83 + * by no more than EPSILON.
  84 + * */
  85 +
  86 + while(diff> EPSILON) {
  87 +
  88 + update(N, M, u, unew,&diff);
  89 +
  90 + if(iter%ITER_PRINT==0)
  91 + printf("Iteration %d, diff = %f\n ", iter,diff);
  92 +
  93 + iter++;
  94 + }
  95 +
  96 + prtdat(N, M, u, "final.dat");
  97 +
  98 + free(u);
  99 + free(unew);
  100 +}
  101 +
  102 +
  103 +
  104 +/****************************************************************************
  105 + * subroutine update
  106 + ****************************************************************************/
  107 +void update(int nx,int ny, float *u, float *unew, float * diff)
  108 +{
  109 + int ix, iy;
  110 + *diff=0.0;
  111 +
  112 + for (ix = 1; ix < nx-1; ix++) {
  113 + for (iy = 1; iy < ny-1; iy++) {
  114 + unew[ix*ny+iy] =
  115 + (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] +
  116 + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0;
  117 +
  118 + if (*diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
  119 + {
  120 + *diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
  121 + }
  122 + }
  123 +
  124 + }
  125 +
  126 +
  127 + for (ix = 1; ix < nx-1; ix++) {
  128 + for (iy = 1; iy < ny-1; iy++) {
  129 + u[ix*ny+iy] = unew[ix*ny+iy];
  130 + }
  131 + }
  132 +
  133 +}
  134 +
  135 +/*****************************************************************************
  136 + * Initialize Data
  137 + *****************************************************************************/
  138 +void inidat(int nx, int ny, float *u, float *unew)
  139 +{
  140 + int ix, iy;
  141 +
  142 + /*
  143 + *Set boundary data and interrior values
  144 + * */
  145 + for (ix = 0; ix < nx; ix++)
  146 + for (iy = 0; iy < ny; iy++) {
  147 +
  148 + if(ix==0)
  149 + {
  150 + u[ix*ny+iy]=0.0;
  151 + }
  152 + else
  153 + if(iy==0 && ix!=0)
  154 + {
  155 + u[ix*ny+iy]=100.0;
  156 + }else
  157 +
  158 + if(ix==nx-1)
  159 + {
  160 + u[ix*ny+iy]=100.0;
  161 + }else
  162 +
  163 + if(iy==ny-1 && ix!=0)
  164 + {
  165 + u[ix*ny+iy]=100.0;
  166 + }else
  167 +
  168 + u[ix*ny+iy]=( float ) ( 2 * nx + 2 * ny - 4 );
  169 + }
  170 +}
  171 +
  172 +/**************************************************************************
  173 + * Print Data to files
  174 + **************************************************************************/
  175 +void prtdat(int nx, int ny, float *u,const char *fnam)
  176 +{
  177 +
  178 + int ix, iy;
  179 + FILE *fp;
  180 +
  181 + if(ITER_PRINT==0)return;
  182 +
  183 + fp = fopen(fnam, "w");
  184 +
  185 + for (ix = 0 ; ix < nx; ix++) {
  186 + for (iy =0; iy < ny; iy++) {
  187 +
  188 + fprintf(fp, "%8.3f", u[ix*ny+iy]);
  189 +
  190 + if(iy!=ny-1)
  191 + {
  192 + fprintf(fp, " ");
  193 + }else
  194 + {
  195 + fprintf(fp, "\n");
  196 + }
  197 + }
  198 + }
  199 +
  200 + fclose(fp);
  201 +}