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)); 94 PetscCall(VecDestroy(&r)); 95 PetscCall(VecDestroy(&U)); 96 PetscCall(VecDestroy(&F)); 97 PetscCall(MatDestroy(&J)); 98 PetscCall(MatDestroy(&B)); 99 PetscCall(SNESDestroy(&snes)); 100 PetscCall(PetscFinalize()); 101 return 0; 102 } 103 /* -------------------- Evaluate Function F(x) --------------------- */ 104 105 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *dummy) 106 { 107 const PetscScalar *xx, *FF; 108 PetscScalar *ff, d; 109 PetscInt i, n; 110 111 PetscCall(VecGetArrayRead(x, &xx)); 112 PetscCall(VecGetArray(f, &ff)); 113 PetscCall(VecGetArrayRead((Vec)dummy, &FF)); 114 PetscCall(VecGetSize(x, &n)); 115 d = (PetscReal)(n - 1); 116 d = d * d; 117 ff[0] = xx[0]; 118 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]; 119 ff[n - 1] = xx[n - 1] - 1.0; 120 PetscCall(VecRestoreArrayRead(x, &xx)); 121 PetscCall(VecRestoreArray(f, &ff)); 122 PetscCall(VecRestoreArrayRead((Vec)dummy, &FF)); 123 return 0; 124 } 125 126 PetscErrorCode FormFunctioni(void *dummy, PetscInt i, Vec x, PetscScalar *s) 127 { 128 const PetscScalar *xx, *FF; 129 PetscScalar d; 130 PetscInt n; 131 SNES snes = (SNES)dummy; 132 Vec F; 133 134 PetscCall(SNESGetFunction(snes, NULL, NULL, (void **)&F)); 135 PetscCall(VecGetArrayRead(x, &xx)); 136 PetscCall(VecGetArrayRead(F, &FF)); 137 PetscCall(VecGetSize(x, &n)); 138 d = (PetscReal)(n - 1); 139 d = d * d; 140 if (i == 0) { 141 *s = xx[0]; 142 } else if (i == n - 1) { 143 *s = xx[n - 1] - 1.0; 144 } else { 145 *s = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i]; 146 } 147 PetscCall(VecRestoreArrayRead(x, &xx)); 148 PetscCall(VecRestoreArrayRead(F, &FF)); 149 return 0; 150 } 151 152 /* 153 154 Example function that when differenced produces the same matrix free Jacobian as FormFunction() 155 this is provided to show how a user can provide a different function 156 */ 157 PetscErrorCode OtherFunctionForDifferencing(void *dummy, Vec x, Vec f) 158 { 159 PetscCall(FormFunction(NULL, x, f, dummy)); 160 PetscCall(VecShift(f, 1.0)); 161 return 0; 162 } 163 164 /* -------------------- Form initial approximation ----------------- */ 165 166 PetscErrorCode FormInitialGuess(SNES snes, Vec x) 167 { 168 PetscScalar pfive = .50; 169 PetscCall(VecSet(x, pfive)); 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 AppCtx *user = (AppCtx *)dummy; 184 185 PetscCall(VecGetArrayRead(x, &xx)); 186 PetscCall(VecGetSize(x, &n)); 187 d = (PetscReal)(n - 1); 188 d = d * d; 189 190 i = 0; 191 A[0] = 1.0; 192 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 193 for (i = 1; i < n - 1; i++) { 194 j[0] = i - 1; 195 j[1] = i; 196 j[2] = i + 1; 197 A[0] = d; 198 A[1] = -2.0 * d + 2.0 * xx[i]; 199 A[2] = d; 200 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 201 } 202 i = n - 1; 203 A[0] = 1.0; 204 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 205 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 206 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 207 PetscCall(VecRestoreArrayRead(x, &xx)); 208 209 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 210 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 211 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 212 return 0; 213 } 214 215 PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 216 { 217 AppCtx *user = (AppCtx *)dummy; 218 219 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 220 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 221 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 222 return 0; 223 } 224 225 /* -------------------- User-defined monitor ----------------------- */ 226 227 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy) 228 { 229 MonitorCtx *monP = (MonitorCtx *)dummy; 230 Vec x; 231 MPI_Comm comm; 232 233 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 234 PetscCall(PetscFPrintf(comm, stdout, "iter = %" PetscInt_FMT ", SNES Function norm %g \n", its, (double)fnorm)); 235 PetscCall(SNESGetSolution(snes, &x)); 236 PetscCall(VecView(x, monP->viewer)); 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