1 2 static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\ 3 matrix-free techniques with user-provided explicit preconditioner matrix.\n\n"; 4 5 #include <petscsnes.h> 6 7 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 8 extern PetscErrorCode FormJacobianNoMatrix(SNES,Vec,Mat,Mat,void*); 9 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 10 extern PetscErrorCode FormFunctioni(void *,PetscInt,Vec,PetscScalar *); 11 extern PetscErrorCode OtherFunctionForDifferencing(void*,Vec,Vec); 12 extern PetscErrorCode FormInitialGuess(SNES,Vec); 13 extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*); 14 15 typedef struct { 16 PetscViewer viewer; 17 } MonitorCtx; 18 19 typedef struct { 20 PetscBool variant; 21 } AppCtx; 22 23 int main(int argc,char **argv) 24 { 25 SNES snes; /* SNES context */ 26 SNESType type = SNESNEWTONLS; /* default nonlinear solution method */ 27 Vec x,r,F,U; /* vectors */ 28 Mat J,B; /* Jacobian matrix-free, explicit preconditioner */ 29 AppCtx user; /* user-defined work context */ 30 PetscScalar h,xp = 0.0,v; 31 PetscInt its,n = 5,i; 32 PetscErrorCode ierr; 33 PetscBool puremf = PETSC_FALSE; 34 35 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 36 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 37 ierr = PetscOptionsHasName(NULL,NULL,"-variant",&user.variant);CHKERRQ(ierr); 38 h = 1.0/(n-1); 39 40 /* Set up data structures */ 41 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&x);CHKERRQ(ierr); 42 ierr = PetscObjectSetName((PetscObject)x,"Approximate Solution");CHKERRQ(ierr); 43 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 44 ierr = VecDuplicate(x,&F);CHKERRQ(ierr); 45 ierr = VecDuplicate(x,&U);CHKERRQ(ierr); 46 ierr = PetscObjectSetName((PetscObject)U,"Exact Solution");CHKERRQ(ierr); 47 48 /* create explicit matrix preconditioner */ 49 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&B);CHKERRQ(ierr); 50 51 /* Store right-hand-side of PDE and exact solution */ 52 for (i=0; i<n; i++) { 53 v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */ 54 ierr = VecSetValues(F,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 55 v = xp*xp*xp; 56 ierr = VecSetValues(U,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 57 xp += h; 58 } 59 60 /* Create nonlinear solver */ 61 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 62 ierr = SNESSetType(snes,type);CHKERRQ(ierr); 63 64 /* Set various routines and options */ 65 ierr = SNESSetFunction(snes,r,FormFunction,F);CHKERRQ(ierr); 66 if (user.variant) { 67 /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */ 68 ierr = MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);CHKERRQ(ierr); 69 ierr = MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);CHKERRQ(ierr); 70 ierr = MatMFFDSetFunctioni(J,FormFunctioni);CHKERRQ(ierr); 71 /* Use the matrix free operator for both the Jacobian used to define the linear system and used to define the preconditioner */ 72 /* This tests MatGetDiagonal() for MATMFFD */ 73 ierr = PetscOptionsHasName(NULL,NULL,"-puremf",&puremf);CHKERRQ(ierr); 74 } else { 75 /* create matrix free matrix for Jacobian */ 76 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 77 /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */ 78 /* note we use the same context for this function as FormFunction, the F vector */ 79 ierr = MatMFFDSetFunction(J,OtherFunctionForDifferencing,F);CHKERRQ(ierr); 80 } 81 82 /* Set various routines and options */ 83 ierr = SNESSetJacobian(snes,J,puremf ? J : B,puremf ? FormJacobianNoMatrix : FormJacobian,&user);CHKERRQ(ierr); 84 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 85 86 /* Solve nonlinear system */ 87 ierr = FormInitialGuess(snes,x);CHKERRQ(ierr); 88 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); 89 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 90 ierr = PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);CHKERRQ(ierr); 91 92 /* Free data structures */ 93 ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); 94 ierr = VecDestroy(&U);CHKERRQ(ierr); ierr = VecDestroy(&F);CHKERRQ(ierr); 95 ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = MatDestroy(&B);CHKERRQ(ierr); 96 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 97 ierr = PetscFinalize(); 98 return ierr; 99 } 100 /* -------------------- Evaluate Function F(x) --------------------- */ 101 102 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy) 103 { 104 const PetscScalar *xx,*FF; 105 PetscScalar *ff,d; 106 PetscInt i,n; 107 PetscErrorCode ierr; 108 109 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 110 ierr = VecGetArray(f,&ff);CHKERRQ(ierr); 111 ierr = VecGetArrayRead((Vec) dummy,&FF);CHKERRQ(ierr); 112 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 113 d = (PetscReal)(n - 1); d = d*d; 114 ff[0] = xx[0]; 115 for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i]; 116 ff[n-1] = xx[n-1] - 1.0; 117 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 118 ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr); 119 ierr = VecRestoreArrayRead((Vec)dummy,&FF);CHKERRQ(ierr); 120 return 0; 121 } 122 123 PetscErrorCode FormFunctioni(void *dummy,PetscInt i,Vec x,PetscScalar *s) 124 { 125 const PetscScalar *xx,*FF; 126 PetscScalar d; 127 PetscInt n; 128 PetscErrorCode ierr; 129 SNES snes = (SNES) dummy; 130 Vec F; 131 132 ierr = SNESGetFunction(snes,NULL,NULL,(void**)&F);CHKERRQ(ierr); 133 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 134 ierr = VecGetArrayRead(F,&FF);CHKERRQ(ierr); 135 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 136 d = (PetscReal)(n - 1); d = d*d; 137 if (i == 0) { 138 *s = xx[0]; 139 } else if (i == n-1) { 140 *s = xx[n-1] - 1.0; 141 } else { 142 *s = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i]; 143 } 144 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 145 ierr = VecRestoreArrayRead(F,&FF);CHKERRQ(ierr); 146 return 0; 147 } 148 149 /* 150 151 Example function that when differenced produces the same matrix free Jacobian as FormFunction() 152 this is provided to show how a user can provide a different function 153 */ 154 PetscErrorCode OtherFunctionForDifferencing(void *dummy,Vec x,Vec f) 155 { 156 PetscErrorCode ierr; 157 158 ierr = FormFunction(NULL,x,f,dummy);CHKERRQ(ierr); 159 ierr = VecShift(f,1.0);CHKERRQ(ierr); 160 return 0; 161 } 162 163 /* -------------------- Form initial approximation ----------------- */ 164 165 PetscErrorCode FormInitialGuess(SNES snes,Vec x) 166 { 167 PetscErrorCode ierr; 168 PetscScalar pfive = .50; 169 ierr = VecSet(x,pfive);CHKERRQ(ierr); 170 return 0; 171 } 172 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 173 /* Evaluates a matrix that is used to precondition the matrix-free 174 jacobian. In this case, the explicit preconditioner matrix is 175 also EXACTLY the Jacobian. In general, it would be some lower 176 order, simplified apprioximation */ 177 178 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 179 { 180 const PetscScalar *xx; 181 PetscScalar A[3],d; 182 PetscInt i,n,j[3]; 183 PetscErrorCode ierr; 184 AppCtx *user = (AppCtx*) dummy; 185 186 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 187 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 188 d = (PetscReal)(n - 1); d = d*d; 189 190 i = 0; A[0] = 1.0; 191 ierr = MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);CHKERRQ(ierr); 192 for (i=1; i<n-1; i++) { 193 j[0] = i - 1; j[1] = i; j[2] = i + 1; 194 A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d; 195 ierr = MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr); 196 } 197 i = n-1; A[0] = 1.0; 198 ierr = MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);CHKERRQ(ierr); 199 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 200 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 201 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 202 203 if (user->variant) { 204 ierr = MatMFFDSetBase(jac,x,NULL);CHKERRQ(ierr); 205 } 206 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 207 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 208 return 0; 209 } 210 211 PetscErrorCode FormJacobianNoMatrix(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 212 { 213 PetscErrorCode ierr; 214 AppCtx *user = (AppCtx*) dummy; 215 216 if (user->variant) { 217 ierr = MatMFFDSetBase(jac,x,NULL);CHKERRQ(ierr); 218 } 219 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 220 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 221 return 0; 222 } 223 224 /* -------------------- User-defined monitor ----------------------- */ 225 226 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy) 227 { 228 PetscErrorCode ierr; 229 MonitorCtx *monP = (MonitorCtx*) dummy; 230 Vec x; 231 MPI_Comm comm; 232 233 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 234 ierr = PetscFPrintf(comm,stdout,"iter = %D, SNES Function norm %g \n",its,(double)fnorm);CHKERRQ(ierr); 235 ierr = SNESGetSolution(snes,&x);CHKERRQ(ierr); 236 ierr = VecView(x,monP->viewer);CHKERRQ(ierr); 237 return 0; 238 } 239 240 /*TEST 241 242 test: 243 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 244 245 test: 246 suffix: 2 247 args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 248 output_file: output/ex7_1.out 249 250 # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning 251 test: 252 suffix: 3 253 args: -variant -pc_type jacobi -snes_view -ksp_monitor 254 255 # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning 256 test: 257 suffix: 4 258 args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor 259 260 TEST*/ 261