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