linalg_petsc.c

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

  else if(V3 == V2){
    ierr = VecAYPX(V2->V, tmp, V1->V); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_SubVectorVector'");  

  GetDP_End;
}

void LinAlg_SubMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3){

  GetDP_Begin("LinAlg_SubMatrixMatrix");

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

  GetDP_End;
}

/* Prod */

void LinAlg_ProdScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3){

  GetDP_Begin("LinAlg_ProdScalarScalar");

  S3->s = S1->s * S2->s;

  GetDP_End;
}
void LinAlg_ProdScalarDouble(gScalar *S1, double d, gScalar *S2){

  GetDP_Begin("LinAlg_ProdScalarDouble");

  S2->s = S1->s * d;

  GetDP_End;
}

void LinAlg_ProdScalarComplex(gScalar *S, double d1, double d2, double *d3, double *d4){
#if PETSC_USE_COMPLEX
  PetscScalar tmp;
#endif

  GetDP_Begin("LinAlg_ProdScalarComplex");

#if PETSC_USE_COMPLEX
  tmp = S->s * (d1 + PETSC_i * d2);
  *d3 = real(tmp);
  *d4 = imag(tmp);
#else
  *d3 = S->s * d1;
  *d4 = S->s * d2;
#endif

  GetDP_End;
}

