Actual source code: ex39.c

  1: /*
  2: mpiexec -n 8 ./ex39 -ksp_type fbcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -ksp_monitor -n1 32 -n2 32 -n3 32

  4:   Contributed by Jie Chen for testing flexible BiCGStab algorithm
  5: */

  7: static char help[] = "Solves the PDE (in 3D) - laplacian(u) + gamma x dot grad(u) + beta u = 1\n\
  8: with zero Dirichlet condition. The discretization is standard centered\n\
  9: difference. Input parameters include:\n\
 10:   -n1        : number of mesh points in 1st dimension (default 32)\n\
 11:   -n2        : number of mesh points in 2nd dimension (default 32)\n\
 12:   -n3        : number of mesh points in 3rd dimension (default 32)\n\
 13:   -h         : spacing between mesh points (default 1/n1)\n\
 14:   -gamma     : gamma (default 4/h)\n\
 15:   -beta      : beta (default 0.01/h^2)\n\n";

 17: #include <petscksp.h>
 18: int main(int argc,char **args)
 19: {
 20:   Vec            x,b,u;                 /* approx solution, RHS, working vector */
 21:   Mat            A;                     /* linear system matrix */
 22:   KSP            ksp;                   /* linear solver context */
 23:   PetscInt       n1, n2, n3;            /* parameters */
 24:   PetscReal      h, gamma, beta;        /* parameters */
 25:   PetscInt       i,j,k,Ii,J,Istart,Iend;
 26:   PetscScalar    v, co1, co2;

 28:   PetscInitialize(&argc,&args,(char*)0,help);
 29:   n1 = 32;
 30:   n2 = 32;
 31:   n3 = 32;
 32:   PetscOptionsGetInt(NULL,NULL,"-n1",&n1,NULL);
 33:   PetscOptionsGetInt(NULL,NULL,"-n2",&n2,NULL);
 34:   PetscOptionsGetInt(NULL,NULL,"-n3",&n3,NULL);

 36:   h     = 1.0/n1;
 37:   gamma = 4.0/h;
 38:   beta  = 0.01/(h*h);
 39:   PetscOptionsGetReal(NULL,NULL,"-h",&h,NULL);
 40:   PetscOptionsGetReal(NULL,NULL,"-gamma",&gamma,NULL);
 41:   PetscOptionsGetReal(NULL,NULL,"-beta",&beta,NULL);

 43:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 44:          Compute the matrix and set right-hand-side vector.
 45:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 46:   MatCreate(PETSC_COMM_WORLD,&A);
 47:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n1*n2*n3,n1*n2*n3);
 48:   MatSetFromOptions(A);
 49:   MatMPIAIJSetPreallocation(A,7,NULL,7,NULL);
 50:   MatSeqAIJSetPreallocation(A,7,NULL);
 51:   MatSetUp(A);
 52:   MatGetOwnershipRange(A,&Istart,&Iend);

 54:   /*
 55:      Set matrix elements for the 3-D, seven-point stencil in parallel.
 56:       - Each processor needs to insert only elements that it owns
 57:         locally (but any non-local elements will be sent to the
 58:         appropriate processor during matrix assembly).
 59:       - Always specify global rows and columns of matrix entries.
 60:    */
 61:   co1  = gamma * h * h / 2.0;
 62:   co2  = beta * h * h;
 63:   for (Ii=Istart; Ii<Iend; Ii++) {
 64:     i = Ii/(n2*n3); j = (Ii - i*n2*n3)/n3; k = Ii - i*n2*n3 - j*n3;
 65:     if (i>0) {
 66:       J    = Ii - n2*n3;  v = -1.0 + co1*(PetscScalar)i;
 67:       MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 68:     }
 69:     if (i<n1-1) {
 70:       J    = Ii + n2*n3;  v = -1.0 + co1*(PetscScalar)i;
 71:       MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 72:     }
 73:     if (j>0) {
 74:       J    = Ii - n3;  v = -1.0 + co1*(PetscScalar)j;
 75:       MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 76:     }
 77:     if (j<n2-1) {
 78:       J    = Ii + n3;  v = -1.0 + co1*(PetscScalar)j;
 79:       MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 80:     }
 81:     if (k>0) {
 82:       J    = Ii - 1;  v = -1.0 + co1*(PetscScalar)k;
 83:       MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 84:     }
 85:     if (k<n3-1) {
 86:       J    = Ii + 1;  v = -1.0 + co1*(PetscScalar)k;
 87:       MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
 88:     }
 89:     v    = 6.0 + co2;
 90:     MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 91:   }
 92:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 93:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 95:   /* Create parallel vectors and Set right-hand side. */
 96:   VecCreate(PETSC_COMM_WORLD,&b);
 97:   VecSetSizes(b,PETSC_DECIDE,n1*n2*n3);
 98:   VecSetFromOptions(b);
 99:   VecDuplicate(b,&x);
100:   VecDuplicate(b,&u);
101:   VecSet(b,1.0);

103:   /* Create linear solver context */
104:   KSPCreate(PETSC_COMM_WORLD,&ksp);
105:   KSPSetOperators(ksp,A,A);
106:   KSPSetTolerances(ksp,1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,200);
107:   KSPSetFromOptions(ksp);

109:   /* Solve the linear system */
110:   KSPSolve(ksp,b,x);

112:   /* Free work space.  */
113:   KSPDestroy(&ksp);
114:   VecDestroy(&u));  PetscCall(VecDestroy(&x);
115:   VecDestroy(&b));  PetscCall(MatDestroy(&A);
116:   PetscFinalize();
117:   return 0;
118: }

120: /*TEST

122:    test:
123:       nsize: 8
124:       args: -ksp_type fbcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 32 -n2 32 -n3 32

126:    test:
127:       suffix: 2
128:       nsize: 8
129:       args: -ksp_type fbcgsr -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 32 -n2 32 -n3 32
130:       output_file: output/ex39_1.out

132:    test:
133:       suffix: 3
134:       nsize: 8
135:       args: -ksp_type qmrcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 32 -n2 32 -n3 32
136:       output_file: output/ex39_1.out

138: TEST*/