solver.c

cfd求解器使用与gmsh网格的求解

#define RCSID "$Id: Solver.c,v 1.28 2006/02/26 00:43:00 geuzaine Exp $"
/*
 * Copyright (C) 1997-2006 P. Dular, C. Geuzaine
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 *
 * Please report all bugs and problems to <[e-mail]>.
 *
 * Contributor(s):
 *   Jean-Francois Remacle
 *   Benoit Meys
 *   Johan Gyselinck
 *   Ruth Sabariego
 */

#include <stdio.h>
#include <string.h>
#include <math.h>

#include "Solver.h"
#include "Solver_F.h"
#include "Message.h"
#include "Malloc.h"
#include "F_FMMOperations.h"


/* ------------------------------------------------------------------------ */
/*  S o l v e _ m a t r i x                                                 */
/* ------------------------------------------------------------------------ */


extern int Flag_FMM ;
extern int Flag_DTA ;

void solve_matrix (Matrix *M, Solver_Params *p, double *b, double *x){
  FILE    *pf;
  double   fpar[17];
  double  *a, *w, *rhs, *sol, *dx ;
  double   res;
  int      i, j, k, nnz, nnz_ilu, ierr, ipar[17];
  int     *ja, *ia, *jw, *mask, *levels;
  int      its, end, do_permute=0;
  int      zero=0, un=1, deux=2, six=6, douze=12, trente=30, trente_et_un=31;
  int      ROW=0, COLUMN=1;
  double   res1=1.;
  int      TrueNnz=0;

  if (!M->N) {
    Msg(WARNING, "No equations in linear system");
    return;
  }
  
  for(i=0 ; i<M->N ; i++){
    if(b[i] != 0.) break;
    if(i == M->N-1) {
      Msg(WARNING, "Null right hand side in linear system");
      /*
	for(i=0 ; i<M->N ; i++)  x[i] = 0. ;    
	return ;
      */
    }
  }

  if(M->T == DENSE){

    if(p->Algorithm == LU){

      Msg(SPARSKIT, "Dense LU decomposition\n");
      print_matrix_info_DENSE(M->N);

      sol = (double*)Malloc(M->N * sizeof(double));
      w   = (double*)Malloc(M->N * sizeof(double));
      dx  = (double*)Malloc(M->N * sizeof(double));
      
      ipar[1] = p->Re_Use_LU;
      ipar[2] = p->Iterative_Improvement;
      ipar[3] = p->Matrix_Printing;
      ipar[4] = p->Nb_Iter_Max;
      
      fpar[1] = p->Stopping_Test;
      
      flu_(&ipar[1], &fpar[1], M->F.a, M->F.lu, &M->N, &M->N, b, x, dx, sol, w);
    
      Free(sol);
      Free(w);
      Free(dx);

      return ;
    }

    Msg(SPARSKIT, "Dense to sparse matrix conversion\n") ;

    nnz = M->N * M->N ;

    M->S.a    = List_Create(1, 1, sizeof(double));
    M->S.a->n = nnz ;
    M->S.a->array = (char*) M->F.a ;

    M->S.jptr = List_Create(M->N+1, M->N, sizeof(int));
    M->S.ai   = List_Create(nnz, M->N, sizeof(int));

    for(i=1 ; i<=nnz ; i+=M->N){ 
      List_Add(M->S.jptr, &i) ;
      for(j=1 ; j<=M->N ; j++){
	List_Add(M->S.ai, &j) ;
      }
    }

    i = nnz + 1 ; List_Add(M->S.jptr, &i) ;

    if(M->changed){
      do_permute = 1 ;
      M->changed = 0 ;
    }

    for (i=0 ; i<nnz ; i++) if (M->F.a[i]) TrueNnz++ ;
    Msg (INFO, "Number of nonzeros %d/%d (%.4f)",TrueNnz, nnz, (double)TrueNnz/(double)nnz);

    Msg(RESOURCES, "");

  } /* if DENSE */
  else{

    nnz = List_Nbr(M->S.a);    
    if(M->changed){
      do_permute = 1 ;
      csr_format (&M->S, M->N);
      restore_format (&M->S);
      M->changed = 0 ;
    }

  } /* if SPARSE */

  a  = (double*) M->S.a->array;
  ia = (int*) M->S.jptr->array;
  ja = (int*) M->S.ai->array; 

  if(p->Scaling != NONE){
    Msg(SPARSKIT, "Scaling system of equations\n") ;
    scale_matrix (p->Scaling, M) ;
    scale_vector (ROW, M, b) ;    
  }
  else{
    Msg(SPARSKIT, "No scaling of system of equations\n") ;
  }

  for (i=1; i<M->N; i++) {
    if(ia[i]-ia[i-1] <= 0)
      Msg(GERROR, "Zero row in matrix");
  }

  rhs = (double*) Malloc(M->N * sizeof(double));
  sol = (double*) Calloc(M->N, sizeof(double));
  
  /* Renumbering */
  
  if (!M->ILU_Exists){
    M->S.permr  = (int*) Malloc(M->N * sizeof(int));
    M->S.rpermr = (int*) Malloc(M->N * sizeof(int));
    M->S.permp  = (int*) Malloc(2 * M->N * sizeof(int));
  }

  if(do_permute || !M->ILU_Exists){
    for(i=0 ; i<M->N ; i++) {
      M->S.permr[i] = M->S.rpermr[i] = M->S.permp[i+M->N] = i+1;
    } 
    switch (p->Renumbering_Technique){ 	
    case NONE:
      Msg(SPARSKIT, "No renumbering\n");
      break;	
    case RCMK: 
      Msg(SPARSKIT, "RCMK algebraic renumbering\n");	
      if(!M->ILU_Exists){
	M->S.a_rcmk  = (double*) Malloc(nnz * sizeof(double));
	M->S.ia_rcmk = (int*) Malloc((M->N + 1) * sizeof(int));
	M->S.ja_rcmk = (int*) Malloc(nnz * sizeof(int));
      }
      mask    = (int*) Malloc(nnz * sizeof(int));
      levels  = (int*) Malloc((M->N + 1) * sizeof(int));      
      i = j = k = 1;
      cmkreord_(&M->N, a, ja, ia, M->S.a_rcmk, M->S.ja_rcmk, M->S.ia_rcmk, 
		&i, M->S.permr, mask, &j, &k, M->S.rpermr, levels);      
      w = (double*) Malloc(nnz * sizeof(double));      
      sortcol_(&M->N, M->S.a_rcmk, M->S.ja_rcmk, M->S.ia_rcmk, mask, w);      
      Free(w); Free(mask); Free(levels);
      break; 	
    default :
      Msg(GERROR, "Unknown renumbering technique");
      break;
    }
    print_matrix_info_CSR(M->N, ia, ja);
    Msg(RESOURCES, "");
  }
  
  if(p->Renumbering_Technique == RCMK){
    if (p->Re_Use_ILU && !M->ILU_Exists && !do_permute){ 
      /*
      This is incorrect if M is to be changed during the process, and
      we still want to keep the same precond. 

      Free(M->S.a->array) ;
      Free(M->S.jptr->array) ;
      Free(M->S.ai->array) ;
      M->S.a->array =  (char*)M->S.a_rcmk;
      M->S.jptr->array = (char*)M->S.ia_rcmk;
      M->S.ai->array = (char*)M->S.ja_rcmk; 
      */
    }
    a = M->S.a_rcmk;
    ia = M->S.ia_rcmk;
    ja = M->S.ja_rcmk;
  }
  

  if (p->Matrix_Printing == 1 || p->Matrix_Printing == 3) {
    Msg(SPARSKIT, "Matrix printing\n");
    skit_(&M->N, a, ja, ia, &douze, &douze, &ierr); 
    pf = fopen("fort.13","w");
    for (i=0 ; i<M->N ; i++) fprintf(pf, "%d %22.15E\n", i+1, b[i]);
    fclose(pf);
    psplot_(&M->N, ja, ia, &trente, &zero);
  }
  
  /* Incomplete factorizations */
  
  if (!M->ILU_Exists) {
    
    if (p->Re_Use_ILU) M->ILU_Exists = 1;
    
#if defined(HAVE_ILU_FLOAT)
#define ILUSTORAGE "Float"
#else
#define ILUSTORAGE "Double"
#endif
  
    end = 0 ;
    
    switch (p->Preconditioner){
      
    case ILUT  : Msg(SPARSKIT, "ILUT (%s, fill-in = %d)\n", ILUSTORAGE, p->Nb_Fill);
      nnz_ilu = 2 * (M->N+1) * (p->Nb_Fill+1); break;
      
    case ILUTP : Msg(SPARSKIT, "ILUTP (%s, fill-in = %d)\n", ILUSTORAGE, p->Nb_Fill);
      nnz_ilu = 2 * (M->N+1) * (p->Nb_Fill+1); break;
      
    case ILUD  : Msg(SPARSKIT, "ILUD (%s)\n", ILUSTORAGE);    
      /* first guess */
      nnz_ilu = List_Nbr(M->S.a); break;
      
    case ILUDP : Msg(SPARSKIT, "ILUDP (%s)\n", ILUSTORAGE);
      /* first guess */
      nnz_ilu = List_Nbr(M->S.a); break ;
      
    case ILUK  : Msg(SPARSKIT, "ILU%d (%s)\n", p->Nb_Fill, ILUSTORAGE);    
      /* exact for nbfill=0, first guess otherwise */
      nnz_ilu = (p->Nb_Fill+1) * List_Nbr(M->S.a) + (M->N+1);
      break;
      
    case ILU0  : Msg(SPARSKIT, "ILU0 (%s)\n", ILUSTORAGE);    
      nnz_ilu = List_Nbr(M->S.a) + (M->N+1); break;
      
    case MILU0 : Msg(SPARSKIT, "MILU0 (%s)\n", ILUSTORAGE);    
      nnz_ilu = List_Nbr(M->S.a) + (M->N+1); break;
      
    case DIAGONAL : 
      Msg(SPARSKIT, "Diagonal scaling (%s)\n", ILUSTORAGE);
      M->S.alu = (scalar*) Malloc((M->N+1) * sizeof(scalar));
      M->S.jlu = (int*) Malloc((M->N+1) * sizeof(int));
      M->S.ju  = (int*) Malloc((M->N+1) * sizeof(int));
      
      for (i=0 ; i<M->N ; i++) {
	M->S.alu[i] = 1.0 ;
	M->S.jlu[i] = M->N+2 ;  
	M->S.ju[i]  = M->N+2 ;
      }
      M->S.alu[M->N] = 0.0 ;
      M->S.jlu[M->N] = M->N+2 ;
      M->S.ju[M->N]  = 0 ;
      end = 1;
      ierr = 0;
      break;
      
    case NONE  : 
      Msg(SPARSKIT, "No ILU\n");
      end = 1;
      ierr = 0;
      break ;

    default :
      Msg(GERROR, "Unknown ILU method");
      break;
    }
    
    if(!end){
      M->S.alu = (scalar*) Malloc(nnz_ilu * sizeof(scalar));
      M->S.jlu = (int*) Malloc(nnz_ilu * sizeof(int));
      M->S.ju  = (int*) Malloc((M->N+1) * sizeof(int));
    }

    reallocate :
    
    switch(p->Preconditioner){
    case ILUT :