Changes between Initial Version and Version 1 of cypress/Programming/HeatMassTransfer/ex20.c


Ignore:
Timestamp:
May 14, 2015 3:22:21 PM (7 years ago)
Author:
cmaggio
Comment:

migrated content from ccs wiki

Legend:

Unmodified
Added
Removed
Modified
  • cypress/Programming/HeatMassTransfer/ex20.c

    v1 v1  
     1= Ex20.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  2D Heat/Mass Transfer
     6  Ex20.c : Numerical Analysis
     7  Solve for
     8  d2c/dx2 + d2c/dy2 + 1 = 0
     9  With the boundary conditions of c = 0
     10  along lines of x=1,-1, and y=1, and -1
     11 
     12 */
     13#include<stdio.h>
     14#include<stdlib.h>
     15#include<math.h>
     16#include"cg.h"
     17 
     18#define NUM 100  /* number of points for x and y*/
     19#define TOL 1.0e-12  /* convergence tolerance */
     20 
     21int index(int ii,int jj);
     22int main(void ){
     23  int i,j,k,n;
     24  double x,y,dx;
     25  double phi[NUM][NUM]; /* define 1D array of length NUM * NUM */
     26  double rhs[NUM * NUM]; /* define 1D array of length NUM * NUM */
     27  double sol[NUM * NUM]; /* define 1D array of length NUM * NUM */
     28  double nzval[5 * NUM * NUM]; /* non-zero elements of A */
     29  int colind[5 * NUM * NUM]; /* column indices of nonzeros */
     30  int rowptr[NUM * NUM + 1]; /**/
     31  int nnz;
     32  FILE *fp;
     33 
     34  /*assuming dx = dy : domain is 2x2 */
     35  dx = 2.0 / (double)(NUM - 1);
     36 
     37  /* initial & boundary conditions */
     38  for (j = 0 ; j < NUM  ; j++){
     39    for (i = 0 ; i < NUM ; i++){
     40      phi[i][j] = 0.0;
     41    }
     42  }
     43 
     44  /* Setup Matrix A and RHS */
     45  nnz = 0;
     46  for (j = 1 ; j < NUM - 1 ; j++){
     47    for (i = 1 ; i < NUM - 1 ; i++){
     48      /* initial guess */
     49      sol[index(i,j)] = phi[i][j];
     50 
     51      /* set general RHS */
     52      rhs[index(i,j)] = -dx * dx;
     53 
     54      /* set first nonzero column of row index(i,j) */
     55      rowptr[index(i,j)] = nnz;
     56 
     57      /* sounth */
     58      if (j == 1){
     59        rhs[index(i,j)] -= phi[i][0];
     60      }else{
     61        nzval[nnz] = 1.0;
     62        colind[nnz] = index(i,j-1);
     63        nnz++;
     64      }
     65 
     66      /* west */
     67      if (i == 1){
     68        rhs[index(i,j)] -= phi[0][j];
     69      }else{
     70        nzval[nnz] = 1.0;
     71        colind[nnz] = index(i-1,j);
     72        nnz++;
     73      }
     74 
     75      /* diagonal Element */
     76      nzval[nnz] = -4.0;
     77      colind[nnz] = index(i,j);
     78      nnz++;
     79 
     80      /* east */
     81      if (i == NUM - 2){
     82        rhs[index(i,j)] -= phi[NUM-1][j];
     83      }else{
     84        nzval[nnz] = 1.0;
     85        colind[nnz] = index(i+1,j);
     86        nnz++;
     87      }
     88 
     89      /* north */
     90      if (j == NUM - 2){
     91        rhs[index(i,j)] -= phi[i][NUM - 1];
     92      }else{
     93        nzval[nnz] = 1.0;
     94        colind[nnz] = index(i,j+1);
     95        nnz++;
     96      }
     97    }
     98  }
     99 
     100  /* last element of rowptr si nnz */
     101  rowptr[(NUM-2) * (NUM-2)] = nnz;
     102 
     103  /* solve with CG */
     104  k = solve_cg((NUM-2) * (NUM-2), nnz, nzval, colind, rowptr, sol, rhs, TOL);
     105  if (k < 0){
     106    printf("calculation failed\n");
     107    exit(-1);
     108  }
     109  printf("Number of Interation=%d\n",k);
     110 
     111  /* store solution to main aray */
     112  for (j = 1 ; j < NUM - 1 ; j++){
     113    for (i = 1 ; i < NUM - 1 ; i++){
     114      phi[i][j] = sol[index(i,j)];
     115    }
     116  }
     117 
     118  /* Output Result */
     119  fp = fopen("res.dat","w");
     120  if (fp == NULL){
     121    printf("File could not create\n");
     122    exit(-1);
     123  }
     124 
     125  for (i = 0 ;i < NUM ; i++){
     126    for (j = 0 ; j < NUM ; j++){
     127      x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
     128      y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
     129      fprintf(fp,"%e %e %e\n",x,y,phi[i][j]);
     130    }
     131  }
     132  fclose(fp);
     133 
     134  printf("Done\n");
     135}
     136 
     137/* index for interior Points */
     138int index(int ii,int jj){
     139  return ii - 1 + (NUM - 2) * (jj - 1);
     140}
     141}}}
     142
     143[[cypress/Programming/Cg|cg.c]]