From 53346b46c15cba98de051c07c0c720a8cf2ab8ca Mon Sep 17 00:00:00 2001 From: raphael couturier Date: Wed, 1 Oct 2014 11:05:37 +0200 Subject: [PATCH 1/1] update ex49 --- code/ex49.c | 100 +++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 83 insertions(+), 17 deletions(-) diff --git a/code/ex49.c b/code/ex49.c index b93b5c6..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-6, 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); @@ -189,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); } @@ -247,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; @@ -266,8 +274,8 @@ int KrylovMinimizeLSQR(Mat A, Vec b, Vec x) { //Variables PetscScalar alpha, beta; - PetscReal norm=20, Eprecision=1e-6, 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; @@ -275,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; @@ -314,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); @@ -381,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) { @@ -479,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; @@ -1477,12 +1494,16 @@ static PetscErrorCode solve_elasticity_2d(PetscInt mx,PetscInt my) 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; @@ -1494,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); -- 2.39.5