Dear all,
Sorry to disturb you. I am a user of Petsc. I am trying to use Fieldsplit in Petsc to do preconditioning for Navier-Stokes problem. I have some problems when I trying to use Fieldsplit function. I am now defining the nest matrix first, then I get the IS from nested matrix. But I find that my code just work for one core. When I change to parallel case, I could only get zero solution. I wonder is there any special requirements for IS definition in Fieldsplit? I include one code here. If you have any idea, hope you reply soon. Thank you for your help. Thank you very much. Yours sincerely, Qiming Zhu,
makefile
Description: makefile
static char help[] = "Test";
#include <petscksp.h>
#include <petscpc.h>
int main(int argc,char **args)
{
Mat Asub[4];
Mat Anest;
Vec x,b,y;
IS isg[2];
PetscInt row,end,i,j,Ist,Ien;
PetscErrorCode ierr;
PetscMPIInt rank,size;
KSP ksp;
PC pc;
PetscInt kk;
// Initialize petsc and mpi here
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
/* My matrix and rhs vec are as followed:
The matrix Anest I define here is:
1 0 0 0 | 0 0 0 0
0 2 0 0 | 0 0 0 0
0 0 3 0 | 0 0 0 0
0 0 0 4 | 0 0 0 0
- - - - - - - - -
0 0 0 0 | 5 0 0 0
0 0 0 0 | 0 6 0 0
0 0 0 0 | 0 0 7 0
0 0 0 0 | 0 0 0 8
Thr rhs I define here is :
1
2
3
4
-
5
6
7
8
So my solution should be all one:
1
1
1
1
1
1
1
1
*/
// kk is submatrix size
kk=4;
MatCreate(PETSC_COMM_WORLD,&Asub[0]);
MatSetSizes(Asub[0],PETSC_DECIDE,PETSC_DECIDE,kk,kk);
MatSetType(Asub[0],MATMPIAIJ);
MatMPIAIJSetPreallocation(Asub[0],5,NULL,5,NULL);
MatGetOwnershipRange(Asub[0], &Ist, &Ien);
for(i=Ist;i<Ien;i++){
MatSetValue(Asub[0],i,i,i+1.0,INSERT_VALUES);
}
ierr = MatAssemblyBegin(Asub[0],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Asub[0],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
MatCreate(PETSC_COMM_WORLD,&Asub[1]);
MatSetSizes(Asub[1],PETSC_DECIDE,PETSC_DECIDE,kk,kk);
MatSetType(Asub[1],MATMPIAIJ);
MatMPIAIJSetPreallocation(Asub[1],5,NULL,5,NULL);
ierr = MatAssemblyBegin(Asub[1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Asub[1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
MatCreate(PETSC_COMM_WORLD,&Asub[2]);
MatSetSizes(Asub[2],PETSC_DECIDE,PETSC_DECIDE,kk,kk);
MatSetType(Asub[2],MATMPIAIJ);
MatMPIAIJSetPreallocation(Asub[2],5,NULL,5,NULL);
ierr = MatAssemblyBegin(Asub[2],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Asub[2],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
MatCreate(PETSC_COMM_WORLD,&Asub[3]);
MatSetSizes(Asub[3],PETSC_DECIDE,PETSC_DECIDE,kk,kk);
MatSetType(Asub[3],MATMPIAIJ);
MatMPIAIJSetPreallocation(Asub[3],5,NULL,5,NULL);
MatGetOwnershipRange(Asub[3], &Ist, &Ien);
for(i=Ist;i<Ien;i++){
MatSetValue(Asub[3],i,i,i+kk*1.0+1.0,INSERT_VALUES);
}
ierr = MatAssemblyBegin(Asub[3],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Asub[3],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
// Define nest matrix here
MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, Asub, &Anest);
ierr = MatAssemblyBegin(Anest,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(Anest,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
// Get IS from nest matrix here
MatNestGetISs(Anest,isg,NULL);
// Create vec here
VecCreate(PETSC_COMM_WORLD, &b);
VecSetSizes(b, PETSC_DECIDE, 2*kk);
VecSetType(b, VECMPI);
// Assign vec values
VecGetOwnershipRange(b, &Ist, &Ien);
for(i=Ist;i<Ien;i++){
VecSetValue(b,i,i+1.0,INSERT_VALUES);
}
VecDuplicate(b, &x);
//VecSetRandom(x, NULL);
// Assemble vec here
ierr = VecAssemblyBegin(b);
ierr = VecAssemblyEnd(b);
ierr = VecAssemblyBegin(x);
ierr = VecAssemblyEnd(x);
//VecView(x,PETSC_VIEWER_STDOUT_WORLD);
// Print out matrix and rhs here:
printf("A00 block print here.\n");
MatView(Asub[0],PETSC_VIEWER_STDOUT_WORLD);
printf("A01 block print here.\n");
MatView(Asub[1],PETSC_VIEWER_STDOUT_WORLD);
printf("A10 block print here.\n");
MatView(Asub[2],PETSC_VIEWER_STDOUT_WORLD);
printf("A11 block print here.\n");
MatView(Asub[3],PETSC_VIEWER_STDOUT_WORLD);
VecView(b,PETSC_VIEWER_STDOUT_WORLD);
// Print out IS here:
ISView(isg[0],PETSC_VIEWER_STDOUT_WORLD);
ISView(isg[1],PETSC_VIEWER_STDOUT_WORLD);
// Define solver here
KSPCreate(PETSC_COMM_WORLD, &ksp);
//KSPSetType(ksp,KSPFGMRES);
KSPSetOperators(ksp, Anest, Anest);
KSPSetFromOptions(ksp);
// I am tring to use basic fieldsplit preconditioner here:
KSPGetPC(ksp, &pc);
PCSetType(pc,PCFIELDSPLIT);
PCFieldSplitSetIS(pc, "0", isg[0]);
PCFieldSplitSetIS(pc, "1", isg[1]);
//PCFieldSplitSetType(pc,PC_COMPOSITE_ADDITIVE);
//ierr = PetscOptionsSetValue(NULL,"-fieldsplit_0_ksp_type","preonly");
//ierr = PetscOptionsSetValue(NULL,"-fieldsplit_0_pc_type","jacobi");
//ierr = PetscOptionsSetValue(NULL,"-fieldsplit_1_ksp_type","preonly");
//ierr = PetscOptionsSetValue(NULL,"-fieldsplit_1_pc_type","jacobi");
PCSetUp(pc);
// Solve linear equation and print solution
KSPSolve(ksp,b,x);
KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
VecView(x,PETSC_VIEWER_STDOUT_WORLD);
ierr = PetscFinalize();
return ierr;
}
