13-10-2014 02

[GMRES2stage.git] / code / ex49.c

diff --git a/code/ex49.c b/code/ex49.c
index b875dd3aecef485cb850f1bbd7f08fd23beda96f..6e5331d9966195065efcd53aa6fb9ba3458c39a9 100644 (file)
--- a/code/ex49.c
+++ b/code/ex49.c
@@ -1,4 +1,7 @@
-//  /home/couturie/work/petsc-3.5.1/arch-linux2-c-debug/bin/mpirun -np 4    ./ex49 -mx 900 -my 900 -ksp_type fgmres
+//  /home/couturie/work/petsc-3.5.1_old/arch-linux2-c-debug/bin/mpirun -np 4  -machinefile archi  ./ex49 -mx 900 -my 900 -ksp_type fgmres -pc_type mg
+
+
+
 
 
 static char help[] =  "   Solves the compressible plane strain elasticity equations in 2d on the unit domain using Q1 finite elements. \n\
@@ -88,8 +91,8 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  gamma, alpha, oldgamma, beta;
-  PetscReal norm=20, Eprecision=1e-8, cgprec=1e-40;     
-  PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+  PetscReal norm=20, Eprecision=1e-3, cgprec=1e-40;     
+  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
@@ -97,7 +100,7 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   PetscInt size;
   PetscInt Istart,Iend;
   PetscInt i,its;
-  Vec       x_old, residu;
+  Vec       residu;
   Mat S, AS;
   PetscScalar *array;
   PetscInt *ind_row;
@@ -141,7 +144,6 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 
-  ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 
 
@@ -152,13 +154,14 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
 
   //GMRES WITH MINIMIZATION
   T1 = MPI_Wtime();
+  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
   while(giter<Emaxiter && norm>Eprecision ){
     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 
@@ -188,7 +191,6 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
 
 
       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
-      ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 
 
     }
@@ -246,6 +248,13 @@ int KrylovMinimize(Mat A, Vec b, Vec x) {
   T2 = MPI_Wtime();
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = VecDestroy(&r);CHKERRQ(ierr);  
+  ierr = VecDestroy(&vect);CHKERRQ(ierr);  
+  ierr = VecDestroy(&p);CHKERRQ(ierr);  
+  ierr = VecDestroy(&ss);CHKERRQ(ierr);  
+  ierr = VecDestroy(&Ax);CHKERRQ(ierr);  
+  ierr = VecDestroy(&residu);CHKERRQ(ierr);  
 
   return 0;
 
@@ -265,8 +274,8 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   //Variables
 
   PetscScalar  alpha, beta;
-  PetscReal norm=20, Eprecision=1e-8, tol=1e-40;     
-  PetscInt giter=0, ColS=8, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
+  PetscReal norm=20, Eprecision=1e-3, tol=1e-40;     
+  PetscInt giter=0, ColS=12, col=0, Emaxiter=50000000, iter=0, iterations=15, Iiter=0;
   PetscErrorCode ierr;
   PetscScalar T1, T2;
   KSP ksp;
@@ -274,7 +283,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   PetscInt size;
   PetscInt Istart,Iend;
   PetscInt i,its;
-  Vec       x_old, residu;
+  Vec       residu;
   Mat S, AS;
   PetscScalar *array;
   PetscInt *ind_row;
@@ -313,7 +322,6 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
 
-  ierr = VecDuplicate(b,&x_old);CHKERRQ(ierr);
   ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
 
 
@@ -342,7 +350,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //Initializations
   //  ierr = KSPGMRESSetRestart(ksp, 16); CHKERRQ(ierr);
-  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 16); CHKERRQ(ierr);
+  ierr = KSPSetTolerances(ksp, 1e-13, 1e-13, PETSC_DEFAULT, 30); CHKERRQ(ierr);
   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
 
 
@@ -350,6 +358,7 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
   //GMRES WITH MINIMIZATION
   T1 = MPI_Wtime();