void LinAlg_ProdVectorScalar(gVector *V1, gScalar *S, gVector *V2){  

  GetDP_Begin("LinAlg_ProdVectorScalar");

  if(V2 == V1){
    ierr = VecScale(V1->V, S->s); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdVectorScalar'");

  GetDP_End;
}

void LinAlg_ProdVectorDouble(gVector *V1, double d, gVector *V2){
  PetscScalar tmp;
  tmp = d;

  GetDP_Begin("LinAlg_ProdVectorDouble");

  if(V2 == V1){
    ierr = VecScale(V1->V, tmp); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdVectorDouble'");

  GetDP_End;
}

void LinAlg_ProdVectorComplex(gVector *V1, double d1, double d2, gVector *V2){

  GetDP_Begin("LinAlg_ProdvectorComplex");

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

  GetDP_End;
}

void LinAlg_ProdVectorVector(gVector *V1, gVector *V2, double *d){
  PetscScalar tmp;

  GetDP_Begin("LinAlg_ProdVectorVector");

  ierr = VecDot(V1->V, V2->V, &tmp); MYCHECK(ierr);
#if PETSC_USE_COMPLEX
  *d = real(tmp);
#else
  *d = tmp;
#endif

  GetDP_End;
}

void LinAlg_ProdMatrixVector(gMatrix *M, gVector *V1, gVector *V2){

  GetDP_Begin("LinAlg_ProdMatrixVector");

  if(V2==V1)
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdMatrixVector'");
  else
    ierr = MatMult(M->M, V1->V, V2->V); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_ProdMatrixScalar(gMatrix *M1, gScalar *S, gMatrix *M2){

  GetDP_Begin("LinAlg_ProdMatrixScalar");

  if(M2 == M1){
    ierr = MatScale(M1->M, S->s); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdMatrixScalar'");

  GetDP_End;
}

void LinAlg_ProdMatrixDouble(gMatrix *M1, double d, gMatrix *M2){
  PetscScalar tmp;
  tmp = d;

  GetDP_Begin("LinAlg_ProdMatrixDouble");

  if(M2 == M1){
    ierr = MatScale(M1->M, tmp); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdMatrixDouble'");

  GetDP_End;
}

void LinAlg_ProdMatrixComplex(gMatrix *M1, double d1, double d2, gMatrix *M2){

  GetDP_Begin("LinAlg_ProdMatrixComplex");

#if PETSC_USE_COMPLEX
  if(M2 == M1){
    PetscScalar tmp = d1 + (PETSC_i * d2);
    ierr = MatScale(M1->M, tmp); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdMatrixDouble'");
#else
  Msg(GERROR, "'LinAlg_ProdMatrixComplex' not yet implemented");
#endif
  
  GetDP_End;
}

/* Div */

void LinAlg_DivScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3){

  GetDP_Begin("LinAlg_DivScalarScalar");

  S3->s = S1->s / S2->s;

  GetDP_End;
}
void LinAlg_DivScalarDouble(gScalar *S1, double d, gScalar *S2){

  GetDP_Begin("LinAlg_DivScalarDouble");

  S2->s = S1->s / d;

  GetDP_End;
}

/* Norm */

void LinAlg_VectorNorm2(gVector *V1, double *norm){
  PetscReal tmp;

  GetDP_Begin("LinAlg_VectorNorm2");
  
  VecNorm(V1->V, NORM_2, &tmp);

  *norm = tmp;

  GetDP_End;
}

void LinAlg_VectorNormInf(gVector *V1, double *norm){
  PetscReal tmp;

  GetDP_Begin("LinAlg_VectorNormInf");
  
  VecNorm(V1->V, NORM_INFINITY, &tmp);

  *norm = tmp;

  GetDP_End;
}

/* Assemble */

void LinAlg_AssembleMatrix(gMatrix *M){

  GetDP_Begin("LinAlg_AssembleMatrix");

  ierr = MatAssemblyBegin(M->M, MAT_FINAL_ASSEMBLY); MYCHECK(ierr);
  ierr = MatAssemblyEnd(M->M, MAT_FINAL_ASSEMBLY); MYCHECK(ierr);  

  GetDP_End;
}

void LinAlg_AssembleVector(gVector *V){

  GetDP_Begin("LinAlg_AssembleVector");

  ierr = VecAssemblyBegin(V->V); MYCHECK(ierr);
  ierr = VecAssemblyEnd(V->V); MYCHECK(ierr);

  GetDP_End;
}

/* FMM interface routine */

int LinAlg_ApplyFMM(void *ptr, Vec xin, Vec xout){

  static int ii = 0;
  static Vec DTAxi;
  PetscScalar *tmpV, tmp, mone = -1.0;
  double *x = NULL, *y = NULL;
  int i, n;

  GetDP_Begin("LinAlg_ApplyFMM");

  ierr = VecCopy(xin, xout); MYCHECK(ierr);
  ierr = VecDuplicate(xin, &DTAxi); MYCHECK(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD,"iteration %d xin vector:\n",ii); MYCHECK(ierr);
  ierr = VecView(xin,PETSC_VIEWER_STDOUT_WORLD); MYCHECK(ierr);

  ierr = VecGetLocalSize(xin, &n);  MYCHECK(ierr);
  ierr = VecGetArray(xin, &tmpV);  MYCHECK(ierr);
 
  if(!(ii%2)){    
#if PETSC_USE_COMPLEX
    x   = (double*)Malloc(2*n*sizeof(double));
    y   = (double*)Calloc(2*n,sizeof(double));
    for(i = 0; i < n; i++){    
      x[2*i]   = real(tmpV[i]);
      x[2*i+1] = imag(tmpV[i]);
    }
#else
    x   = (double*)Malloc(n*sizeof(double));
    y   = (double*)Calloc(n,sizeof(double));
    for(i = 0; i < n; i++) x[i] = tmpV[i];
#endif
   
    FMM_MatVectorProd( x, y );
    
    for (i = 0; i < n; i++){
#if PETSC_USE_COMPLEX   
      tmp = y[2*i] + PETSC_i* y[2*i+1]; 
#else     
      tmp = y[i];
#endif
      VecSetValues(DTAxi, 1, &i, &tmp, INSERT_VALUES);
    }    
  }
  else{
    ierr = VecAXPY(xout, mone,DTAxi); MYCHECK(ierr);
  }

  ierr = VecRestoreArray(xin, &tmpV); MYCHECK(ierr);
  ierr = VecAssemblyEnd(xin); MYCHECK(ierr);
  ierr = VecAssemblyEnd(xout); MYCHECK(ierr);
  Free(x);
  Free(y);
 
  ii++;
  printf("ITERATIONS! %d %d\n",ii,n);
  GetDP_Return(0);
}

void LinAlg_AddFMMDTAx( Vec xin, Vec xout ){

  PetscScalar *V, tmp;
  double *x, *y;
  int i, n;

  GetDP_Begin("LinAlg_AddFMMDTAx");

  ierr = VecGetLocalSize(xin, &n);  MYCHECK(ierr);
  ierr = VecGetArray(xin, &V);  MYCHECK(ierr);
  ierr = VecRestoreArray(xin, &V); MYCHECK(ierr);

#if PETSC_USE_COMPLEX
  x   = (double*)Malloc(2*n*sizeof(double));
  y   = (double*)Calloc(2*n,sizeof(double));
  for(i = 0; i < n; i++){    
    x[2*i]   = real(V[i]);
    x[2*i+1] = imag(V[i]);
  }
#else
  x   = (double*)Malloc(n*sizeof(double));
  y   = (double*)Calloc(n,sizeof(double));
  for(i = 0; i < n; i++) x[i] = V[i];
#endif
  
  FMM_MatVectorProd( x, y ); /* y contains the DTAxi values */
  
  for (i = 0; i < n; i++){
#if PETSC_USE_COMPLEX   
    tmp = y[2*i] + PETSC_i* y[2*i+1]; 
#else     
    tmp = y[i];
#endif
    ierr = VecSetValues(xout, 1, &i, &tmp, ADD_VALUES); MYCHECK(ierr);
  }    

  ierr = VecRestoreArray(xin, &V); MYCHECK(ierr);
  ierr = VecAssemblyEnd(xin); MYCHECK(ierr);
  ierr = VecAssemblyEnd(xout); MYCHECK(ierr);
  Free(x);
  Free(y);
 
  GetDP_End;
}

int LinAlg_ApplyFMMMonitor(KSP ksp, int it,double rnorm,void *dummy){

  Vec      x, rhs, DTAxi;
  /* Vec      Residual, t, v; */
  PetscScalar *tmpV, tmp, mone = -1.0, zero = 0.0;
  double *x_it, *y_it;
  int i, n;

  GetDP_Begin("LinAlg_ApplyFMMMonitor");
  
  /* Build the solution vector */
  
  ierr = KSPBuildSolution(ksp,PETSC_NULL, &x); MYCHECK(ierr);
  
  ierr = VecDuplicate(x, &DTAxi); MYCHECK(ierr);
  ierr = VecSet(DTAxi, zero); MYCHECK(ierr);
  
  ierr = VecGetLocalSize(x, &n);  MYCHECK(ierr);
  ierr = VecGetArray(x, &tmpV);  MYCHECK(ierr);
  
#if PETSC_USE_COMPLEX
  x_it   = (double*)Malloc(2*n*sizeof(double));
  y_it   = (double*)Calloc(2*n,sizeof(double));
  for(i = 0; i < n; i++){    
    x_it[2*i]   = real(tmpV[i]);
    x_it[2*i+1] = imag(tmpV[i]);
  }
#else
  x_it = (double*)Malloc(n*sizeof(double));
  y_it   = (double*)Calloc(n, sizeof(double));
  for(i = 0; i < n; i++) x_it[i] = tmpV[i];
#endif
  ierr = VecRestoreArray(x, &tmpV); MYCHECK(ierr);

  FMM_MatVectorProd( x_it, y_it );

  for(i = 0; i < n; i++){
#if PETSC_USE_COMPLEX   
    tmp = y_it[2*i] + PETSC_i* y_it[2*i+1]; 
#else     
    tmp = y_it[i];
#endif
    ierr = VecSetValues(DTAxi, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);
  }

  ierr = KSPGetRhs(ksp, &rhs); MYCHECK(ierr);

  ierr = VecAYPX(rhs, mone, DTAxi); MYCHECK(ierr);

  /* ierr = VecAssemblyEnd(rhs); MYCHECK(ierr); */
  ierr = VecDestroy(DTAxi); MYCHECK(ierr);

  /*
  ierr = PetscPrintf(PETSC_COMM_WORLD,"iteration %d rhs vector:\n",it); MYCHECK(ierr);
  ierr = VecView(rhs,PETSC_VIEWER_STDOUT_WORLD); MYCHECK(ierr);
  ierr = VecDuplicate(x, &t); MYCHECK(ierr);
  ierr = VecSet(&zero,t); MYCHECK(ierr);
  ierr = VecDuplicate(x, &v); MYCHECK(ierr);
  ierr = VecSet(&zero,v); MYCHECK(ierr);
  ierr = KSPBuildResidual(ksp, t, v, &Residual); MYCHECK(ierr);
  */

  ierr = PetscPrintf(PETSC_COMM_WORLD,"iteration %d KSP Residual norm %14.12e \n",
		     it, rnorm); MYCHECK(ierr);

  GetDP_Return(0);
}

/* FMM */

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

  GetDP_Begin("LinAlg_MatVectorProdVector");

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

  GetDP_End;
}

/* Solve */

static void _solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, int precond){
  int its, RankCpu, i, j, view = 0;

  MPI_Comm_rank(PETSC_COMM_WORLD, &RankCpu);

  if(!Solver->ksp && !RankCpu) view = 1;

  if(view){
    ierr = MatGetSize(A->M, &i, &j); MYCHECK(ierr);
    Msg(PETSC, "N: %d", i);
  }

  if(!Solver->ksp) {
    ierr = KSPCreate(PETSC_COMM_WORLD, &Solver->ksp); MYCHECK(ierr);
    ierr = KSPSetOperators(Solver->ksp, A->M, A->M, DIFFERENT_NONZERO_PATTERN); MYCHECK(ierr);
    ierr = KSPGetPC(Solver->ksp, &Solver->pc); MYCHECK(ierr);
    if(Flag_FMM){
      ierr = PCSetType(Solver->pc, PCSHELL); MYCHECK(ierr);
      ierr = PCShellSetName(Solver->pc, "FMM"); MYCHECK(ierr);
      ierr = PCShellSetApply(Solver->pc, LinAlg_ApplyFMM); MYCHECK(ierr);
      ierr = KSPSetMonitor(Solver->ksp, LinAlg_ApplyFMMMonitor, PETSC_NULL, 0); MYCHECK(ierr);    
      ierr = KSPSetFromOptions(Solver->ksp); MYCHECK(ierr);
    }
    else{
      /* set some default options */
      ierr = PCSetType(Solver->pc, PCILU); MYCHECK(ierr);
      ierr = PCILUSetMatOrdering(Solver->pc, MATORDERING_RCM); MYCHECK(ierr);
      ierr = PCILUSetLevels(Solver->pc, 6); MYCHECK(ierr);
      ierr = KSPSetTolerances(Solver->ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, 
			      PETSC_DEFAULT); MYCHECK(ierr);
      
      /* override the default options with the ones from the option
	 database (if any) */
      ierr = KSPSetFromOptions(Solver->ksp); MYCHECK(ierr);
    }
  }
  else if(precond){
    ierr = KSPSetOperators(Solver->ksp, A->M, A->M, DIFFERENT_NONZERO_PATTERN); MYCHECK(ierr);
  }
  
  ierr = KSPSolve(Solver->ksp, B->V, X->V); MYCHECK(ierr);

  if(view){
    ierr = KSPView(Solver->ksp,PETSC_VIEWER_STDOUT_WORLD); MYCHECK(ierr);
  }

  if(!RankCpu){
    ierr = KSPGetIterationNumber(Solver->ksp, &its); MYCHECK(ierr);
    Msg(PETSC, "%d iterations", its);
  }

  GetDP_End;

}

void LinAlg_Solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X){
  GetDP_Begin("LinAlg_Solve");

  _solve(A, B, Solver, X, 1);

  GetDP_End;
}

void LinAlg_SolveAgain(gMatrix *A, gVector *B, gSolver *Solver, gVector *X){
  GetDP_Begin("LinAlg_SolveAgain");

  _solve(A, B, Solver, X, 0);

  GetDP_End;
}

#endif