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 PetscFunctionBeginUser; 35 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 36 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 37 PetscCall(PetscOptionsHasName(NULL,NULL,"-variant",&user.variant)); 38 h = 1.0/(n-1); 39 40 /* Set up data structures */ 41 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,&x)); 42 PetscCall(PetscObjectSetName((PetscObject)x,"Approximate Solution")); 43 PetscCall(VecDuplicate(x,&r)); 44 PetscCall(VecDuplicate(x,&F)); 45 PetscCall(VecDuplicate(x,&U)); 46 PetscCall(PetscObjectSetName((PetscObject)U,"Exact Solution")); 47 48 /* create explicit matrix preconditioner */ 49 PetscCall(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 PetscCall(VecSetValues(F,1,&i,&v,INSERT_VALUES)); 55 v = xp*xp*xp; 56 PetscCall(VecSetValues(U,1,&i,&v,INSERT_VALUES)); 57 xp += h; 58 } 59 60 /* Create nonlinear solver */ 61 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 62 PetscCall(SNESSetType(snes,type)); 63 64 /* Set various routines and options */ 65 PetscCall(SNESSetFunction(snes,r,FormFunction,F)); 66 if (user.variant) { 67 /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */ 68 PetscCall(MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J)); 69 PetscCall(MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes)); 70 PetscCall(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 PetscCall(PetscOptionsHasName(NULL,NULL,"-puremf",&puremf)); 74 } else { 75 /* create matrix free matrix for Jacobian */ 76 PetscCall(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 PetscCall(MatMFFDSetFunction(J,OtherFunctionForDifferencing,F)); 80 } 81 82 /* Set various routines and options */ 83 PetscCall(SNESSetJacobian(snes,J,puremf ? J : B,puremf ? FormJacobianNoMatrix : FormJacobian,&user)); 84 PetscCall(SNESSetFromOptions(snes)); 85 86 /* Solve nonlinear system */ 87 PetscCall(FormInitialGuess(snes,x)); 88 PetscCall(SNESSolve(snes,NULL,x)); 89 PetscCall(SNESGetIterationNumber(snes,&its)); 90 PetscCall(PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %" PetscInt_FMT "\n\n",its)); 91 92 /* Free data structures */ 93 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 94 PetscCall(VecDestroy(&U)); PetscCall(VecDestroy(&F)); 95 PetscCall(MatDestroy(&J)); PetscCall(MatDestroy(&B)); 96 PetscCall(SNESDestroy(&snes)); 97 PetscCall(PetscFinalize()); 98 return 0; 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 PetscCall(VecGetArrayRead(x,&xx)); 109 PetscCall(VecGetArray(f,&ff)); 110 PetscCall(VecGetArrayRead((Vec) dummy,&FF)); 111 PetscCall(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 PetscCall(VecRestoreArrayRead(x,&xx)); 117 PetscCall(VecRestoreArray(f,&ff)); 118 PetscCall(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 PetscCall(SNESGetFunction(snes,NULL,NULL,(void**)&F)); 131 PetscCall(VecGetArrayRead(x,&xx)); 132 PetscCall(VecGetArrayRead(F,&FF)); 133 PetscCall(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 PetscCall(VecRestoreArrayRead(x,&xx)); 143 PetscCall(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 PetscCall(FormFunction(NULL,x,f,dummy)); 156 PetscCall(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 PetscCall(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 PetscCall(VecGetArrayRead(x,&xx)); 182 PetscCall(VecGetSize(x,&n)); 183 d = (PetscReal)(n - 1); d = d*d; 184 185 i = 0; A[0] = 1.0; 186 PetscCall(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 PetscCall(MatSetValues(B,1,&i,3,j,A,INSERT_VALUES)); 191 } 192 i = n-1; A[0] = 1.0; 193 PetscCall(MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES)); 194 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 195 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 196 PetscCall(VecRestoreArrayRead(x,&xx)); 197 198 if (user->variant) PetscCall(MatMFFDSetBase(jac,x,NULL)); 199 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 200 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 201 return 0; 202 } 203 204 PetscErrorCode FormJacobianNoMatrix(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 205 { 206 AppCtx *user = (AppCtx*) dummy; 207 208 if (user->variant) PetscCall(MatMFFDSetBase(jac,x,NULL)); 209 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 210 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 211 return 0; 212 } 213 214 /* -------------------- User-defined monitor ----------------------- */ 215 216 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy) 217 { 218 MonitorCtx *monP = (MonitorCtx*) dummy; 219 Vec x; 220 MPI_Comm comm; 221 222 PetscCall(PetscObjectGetComm((PetscObject)snes,&comm)); 223 PetscCall(PetscFPrintf(comm,stdout,"iter = %" PetscInt_FMT ", SNES Function norm %g \n",its,(double)fnorm)); 224 PetscCall(SNESGetSolution(snes,&x)); 225 PetscCall(VecView(x,monP->viewer)); 226 return 0; 227 } 228 229 /*TEST 230 231 test: 232 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 233 234 test: 235 suffix: 2 236 args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 237 output_file: output/ex7_1.out 238 239 # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning 240 test: 241 suffix: 3 242 args: -variant -pc_type jacobi -snes_view -ksp_monitor 243 244 # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning 245 test: 246 suffix: 4 247 args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor 248 249 TEST*/ 250