linalg_petsc.c

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

#define RCSID "$Id: LinAlg_PETSC.c,v 1.63 2006/02/26 00:42:52 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):
 *   David Colignon
 *   Ruth Sabariego
 *   Jose Geraldo A. Brito Neto
 */

/* 
   This is the interface library for the PETSC solvers

   Options for PETSc can be provided on the command line, or in the file
   ~/.petscrc. 

   By default, we use the following options (GMRES iterative solver,
   1.e-8 relative tolerance with ILU(6) preconditioner and RCMK
   renumbering):

   -pc_type ilu
   -pc_ilu_levels 6
   -pc_ilu_mat_ordering_type rcm
   -ksp_rtol 1.e-8

   Other useful options include:

   -ksp_gmres_restart 100
   -ksp_monitor
   ...

   To use a direct solver (a sparse lu) instead of an iterative
   solver, use

   -ksp_type preonly -pc_type lu

   When PETSc is compiled with external solvers (UMFPACK, SuperLU,
   etc.), they can be accessed simply by changing the matrix type. For
   example, with umfpack:

   -mat_type umfpack -ksp_type preonly -pc_type lu 
*/

#include "GetDP.h"

#if defined(HAVE_PETSC)

#include "LinAlg.h"
#include "F_FMMOperations.h"
#include "SafeIO.h"
#include "CurrentData.h"

extern int Flag_FMM;

/* error checking, petsc-style */

#define MYCHECK(ierr) CHKERRABORT(PETSC_COMM_WORLD,ierr)
static int  ierr, SolverInitialized = 0;

/* stuff for matrix-free matrices */

typedef struct{
  void *ctx;
  void (*fct)(gMatrix *A, gVector *x, gVector *y);
} petscCtx;

int petscMult(Mat A, Vec x, Vec y){
  gMatrix A1;
  gVector x1, y1;
  petscCtx *ctx;

  A1.M = A; x1.V = x; y1.V = y;
  MatShellGetContext(A,(void**)&ctx);
  ctx->fct(&A1,&x1,&y1);
  return 0;
}

/* Initialize */

void LinAlg_Initialize(int* argc, char*** argv, int *NbrCpu, int *RankCpu){

  GetDP_Begin("LinAlg_Initialize");

  MPI_Init(argc, argv);
  MPI_Comm_size(MPI_COMM_WORLD, NbrCpu);
  MPI_Comm_rank(MPI_COMM_WORLD, RankCpu);

  GetDP_End;
}

void LinAlg_InitializeSolver(int* argc, char*** argv, int *NbrCpu, int *RankCpu){

  GetDP_Begin("LinAlg_InitializeSolver");

  /* This function detects if MPI is initialized */
  PetscInitialize(argc, argv, PETSC_NULL, PETSC_NULL);
  SolverInitialized = 1;
  /* just in case... */
  MPI_Comm_size(PETSC_COMM_WORLD, NbrCpu);
  MPI_Comm_rank(PETSC_COMM_WORLD, RankCpu);

  GetDP_End;
}

/* Finalize */

void LinAlg_Finalize(void){
  int flag;

  GetDP_Begin("LinAlg_Finalize");

  MPI_Initialized(&flag);
  if(flag) MPI_Finalize();

  GetDP_End;
}

void LinAlg_FinalizeSolver(void){

  GetDP_Begin("LinAlg_Finalize");

  if(SolverInitialized) PetscFinalize();

  GetDP_End;
}

/* Barrier */

void LinAlg_Barrier(void){

  GetDP_Begin("LinAlg_Barrier");

  MPI_Barrier(PETSC_COMM_WORLD);

  GetDP_End;
}

/* Sequential */

