X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/GMRES2stage.git/blobdiff_plain/1ecb79e1234256ff879d1a787113c41fca16a7cc..9fe58147d9a4e7a4182385b4a397fe96695b7862:/code/ex49.c diff --git a/code/ex49.c b/code/ex49.c index b875dd3..6e5331d 100644 --- 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(giterEprecision ){ for(col=0; colEprecision; 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(giterEprecision ){ for(col=0; colEprecision; 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);