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 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 37 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-variant",&user.variant)); 38 h = 1.0/(n-1); 39 40 /* Set up data structures */ 41 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&x)); 42 CHKERRQ(PetscObjectSetName((PetscObject)x,"Approximate Solution")); 43 CHKERRQ(VecDuplicate(x,&r)); 44 CHKERRQ(VecDuplicate(x,&F)); 45 CHKERRQ(VecDuplicate(x,&U)); 46 CHKERRQ(PetscObjectSetName((PetscObject)U,"Exact Solution")); 47 48 /* create explicit matrix preconditioner */ 49 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&B)); 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 CHKERRQ(VecSetValues(F,1,&i,&v,INSERT_VALUES)); 55 v = xp*xp*xp; 56 CHKERRQ(VecSetValues(U,1,&i,&v,INSERT_VALUES)); 57 xp += h; 58 } 59 60 /* Create nonlinear solver */ 61 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 62 CHKERRQ(SNESSetType(snes,type)); 63 64 /* Set various routines and options */ 65 CHKERRQ(SNESSetFunction(snes,r,FormFunction,F)); 66 if (user.variant) { 67 /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */ 68 CHKERRQ(MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J)); 69 CHKERRQ(MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes)); 70 CHKERRQ(MatMFFDSetFunctioni(J,FormFunctioni)); 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 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-puremf",&puremf)); 74 } else { 75 /* create matrix free matrix for Jacobian */ 76 CHKERRQ(MatCreateSNESMF(snes,&J)); 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 CHKERRQ(MatMFFDSetFunction(J,OtherFunctionForDifferencing,F)); 80 } 81 82 /* Set various routines and options */ 83 CHKERRQ(SNESSetJacobian(snes,J,puremf ? J : B,puremf ? FormJacobianNoMatrix : FormJacobian,&user)); 84 CHKERRQ(SNESSetFromOptions(snes)); 85 86 /* Solve nonlinear system */ 87 CHKERRQ(FormInitialGuess(snes,x)); 88 CHKERRQ(SNESSolve(snes,NULL,x)); 89 CHKERRQ(SNESGetIterationNumber(snes,&its)); 90 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its)); 91 92 /* Free data structures */ 93 CHKERRQ(VecDestroy(&x)); CHKERRQ(VecDestroy(&r)); 94 CHKERRQ(VecDestroy(&U)); CHKERRQ(VecDestroy(&F)); 95 CHKERRQ(MatDestroy(&J)); CHKERRQ(MatDestroy(&B)); 96 CHKERRQ(SNESDestroy(&snes)); 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 108 CHKERRQ(VecGetArrayRead(x,&xx)); 109 CHKERRQ(VecGetArray(f,&ff)); 110 CHKERRQ(VecGetArrayRead((Vec) dummy,&FF)); 111 CHKERRQ(VecGetSize(x,&n)); 112 d = (PetscReal)(n - 1); d = d*d; 113 ff[0] = xx[0]; 114 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]; 115 ff[n-1] = xx[n-1] - 1.0; 116 CHKERRQ(VecRestoreArrayRead(x,&xx)); 117 CHKERRQ(VecRestoreArray(f,&ff)); 118 CHKERRQ(VecRestoreArrayRead((Vec)dummy,&FF)); 119 return 0; 120 } 121 122 PetscErrorCode FormFunctioni(void *dummy,PetscInt i,Vec x,PetscScalar *s) 123 { 124 const PetscScalar *xx,*FF; 125 PetscScalar d; 126 PetscInt n; 127 SNES snes = (SNES) dummy; 128 Vec F; 129 130 CHKERRQ(SNESGetFunction(snes,NULL,NULL,(void**)&F)); 131 CHKERRQ(VecGetArrayRead(x,&xx)); 132 CHKERRQ(VecGetArrayRead(F,&FF)); 133 CHKERRQ(VecGetSize(x,&n)); 134 d = (PetscReal)(n - 1); d = d*d; 135 if (i == 0) { 136 *s = xx[0]; 137 } else if (i == n-1) { 138 *s = xx[n-1] - 1.0; 139 } else { 140 *s = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i]; 141 } 142 CHKERRQ(VecRestoreArrayRead(x,&xx)); 143 CHKERRQ(VecRestoreArrayRead(F,&FF)); 144 return 0; 145 } 146 147 /* 148 149 Example function that when differenced produces the same matrix free Jacobian as FormFunction() 150 this is provided to show how a user can provide a different function 151 */ 152 PetscErrorCode OtherFunctionForDifferencing(void *dummy,Vec x,Vec f) 153 { 154 155 CHKERRQ(FormFunction(NULL,x,f,dummy)); 156 CHKERRQ(VecShift(f,1.0)); 157 return 0; 158 } 159 160 /* -------------------- Form initial approximation ----------------- */ 161 162 PetscErrorCode FormInitialGuess(SNES snes,Vec x) 163 { 164 PetscScalar pfive = .50; 165 CHKERRQ(VecSet(x,pfive)); 166 return 0; 167 } 168 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 169 /* Evaluates a matrix that is used to precondition the matrix-free 170 jacobian. In this case, the explicit preconditioner matrix is 171 also EXACTLY the Jacobian. In general, it would be some lower 172 order, simplified apprioximation */ 173 174 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 175 { 176 const PetscScalar *xx; 177 PetscScalar A[3],d; 178 PetscInt i,n,j[3]; 179 AppCtx *user = (AppCtx*) dummy; 180 181 CHKERRQ(VecGetArrayRead(x,&xx)); 182 CHKERRQ(VecGetSize(x,&n)); 183 d = (PetscReal)(n - 1); d = d*d; 184 185 i = 0; A[0] = 1.0; 186 CHKERRQ(MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES)); 187 for (i=1; i<n-1; i++) { 188 j[0] = i - 1; j[1] = i; j[2] = i + 1; 189 A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d; 190 CHKERRQ(MatSetValues(B,1,&i,3,j,A,INSERT_VALUES)); 191 } 192 i = n-1; A[0] = 1.0; 193 CHKERRQ(MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES)); 194 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 195 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 196 CHKERRQ(VecRestoreArrayRead(x,&xx)); 197 198 if (user->variant) { 199 CHKERRQ(MatMFFDSetBase(jac,x,NULL)); 200 } 201 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 202 CHKERRQ(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 203 return 0; 204 } 205 206 PetscErrorCode FormJacobianNoMatrix(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 207 { 208 AppCtx *user = (AppCtx*) dummy; 209 210 if (user->variant) { 211 CHKERRQ(MatMFFDSetBase(jac,x,NULL)); 212 } 213 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 214 CHKERRQ(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 215 return 0; 216 } 217 218 /* -------------------- User-defined monitor ----------------------- */ 219 220 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy) 221 { 222 MonitorCtx *monP = (MonitorCtx*) dummy; 223 Vec x; 224 MPI_Comm comm; 225 226 CHKERRQ(PetscObjectGetComm((PetscObject)snes,&comm)); 227 CHKERRQ(PetscFPrintf(comm,stdout,"iter = %D, SNES Function norm %g \n",its,(double)fnorm)); 228 CHKERRQ(SNESGetSolution(snes,&x)); 229 CHKERRQ(VecView(x,monP->viewer)); 230 return 0; 231 } 232 233 /*TEST 234 235 test: 236 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 237 238 test: 239 suffix: 2 240 args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 241 output_file: output/ex7_1.out 242 243 # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning 244 test: 245 suffix: 3 246 args: -variant -pc_type jacobi -snes_view -ksp_monitor 247 248 # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning 249 test: 250 suffix: 4 251 args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor 252 253 TEST*/ 254