kmazouzi / FormationMeso

Compare View

switch
from
...
to
 
Commits (4)

Diff

Showing 4 changed files Side-by-side Diff

lab3/Makefile View file @ 2b4ce20
... ... @@ -0,0 +1,29 @@
  1 +GCC = gcc
  2 +CFLAGS = -O3 -fopenmp
  3 +OMP_FLAG = -fopenmp
  4 +RM = rm -rf
  5 +
  6 +
  7 +EXE = omp_heat2D ser_heat2D
  8 +
  9 +all : $(EXE)
  10 +
  11 +#.PHONY: all clean purge
  12 +
  13 +
  14 +pi_ser: ser_heat2D.o
  15 + $(GCC) $(CFLAGS) -o $@ $^
  16 +
  17 +pi_task: omp_heat2D.o
  18 + $(GCC) $(CFLAGS) -o $@ $^
  19 +
  20 +
  21 +%.o :%.c
  22 + $(GCC) $(CFLAGS) -c -o $@ $<
  23 +
  24 +
  25 +clean:
  26 + $(RM) *.o *.dat
  27 +
  28 +purge: clean
  29 + $(RM) $(EXE)
lab3/omp_heat2D.c View file @ 2b4ce20
... ... @@ -0,0 +1,240 @@
  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 NN 50
  32 +#define MM 50
  33 +
  34 +#define ITER_PRINT 100
  35 +#define PRINT_DATA 1
  36 +
  37 +#define _EPSILON 0.001
  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 + float EPSILON=_EPSILON;
  52 + int N=NN,M=MM;
  53 +
  54 + if(argc !=3)
  55 + {
  56 + fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
  57 + fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
  58 + fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
  59 + return -1;
  60 + }
  61 +
  62 + N = M = atoi(argv[1]);
  63 + EPSILON = atof(argv[2]);
  64 +
  65 + float *u = (float *)malloc(N * M * sizeof(float));
  66 + float *unew = (float *)malloc(N * M * sizeof(float));
  67 +
  68 + if(u==0 || unew ==0)
  69 + {
  70 + perror("Can't allocated data\n");
  71 + return -1;
  72 + }
  73 +
  74 + printf ( "\n" );
  75 + printf ( "HEATED_PLATE\n" );
  76 + printf ( " Parallel OpenMP version, using %d Threads\n",omp_get_max_threads() );
  77 + printf ( " A program to solve for the steady state temperature distribution\n" );
  78 + printf ( " over a rectangular plate.\n" );
  79 + printf ( " Spatial grid of %d by %d points.\n\n", M, N );
  80 +
  81 +
  82 + /* Initialize grid and create input file */
  83 + printf("Initializing grid\n");
  84 +
  85 + inidat(N, M,u,unew);
  86 +
  87 + prtdat(N, M,u, "initial.dat");
  88 +
  89 + printf("Start computing\n");
  90 +
  91 + int iter=0;
  92 +
  93 + /*
  94 + * iterate until the new solution unew differs from the old solution u
  95 + * by no more than EPSILON.
  96 + * */
  97 +
  98 + while(diff> EPSILON) {
  99 +
  100 + update(N, M, u, unew,&diff);
  101 +
  102 + if(iter%ITER_PRINT==0)
  103 + printf("Iteration %d, diff = %f\n ", iter,diff);
  104 +
  105 + iter++;
  106 + }
  107 +
  108 + prtdat(N, M, u, "final.dat");
  109 +
  110 + free(u);
  111 + free(unew);
  112 +}
  113 +
  114 +
  115 +
  116 +/****************************************************************************
  117 + * subroutine update
  118 + ****************************************************************************/
  119 +void update(int nx,int ny, float *u, float *unew, float * diff)
  120 +{
  121 + int ix, iy;
  122 + *diff=0.0;
  123 +
  124 +#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy)
  125 + for (ix = 1; ix < nx-1; ix++) {
  126 + for (iy = 1; iy < ny-1; iy++) {
  127 + unew[ix*ny+iy] =
  128 + (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] +
  129 + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0;
  130 +
  131 + }
  132 + }
  133 +
  134 + //compute reduction
  135 +
  136 +
  137 + float mydiff;
  138 +
  139 +#pragma omp parallel shared(nx,ny,u,unew, diff) private (ix,iy,mydiff)
  140 + {
  141 + mydiff=0.0;
  142 +#pragma omp for
  143 + for (ix = 1; ix < nx-1; ix++) {
  144 + for (iy = 1; iy < ny-1; iy++) {
  145 + if (mydiff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
  146 + {
  147 + mydiff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
  148 + }
  149 + }
  150 + }
  151 +
  152 +
  153 +# pragma omp critical
  154 + {
  155 + if (*diff < mydiff )
  156 + {
  157 + *diff = mydiff;
  158 + }
  159 + }
  160 +
  161 +
  162 +#pragma omp for
  163 + for (ix = 1; ix < nx-1; ix++) {
  164 + for (iy = 1; iy < ny-1; iy++) {
  165 + u[ix*ny+iy] = unew[ix*ny+iy];
  166 + }
  167 + }
  168 + }
  169 +}
  170 +
  171 +/*****************************************************************************
  172 + * Initialize Data
  173 + *****************************************************************************/
  174 +void inidat(int nx, int ny, float *u, float *unew)
  175 +{
  176 + int ix, iy;
  177 +
  178 + /*
  179 + *Set boundary data and interrior values
  180 + * */
  181 + for (ix = 0; ix < nx; ix++)
  182 + for (iy = 0; iy < ny; iy++) {
  183 +
  184 + if(ix==0)
  185 + {
  186 + u[ix*ny+iy]=0.0;
  187 + }
  188 + else
  189 + if(iy==0 && ix!=0)
  190 + {
  191 + u[ix*ny+iy]=100.0;
  192 + }else
  193 +
  194 + if(ix==nx-1)
  195 + {
  196 + u[ix*ny+iy]=100.0;
  197 + }else
  198 +
  199 + if(iy==ny-1 && ix!=0)
  200 + {
  201 + u[ix*ny+iy]=100.0;
  202 + }else
  203 +
  204 + u[ix*ny+iy]=0.0;
  205 + }
  206 +}
  207 +
  208 +/**************************************************************************
  209 + * Print Data to files
  210 + **************************************************************************/
  211 +void prtdat(int nx, int ny, float *u,const char *fnam)
  212 +{
  213 +
  214 + int ix, iy;
  215 + FILE *fp;
  216 +
  217 + if(ITER_PRINT==0)return;
  218 +
  219 + fp = fopen(fnam, "w");
  220 +
  221 + for (ix = 0 ; ix < nx; ix++) {
  222 + for (iy =0; iy < ny; iy++) {
  223 +
  224 + fprintf(fp, "%8.3f", u[ix*ny+iy]);
  225 +
  226 + if(iy!=ny-1)
  227 + {
  228 + fprintf(fp, " ");
  229 + }else
  230 + {
  231 + fprintf(fp, "\n");
  232 + }
  233 + }
  234 + }
  235 +
  236 + fclose(fp);
  237 +}
  238 +
  239 +
  240 +
lab3/plot.py View file @ 2b4ce20
... ... @@ -0,0 +1,18 @@
  1 +#!/usr/bin/python
  2 +#-*-coding:utf8-*-
  3 +from sys import argv
  4 +from pylab import *
  5 +
  6 +if __name__ == "__main__":
  7 + if len(argv) != 2:
  8 + print "Usage: python", argv[0], " matrix.dat"
  9 + exit (0)
  10 +
  11 +
  12 +f = open ( argv[1] , 'r')
  13 +A = [ map(float,line.split()) for line in f ]
  14 +
  15 +figure(1)
  16 +imshow(A, interpolation='nearest')
  17 +grid(True)
  18 +show()
lab3/ser_heat2D.c View file @ 2b4ce20
... ... @@ -0,0 +1,213 @@
  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 NN 50
  31 +#define MM 50
  32 +
  33 +#define ITER_PRINT 100
  34 +#define PRINT_DATA 1
  35 +
  36 +#define _EPSILON 0.01
  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 + int N=NN,M=MM;
  49 +
  50 + float EPSILON=_EPSILON;
  51 +
  52 + if(argc !=3)
  53 + {
  54 + fprintf(stderr,"usage %s N EPSILON\n ", argv[0]);
  55 + fprintf(stderr,"\t\twhere N is GRID size, EPSILON is Tolerance\n");
  56 + fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1\n");
  57 + return -1;
  58 + }
  59 +
  60 + N = M = atoi(argv[1]);
  61 + EPSILON = atof(argv[2]);
  62 +
  63 + float diff=1.0;
  64 +
  65 + float *u = (float *)malloc(N * M * sizeof(float));
  66 + float *unew = (float *)malloc(N * M * sizeof(float));
  67 +
  68 + if(u==0 || unew ==0)
  69 + {
  70 + perror("Can't allocated data\n");
  71 + return -1;
  72 + }
  73 +
  74 + printf ( "\n" );
  75 + printf ( "HEATED_PLATE\n" );
  76 + printf ( " Serial version\n" );
  77 + printf ( " A program to solve for the steady state temperature distribution\n" );
  78 + printf ( " over a rectangular plate.\n" );
  79 + printf ( "\n" );
  80 + printf ( " Spatial grid of %d by %d points.\n", M, N );
  81 + printf ( " The iteration will end until tolerance <= %f\n\n",EPSILON);
  82 +
  83 + /* Initialize grid and create input file */
  84 + printf("Initializing grid\n");
  85 +
  86 + inidat(N, M,u,unew);
  87 +
  88 + prtdat(N, M,u, "initial.dat");
  89 +
  90 + printf("Start computing\n\n");
  91 +
  92 + int iter=0;
  93 +
  94 + /*
  95 + * iterate until the new solution unew differs from the old solution u
  96 + * by no more than EPSILON.
  97 + * */
  98 +
  99 + while(diff> EPSILON) {
  100 +
  101 + update(N, M, u, unew,&diff);
  102 +
  103 + if(iter%ITER_PRINT==0)
  104 +
  105 + printf("Iteration %d, diff = %f\n ", iter,diff);
  106 +
  107 + iter++;
  108 + }
  109 +
  110 + prtdat(N, M, u, "final.dat");
  111 +
  112 + free(u);
  113 + free(unew);
  114 +}
  115 +
  116 +
  117 +
  118 +/****************************************************************************
  119 + * subroutine update
  120 + ****************************************************************************/
  121 +void update(int nx,int ny, float *u, float *unew, float * diff)
  122 +{
  123 + int ix, iy;
  124 + *diff=0.0;
  125 +
  126 + for (ix = 1; ix < nx-1; ix++) {
  127 + for (iy = 1; iy < ny-1; iy++) {
  128 + unew[ix*ny+iy] = (u[(ix+1)*ny+iy] + u[(ix-1)*ny+iy] + u[ix*ny+iy+1] + u[ix*ny+iy-1] )/4.0
  129 + ;
  130 + if (*diff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
  131 + {
  132 + *diff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
  133 + }
  134 + }
  135 +
  136 + }
  137 +
  138 +
  139 + for (ix = 1; ix < nx-1; ix++) {
  140 + for (iy = 1; iy < ny-1; iy++) {
  141 + u[ix*ny+iy] = unew[ix*ny+iy];
  142 + }
  143 + }
  144 +
  145 +}
  146 +
  147 +/*****************************************************************************
  148 + * Initialize Data
  149 + *****************************************************************************/
  150 +void inidat(int nx, int ny, float *u, float *unew)
  151 +{
  152 + int ix, iy;
  153 +
  154 + /*
  155 + *Set boundary data and interrior values
  156 + * */
  157 + for (ix = 0; ix < nx; ix++)
  158 + for (iy = 0; iy < ny; iy++) {
  159 +
  160 + if(ix==0)
  161 + {
  162 + u[ix*ny+iy]=0.0;
  163 + }
  164 + else
  165 + if(iy==0 && ix!=0)
  166 + {
  167 + u[ix*ny+iy]=100.0;
  168 + }else
  169 +
  170 + if(ix==nx-1)
  171 + {
  172 + u[ix*ny+iy]=100.0;
  173 + }else
  174 +
  175 + if(iy==ny-1 && ix!=0)
  176 + {
  177 + u[ix*ny+iy]=100.0;
  178 + }else
  179 +
  180 + u[ix*ny+iy]=0.0;
  181 + }
  182 +}
  183 +
  184 +/**************************************************************************
  185 + * Print Data to files
  186 + **************************************************************************/
  187 +void prtdat(int nx, int ny, float *u,const char *fname)
  188 +{
  189 +
  190 + int ix, iy;
  191 + FILE *fp;
  192 +
  193 + if(ITER_PRINT==0)return;
  194 +
  195 + fp = fopen(fname, "w");
  196 +
  197 + // fprintf ( fp, "%d\n", M );
  198 + // fprintf ( fp, "%d\n", N );
  199 +
  200 + for (ix = 0 ; ix < nx; ix++) {
  201 + for (iy =0; iy < ny; iy++) {
  202 +
  203 + fprintf(fp, "%6.2f ", u[ix*ny+iy]);
  204 + }
  205 + fputc ( '\n', fp);
  206 + }
  207 +
  208 + printf (" Data written to the output file %s\n", fname);
  209 + fclose(fp);
  210 +}
  211 +
  212 +
  213 +