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


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

Legend:

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

    v1 v1  
     1= Ex21.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  1D Heat/Mass Transfer
     6  ex21.c : use cg.c
     7 */
     8#include<stdio.h>
     9#include<stdlib.h>
     10#include<math.h>
     11#include"cg.h"
     12 
     13int main(int argc, char **argv){
     14  double a,b,d;
     15  double *phi;
     16  double *rhs;
     17  double dx;
     18  double tol = 1.0e-12;
     19  int i,num,itc;
     20  double *nzval; /* non-zero elements of A */
     21  int *colind; /* column indices of nonzeros */
     22  int *rowptr; /* row info */
     23  int nnz;
     24  FILE *fp;
     25 
     26  if (argc < 5){
     27    printf("Usage:%s [NUM] [A] [B] [D]\n",argv[0]);
     28    exit(-1);
     29  }
     30 
     31  /* set parameters */
     32  num = atoi(argv[1]);
     33  a = atof(argv[2]);
     34  b = atof(argv[3]);
     35  d = atof(argv[4]);
     36 
     37  printf("num=%d A=%e B=%e D=%e\n",num,a,b,d);
     38 
     39  /* Memory Allocation and Check*/
     40  phi = (double *)calloc(num, sizeof(double));
     41  rhs = (double *)calloc(num, sizeof(double));
     42  nzval = (double *)calloc(3 * num, sizeof(double));
     43  colind = (int *)calloc(3 * num, sizeof(int));
     44  rowptr = (int *)calloc(num + 1, sizeof(int));
     45  if ((phi == NULL) || (rhs == NULL) ||
     46      (nzval == NULL) || (colind == NULL) || (rowptr == NULL)){
     47    printf("Memory Allocation Failed\n");
     48    exit(-1);   
     49  } 
     50 
     51  /** Setup **/
     52  dx = 1.0 / (double)(num - 1);
     53 
     54  /* Boundary Condition*/
     55  phi[0] = a;
     56  phi[num-1] = b;
     57 
     58  /* Setup Matrix A and RHS */
     59  nnz = 0;
     60  for (i = 1 ; i < num - 1 ; i++){
     61 
     62    /* set general RHS */
     63    rhs[i] = -dx * dx / d * (dx * (double)i);
     64 
     65    /* set first nonzero column of row index(i,j) */
     66    rowptr[i - 1] = nnz;
     67 
     68    /* west */
     69    if (i == 1){
     70      rhs[i] -= phi[0];
     71    }else{
     72      nzval[nnz] = 1.0;
     73      colind[nnz] = i-2;
     74      nnz++;
     75    }
     76 
     77    /* diagonal Element */
     78    nzval[nnz] = -2.0;
     79    colind[nnz] = i - 1;
     80    nnz++;
     81 
     82    /* east */
     83    if (i == num - 2){
     84      rhs[i] -= phi[num-1];
     85    }else{
     86      nzval[nnz] = 1.0;
     87      colind[nnz] = i;
     88      nnz++;
     89    }
     90  }
     91 
     92  /* last element of rowptr si nnz */
     93  rowptr[num-2] = nnz;
     94 
     95  /* solve with CG */
     96  itc = solve_cg(num-2, nnz, nzval, colind, rowptr, &phi[1], &rhs[1], tol);
     97  if (itc < 0){
     98    printf("calculation failed\n");
     99    exit(-1);
     100  }
     101 
     102  printf("Number of Iteration=%d\n",itc);
     103 
     104  /* Output Result */
     105  fp = fopen("res.dat","w");
     106  if (fp == NULL){
     107    printf("File could not create\n");
     108    exit(-1);
     109  }
     110 
     111  for (i = 0 ; i < num ; i++){
     112    fprintf(fp,"%e %e\n",dx * (double)i,phi[i]);
     113  }
     114  fclose(fp);
     115 
     116  /* END : Deallocate Memory */
     117  free(phi);
     118  free(rhs);
     119  free(nzval);
     120  free(colind);
     121  free(rowptr);
     122 
     123  printf("Done\n");
     124}
     125}}}
     126
     127[[cypress/Programming/Cg|cg.c]]