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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 35 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 36 PetscCall(PetscOptionsHasName(NULL,NULL,"-variant",&user.variant)); 37 h = 1.0/(n-1); 38 39 /* Set up data structures */ 40 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,&x)); 41 PetscCall(PetscObjectSetName((PetscObject)x,"Approximate Solution")); 42 PetscCall(VecDuplicate(x,&r)); 43 PetscCall(VecDuplicate(x,&F)); 44 PetscCall(VecDuplicate(x,&U)); 45 PetscCall(PetscObjectSetName((PetscObject)U,"Exact Solution")); 46 47 /* create explicit matrix preconditioner */ 48 PetscCall(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 PetscCall(VecSetValues(F,1,&i,&v,INSERT_VALUES)); 54 v = xp*xp*xp; 55 PetscCall(VecSetValues(U,1,&i,&v,INSERT_VALUES)); 56 xp += h; 57 } 58 59 /* Create nonlinear solver */ 60 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 61 PetscCall(SNESSetType(snes,type)); 62 63 /* Set various routines and options */ 64 PetscCall(SNESSetFunction(snes,r,FormFunction,F)); 65 if (user.variant) { 66 /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */ 67 PetscCall(MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J)); 68 PetscCall(MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes)); 69 PetscCall(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 PetscCall(PetscOptionsHasName(NULL,NULL,"-puremf",&puremf)); 73 } else { 74 /* create matrix free matrix for Jacobian */ 75 PetscCall(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 PetscCall(MatMFFDSetFunction(J,OtherFunctionForDifferencing,F)); 79 } 80 81 /* Set various routines and options */ 82 PetscCall(SNESSetJacobian(snes,J,puremf ? J : B,puremf ? FormJacobianNoMatrix : FormJacobian,&user)); 83 PetscCall(SNESSetFromOptions(snes)); 84 85 /* Solve nonlinear system */ 86 PetscCall(FormInitialGuess(snes,x)); 87 PetscCall(SNESSolve(snes,NULL,x)); 88 PetscCall(SNESGetIterationNumber(snes,&its)); 89 PetscCall(PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %" PetscInt_FMT "\n\n",its)); 90 91 /* Free data structures */ 92 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 93 PetscCall(VecDestroy(&U)); PetscCall(VecDestroy(&F)); 94 PetscCall(MatDestroy(&J)); PetscCall(MatDestroy(&B)); 95 PetscCall(SNESDestroy(&snes)); 96 PetscCall(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 PetscCall(VecGetArrayRead(x,&xx)); 108 PetscCall(VecGetArray(f,&ff)); 109 PetscCall(VecGetArrayRead((Vec) dummy,&FF)); 110 PetscCall(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 PetscCall(VecRestoreArrayRead(x,&xx)); 116 PetscCall(VecRestoreArray(f,&ff)); 117 PetscCall(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 PetscCall(SNESGetFunction(snes,NULL,NULL,(void**)&F)); 130 PetscCall(VecGetArrayRead(x,&xx)); 131 PetscCall(VecGetArrayRead(F,&FF)); 132 PetscCall(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 PetscCall(VecRestoreArrayRead(x,&xx)); 142 PetscCall(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 PetscCall(FormFunction(NULL,x,f,dummy)); 155 PetscCall(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 PetscCall(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 PetscCall(VecGetArrayRead(x,&xx)); 181 PetscCall(VecGetSize(x,&n)); 182 d = (PetscReal)(n - 1); d = d*d; 183 184 i = 0; A[0] = 1.0; 185 PetscCall(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 PetscCall(MatSetValues(B,1,&i,3,j,A,INSERT_VALUES)); 190 } 191 i = n-1; A[0] = 1.0; 192 PetscCall(MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES)); 193 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 194 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 195 PetscCall(VecRestoreArrayRead(x,&xx)); 196 197 if (user->variant) { 198 PetscCall(MatMFFDSetBase(jac,x,NULL)); 199 } 200 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 201 PetscCall(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 PetscCall(MatMFFDSetBase(jac,x,NULL)); 211 } 212 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 213 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject)snes,&comm)); 226 PetscCall(PetscFPrintf(comm,stdout,"iter = %" PetscInt_FMT ", SNES Function norm %g \n",its,(double)fnorm)); 227 PetscCall(SNESGetSolution(snes,&x)); 228 PetscCall(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