void LinAlg_SequentialBegin(void){

  GetDP_Begin("LinAlg_SequentialBegin");

  ierr = PetscSequentialPhaseBegin(MPI_COMM_WORLD, 1); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_SequentialEnd(void){

  GetDP_Begin("LinAlg_SequentialEnd");

  ierr = PetscSequentialPhaseEnd(MPI_COMM_WORLD, 1); MYCHECK(ierr);

  GetDP_End;
}

/* Create */

void LinAlg_CreateSolver(gSolver *Solver, char * SolverDataFileName){
  Solver->ksp = NULL;
}

void LinAlg_CreateVector(gVector *V, gSolver *Solver, int n, int NbrPart, int *Part){
#if defined(HAVE_METIS)
  int NbrCpu, RankCpu;
#endif

  GetDP_Begin("LinAlg_CreateVector");

#if defined(HAVE_METIS)
  MPI_Comm_size(MPI_COMM_WORLD, &NbrCpu);
  MPI_Comm_rank(PETSC_COMM_WORLD, &RankCpu);
  if(NbrPart != NbrCpu){
    Msg(WARNING, "%d partitions on %d CPU", NbrPart, NbrCpu);
    ierr = VecCreate(PETSC_COMM_WORLD, &V->V); MYCHECK(ierr);
    ierr = VecSetSizes(V->V,PETSC_DECIDE,n); MYCHECK(ierr);
    ierr = VecSetFromOptions(V->V); MYCHECK(ierr);
  }
  else{
    ierr = VecCreateMPI(PETSC_COMM_WORLD, Part[RankCpu+1]-Part[RankCpu], 
			n, &V->V); MYCHECK(ierr);
  }
#else
  ierr = VecCreate(PETSC_COMM_WORLD, &V->V); MYCHECK(ierr);
  ierr = VecSetSizes(V->V, PETSC_DECIDE, n); MYCHECK(ierr);

  /* set some default options */
  ierr = VecSetType(V->V, VECSEQ); MYCHECK(ierr);
  
  /* override the default options with the ones from the option
     database (if any) */
  ierr = VecSetFromOptions(V->V); MYCHECK(ierr);
#endif

  GetDP_End;
}

void LinAlg_CreateMatrix(gMatrix *M, gSolver *Solver, int n, int m, 
			 int NbrPart, int *Part, int *Nnz){
#if defined(HAVE_METIS)
  int NbrCpu, RankCpu, i_Start, i_End;
#endif

  GetDP_Begin("LinAlg_CreateMatrix");

#if defined(HAVE_METIS)
  MPI_Comm_size(PETSC_COMM_WORLD, &NbrCpu);
  MPI_Comm_rank(PETSC_COMM_WORLD, &RankCpu);

  if(NbrPart != NbrCpu) Msg(GERROR, "%d partitions on %d CPU", NbrPart, NbrCpu);

  i_Start = Part[RankCpu]-1;
  i_End   = Part[RankCpu+1]-1;
  ierr = MatCreateMPIAIJ(PETSC_COMM_WORLD,
			 i_End-i_Start, i_End-i_Start, 
			 n, m, 
			 1, &Nnz[i_Start], 
			 1, &Nnz[i_Start], 
			 &M->M); MYCHECK(ierr); 
#else
  ierr = MatCreate(PETSC_COMM_WORLD, &M->M); MYCHECK(ierr); 
  ierr = MatSetSizes(M->M, PETSC_DECIDE, PETSC_DECIDE, n, m); MYCHECK(ierr); 

  /* set some default options */
  ierr = MatSetType(M->M, MATSEQAIJ); MYCHECK(ierr); 
  ierr = MatSeqAIJSetPreallocation(M->M, 50, PETSC_NULL); MYCHECK(ierr); 

  /* override the default options with the ones from the option
     database (if any) */
  ierr = MatSetFromOptions(M->M); MYCHECK(ierr);
#endif

  GetDP_End;
}

void LinAlg_CreateMatrixShell(gMatrix *A, gSolver *Solver, int n, int m, void *myCtx, 
			      void (*myMult)(gMatrix *A, gVector *x, gVector *y)){
  petscCtx *ctx;

  GetDP_Begin("LinAlg_CreateMatrixShell");

  ctx = (petscCtx*)Malloc(sizeof(petscCtx)); /* memory leak! to change... */
  ctx->ctx = myCtx;
  ctx->fct = myMult;
  MatCreateShell(PETSC_COMM_WORLD,n,m,PETSC_DETERMINE,
		 PETSC_DETERMINE,(void*)ctx,&A->M);
  MatShellSetOperation(A->M,MATOP_MULT,(void(*)())petscMult);

  GetDP_End;
}

/* Destroy */

void LinAlg_DestroySolver(gSolver *Solver){

  GetDP_Begin("LinAlg_DestroySolver");

  ierr = KSPDestroy(Solver->ksp); MYCHECK(ierr);
  Solver->ksp = NULL;

  GetDP_End;
}

void LinAlg_DestroyVector(gVector *V){

  GetDP_Begin("LinAlg_DestroyVector");

  ierr = VecDestroy(V->V); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_DestroyMatrix(gMatrix *M){

  GetDP_Begin("LinAlg_DestroyMatrix");

  ierr = MatDestroy(M->M); MYCHECK(ierr);

  GetDP_End;
}

/* Copy */

void LinAlg_CopyScalar(gScalar *S1, gScalar *S2){

  GetDP_Begin("LinAlg_CopyScalar");

  S1->s = S2->s;

  GetDP_End;
}

void LinAlg_CopyVector(gVector *V1, gVector *V2){

  GetDP_Begin("LinAlg_CopyVector");

  ierr = VecCopy(V1->V, V2->V); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_CopyMatrix(gMatrix *M1, gMatrix *M2){

  GetDP_Begin("LinAlg_CopyMatrix");

  ierr = MatCopy(M1->M, M2->M, DIFFERENT_NONZERO_PATTERN); MYCHECK(ierr);

  GetDP_End;
}

/* Zero */

void LinAlg_ZeroScalar(gScalar *S){

  GetDP_Begin("LinAlg_ZeroScalar");

  S->s = 0.;

  GetDP_End;
}

void LinAlg_ZeroVector(gVector *V){
  PetscScalar zero = 0.0;

  GetDP_Begin("LinAlg_ZeroVector");

  ierr = VecSet(V->V, zero); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_ZeroMatrix(gMatrix *M){

  GetDP_Begin("LinAlg_ZeroMatrix");

  ierr = MatZeroEntries(M->M); MYCHECK(ierr);

  GetDP_End;
}

/* Scan */

void LinAlg_ScanScalar(FILE *file, gScalar *S){
#if PETSC_USE_COMPLEX
  double a, b;
#endif

  GetDP_Begin("LinAlg_ScanScalar");

#if PETSC_USE_COMPLEX
  fscanf(file, "%lf %lf", &a, &b);
  S->s = a + PETSC_i * b;
#else
  fscanf(file, "%lf", &S->s);
#endif

  GetDP_End;
}

void LinAlg_ScanVector(FILE *file, gVector *V) {
  int i, n;
  PetscScalar tmp;
  double a, b;  

  GetDP_Begin("LinAlg_ScanVector");
  
  ierr = VecGetSize(V->V, &n); MYCHECK(ierr);
  for(i=0; i<n; i++){
#if PETSC_USE_COMPLEX
    fscanf(file, "%lf %lf", &a, &b);
    tmp = a + PETSC_i * b;
#else
    fscanf(file, "%lf", &a);
    tmp = a;
#endif
    ierr = VecSetValues(V->V, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);
  }

  GetDP_End;
}

void LinAlg_ScanMatrix(FILE *file, gMatrix *M){

  GetDP_Begin("LinAlg_ScanMatrix");

  Msg(GERROR, "'LinAlg_ScanMatrix' not yet implemented");  

  GetDP_End;
}

/* Read */

void LinAlg_ReadScalar(FILE *file, gScalar *S){

  GetDP_Begin("LinAlg_ReadScalar");

  Msg(GERROR, "'LinAlg_ReadScalar' not yet implemented");

  GetDP_End;
}

void LinAlg_ReadVector(FILE *file, gVector *V) {
  int i, n;
  PetscScalar *tmp;

  GetDP_Begin("LinAlg_ReadVector");

  ierr = VecGetSize(V->V, &n); MYCHECK(ierr);
  tmp = (PetscScalar*)Malloc(n*sizeof(PetscScalar));
  fread(tmp, sizeof(PetscScalar), n, file);

  for(i=0; i<n; i++){
    ierr = VecSetValues(V->V, 1, &i, &tmp[i], INSERT_VALUES); MYCHECK(ierr);
  }

  Free(tmp);

  GetDP_End;
}

void LinAlg_ReadMatrix(FILE *file, gMatrix *M){

  GetDP_Begin("LinAlg_ReadMatrix");

  Msg(GERROR, "'LinAlg_ReadMatrix' not yet implemented");  

  GetDP_End;
}

/* Print */

void LinAlg_PrintScalar(FILE *file, gScalar *S){

  GetDP_Begin("LinAlg_PrintScalar");

#if PETSC_USE_COMPLEX
  fprintf(file, "%.16g %.16g", real(S->s), imag(S->s));
#else
  fprintf(file, "%.16g", S->s);
#endif

  GetDP_End;
}

void LinAlg_PrintVector(FILE *file, gVector *V){
  PetscScalar *tmp;
  int     i, n;

  GetDP_Begin("LinAlg_PrintVector");

  ierr = VecGetLocalSize(V->V, &n); MYCHECK(ierr);
  ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);

  for (i=0; i<n; i++){
#if PETSC_USE_COMPLEX
    fprintf(file, "%.16g %.16g\n", real(tmp[i]), imag(tmp[i]));
#else
    fprintf(file, "%.16g\n", tmp[i]);
#endif
  }
  fflush(file);

  ierr = VecRestoreArray(V->V, &tmp); MYCHECK(ierr);

  GetDP_End;
} 

void LinAlg_PrintMatrix(FILE *file, gMatrix *M){

  GetDP_Begin("LinAlg_PrintMatrix");
  
  ierr = MatView(M->M,PETSC_VIEWER_STDOUT_WORLD); MYCHECK(ierr);

  GetDP_End;
}