+  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
   while(giter<Emaxiter && norm>Eprecision ){
     for(col=0; col<ColS  &&  norm>Eprecision; col++){
 
@@ -379,12 +388,12 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
 
 
       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D\n", norm, giter); CHKERRQ(ierr);
-      ierr = VecCopy(x, x_old); CHKERRQ(ierr);
 
 
     }
 
 
+
     //minimization step
     if( norm>Eprecision) {
 
@@ -477,6 +486,16 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) {
   T2 = MPI_Wtime();
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time LSQR : %g (s)\n", T2-T1); CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations LSQR : %D\n", total); CHKERRQ(ierr);
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = VecDestroy(&Ax);CHKERRQ(ierr);  
+  ierr = VecDestroy(&u);CHKERRQ(ierr);  
+  ierr = VecDestroy(&uu);CHKERRQ(ierr);  
+  ierr = VecDestroy(&zero_long);CHKERRQ(ierr);  
+  ierr = VecDestroy(&d);CHKERRQ(ierr);  
+  ierr = VecDestroy(&residu);CHKERRQ(ierr);  
+  ierr = VecDestroy(&vt);CHKERRQ(ierr);  
+  ierr = VecDestroy(&x_lsqr);CHKERRQ(ierr);  
+
 
   return 0;
 
@@ -1420,11 +1439,11 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
   ierr = DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);CHKERRQ(ierr);
   ierr = DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);CHKERRQ(ierr);
 
-  ierr = PetscOptionsGetBool(NULL,"-no_view",&no_view,NULL);CHKERRQ(ierr);
+  /*  ierr = PetscOptionsGetBool(NULL,"-no_view",&no_view,NULL);CHKERRQ(ierr);
   if (!no_view) {
     ierr = DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for elasticity eqn.","properties");CHKERRQ(ierr);
     ierr = DMDACoordViewGnuplot2d(elas_da,"mesh");CHKERRQ(ierr);
-  }
+   }*/
 
   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
   ierr = DMSetMatType(elas_da,MATAIJ);CHKERRQ(ierr);
@@ -1445,9 +1464,9 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
 
 
   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_E);CHKERRQ(ierr);
-  ierr = KSPSetOptionsPrefix(ksp_E,"elas_");CHKERRQ(ierr);  /* elasticity */
+  // ierr = KSPSetOptionsPrefix(ksp_E,"elas_");CHKERRQ(ierr);  /* elasticity */
 
-  ierr = PetscOptionsGetBool(NULL,"-use_nonsymbc",&use_nonsymbc,&flg);CHKERRQ(ierr);
+  //ierr = PetscOptionsGetBool(NULL,"-use_nonsymbc",&use_nonsymbc,&flg);CHKERRQ(ierr);
   /* solve */
   if (!use_nonsymbc) {
     Mat        AA;
@@ -1461,15 +1480,30 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
     ierr = KSPSetOperators(ksp_E,AA,AA);CHKERRQ(ierr);
     ierr = KSPSetFromOptions(ksp_E);CHKERRQ(ierr);
 
- ierr = KSPSetFromOptions(ksp_E);CHKERRQ(ierr);
 
    PetscScalar T1,T2;
-    ierr = KSPSetTolerances(ksp_E, 1e-9, 1e-9, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
-    T1 = MPI_Wtime();
+    ierr = KSPSetTolerances(ksp_E, 1e-7, 1e-7, PETSC_DEFAULT, 50000000); CHKERRQ(ierr);
+
+
+    PC             pc;
+    KSPGetPC(ksp_E, &pc);
+    PCType         type;
+    PCGetType(pc, &type);
 
+    PetscPrintf(PETSC_COMM_WORLD, "PC TYPE %s  \n", type);
+    KSPGetType(ksp_E,&type);
+    PetscPrintf(PETSC_COMM_WORLD, "SOLVER TYPE %s  \n", type);
+
+    /*
+    T1 = MPI_Wtime();
+ ierr = KSPSetUp(ksp_E);CHKERRQ(ierr);
     ierr = KSPSolve(ksp_E,ff,XX);CHKERRQ(ierr);
     T2 = MPI_Wtime();
-    
+     */    
+
+
+
+
     Mat A;
     Vec sol;
     PetscScalar norm;
@@ -1481,20 +1515,65 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
 
 
     KSPGetOperators(ksp_E,&A,NULL);
-    MatMult(A,XX,sol);
+    /* MatMult(A,XX,sol);
     VecAXPY(sol,-1,ff);
     VecNorm(sol, NORM_2, &norm);
-
+    
 
     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Norm of error : %g\n", (double)norm); CHKERRQ(ierr);
     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
+
+     */
+
+
+        //
+    //
+    //version to control the error
+ {
+
+    Vec x2;
+    Vec sol;
+    VecDuplicate(ff,&x2);
+    VecDuplicate(ff,&sol);
+    
+    PetscScalar norm=100;
+    PetscScalar T1,T2;
+    PetscInt total,its;
+    ierr = KSPSetTolerances(ksp_E,1e-10,1e-10,PETSC_DEFAULT,30);CHKERRQ(ierr);
+    ierr = KSPSetInitialGuessNonzero(ksp_E, PETSC_TRUE); CHKERRQ(ierr);
+    T1 = MPI_Wtime();
+    while(norm>1e-3) {
+      ierr = KSPSolve(ksp_E,ff,x2);CHKERRQ(ierr);
+      KSPGetResidualNorm(ksp_E,&norm);
+      ierr = KSPGetIterationNumber(ksp_E, &its); CHKERRQ(ierr);
+      total += its;
+       ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g\n", norm); CHKERRQ(ierr);
+    }
+
+    T2 = MPI_Wtime();
+
+    MatMult(A,x2,sol);
+    VecAXPY(sol,-1,ff);
+    VecNorm(sol, NORM_2, &norm);
+     ierr = PetscPrintf(PETSC_COMM_WORLD,"Computed norm of error %g iterations %D\n",(double)norm,total);CHKERRQ(ierr);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time NORMAL GMRES : %g (s)\n\n\n", T2-T1); CHKERRQ(ierr);
+
+    ierr = KSPDestroy(&ksp_E);CHKERRQ(ierr);
+    ierr = VecDestroy(&x2);CHKERRQ(ierr);
+    ierr = VecDestroy(&sol);CHKERRQ(ierr);
+  }
+     
+
+
+
+ /*
     PetscInt total;
     ierr = KSPGetIterationNumber(ksp_E, &total); CHKERRQ(ierr);
     ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
+  */
 
 
-
-
+  VecCopy(ff,X2);
     KrylovMinimize(A, ff, X2);
     MatMult(A,X2,sol);
 
@@ -1538,7 +1617,7 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my)
     ierr = KSPSolve(ksp_E,f,X);CHKERRQ(ierr);
   }
 
-  if (!no_view) {ierr = DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");CHKERRQ(ierr);}
+  //  if (!no_view) {ierr = DMDAViewGnuplot2d(elas_da,X,"Displacement solution for elasticity eqn.","X");CHKERRQ(ierr);}
   ierr = KSPDestroy(&ksp_E);CHKERRQ(ierr);
 
   ierr = VecDestroy(&X);CHKERRQ(ierr);
@@ -1566,6 +1645,12 @@ int main(int argc,char **args)
   ierr = PetscOptionsGetInt(NULL,"-mx",&mx,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsGetInt(NULL,"-my",&my,NULL);CHKERRQ(ierr);
 
+
+  PetscMPIInt size;
+  MPI_Comm_size(PETSC_COMM_WORLD,&size);
+  PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d\n",size);
+
+
   ierr = solve_elasticity_2d(mx,my);CHKERRQ(ierr);
 
   ierr = PetscFinalize();CHKERRQ(ierr);