1 static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\ 2 matrix-free techniques with user-provided explicit matrix for computing the preconditioner.\n\n"; 3 4 #include <petscsnes.h> 5 6 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 7 extern PetscErrorCode FormJacobianNoMatrix(SNES, Vec, Mat, Mat, void *); 8 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 9 extern PetscErrorCode FormFunctioni(void *, PetscInt, Vec, PetscScalar *); 10 extern PetscErrorCode OtherFunctionForDifferencing(void *, Vec, Vec); 11 extern PetscErrorCode FormInitialGuess(SNES, Vec); 12 extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *); 13 14 typedef struct { 15 PetscViewer viewer; 16 } MonitorCtx; 17 18 typedef struct { 19 PetscBool variant; 20 } AppCtx; 21 22 int main(int argc, char **argv) 23 { 24 SNES snes; /* SNES context */ 25 SNESType type = SNESNEWTONLS; /* default nonlinear solution method */ 26 Vec x, r, F, U; /* vectors */ 27 Mat J, B; /* Jacobian matrix-free, explicit preconditioner */ 28 AppCtx user; /* user-defined work context */ 29 PetscScalar h, xp = 0.0, v; 30 PetscInt its, n = 5, i; 31 PetscBool puremf = PETSC_FALSE; 32 33 PetscFunctionBeginUser; 34 PetscCall(PetscInitialize(&argc, &argv, NULL, 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)); 93 PetscCall(VecDestroy(&r)); 94 PetscCall(VecDestroy(&U)); 95 PetscCall(VecDestroy(&F)); 96 PetscCall(MatDestroy(&J)); 97 PetscCall(MatDestroy(&B)); 98 PetscCall(SNESDestroy(&snes)); 99 PetscCall(PetscFinalize()); 100 return 0; 101 } 102 /* -------------------- Evaluate Function F(x) --------------------- */ 103 104 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *dummy) 105 { 106 const PetscScalar *xx, *FF; 107 PetscScalar *ff, d; 108 PetscInt i, n; 109 110 PetscFunctionBeginUser; 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 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 135 PetscCall(SNESGetFunction(snes, NULL, NULL, (void **)&F)); 136 PetscCall(VecGetArrayRead(x, &xx)); 137 PetscCall(VecGetArrayRead(F, &FF)); 138 PetscCall(VecGetSize(x, &n)); 139 d = (PetscReal)(n - 1); 140 d = d * d; 141 if (i == 0) { 142 *s = xx[0]; 143 } else if (i == n - 1) { 144 *s = xx[n - 1] - 1.0; 145 } else { 146 *s = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i]; 147 } 148 PetscCall(VecRestoreArrayRead(x, &xx)); 149 PetscCall(VecRestoreArrayRead(F, &FF)); 150 PetscFunctionReturn(PETSC_SUCCESS); 151 } 152 153 /* 154 155 Example function that when differenced produces the same matrix-free Jacobian as FormFunction() 156 this is provided to show how a user can provide a different function 157 */ 158 PetscErrorCode OtherFunctionForDifferencing(void *dummy, Vec x, Vec f) 159 { 160 PetscFunctionBeginUser; 161 PetscCall(FormFunction(NULL, x, f, dummy)); 162 PetscCall(VecShift(f, 1.0)); 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 /* -------------------- Form initial approximation ----------------- */ 167 168 PetscErrorCode FormInitialGuess(SNES snes, Vec x) 169 { 170 PetscFunctionBeginUser; 171 PetscCall(VecSet(x, 0.5)); 172 PetscFunctionReturn(PETSC_SUCCESS); 173 } 174 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 175 /* Evaluates a matrix that is used to precondition the matrix-free 176 jacobian. In this case, the explicit matrix used to compute the preconditioner is 177 also EXACTLY the Jacobian. In general, it would be some lower 178 order, simplified apprioximation */ 179 180 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 181 { 182 const PetscScalar *xx; 183 PetscScalar A[3], d; 184 PetscInt i, n, j[3]; 185 AppCtx *user = (AppCtx *)dummy; 186 187 PetscFunctionBeginUser; 188 PetscCall(VecGetArrayRead(x, &xx)); 189 PetscCall(VecGetSize(x, &n)); 190 d = (PetscReal)(n - 1); 191 d = d * d; 192 193 i = 0; 194 A[0] = 1.0; 195 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 196 for (i = 1; i < n - 1; i++) { 197 j[0] = i - 1; 198 j[1] = i; 199 j[2] = i + 1; 200 A[0] = d; 201 A[1] = -2.0 * d + 2.0 * xx[i]; 202 A[2] = d; 203 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 204 } 205 i = n - 1; 206 A[0] = 1.0; 207 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 208 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 209 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 210 PetscCall(VecRestoreArrayRead(x, &xx)); 211 212 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 213 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 214 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 215 PetscFunctionReturn(PETSC_SUCCESS); 216 } 217 218 PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 219 { 220 AppCtx *user = (AppCtx *)dummy; 221 222 PetscFunctionBeginUser; 223 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 224 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 225 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 226 PetscFunctionReturn(PETSC_SUCCESS); 227 } 228 229 /* -------------------- User-defined monitor ----------------------- */ 230 231 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy) 232 { 233 MonitorCtx *monP = (MonitorCtx *)dummy; 234 Vec x; 235 MPI_Comm comm; 236 237 PetscFunctionBeginUser; 238 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 239 PetscCall(PetscFPrintf(comm, stdout, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm)); 240 PetscCall(SNESGetSolution(snes, &x)); 241 PetscCall(VecView(x, monP->viewer)); 242 PetscFunctionReturn(PETSC_SUCCESS); 243 } 244 245 /*TEST 246 247 test: 248 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 249 250 test: 251 suffix: 2 252 args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 253 output_file: output/ex7_1.out 254 255 # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning 256 test: 257 requires: !single 258 suffix: 3 259 args: -variant -pc_type jacobi -snes_view -ksp_monitor 260 261 # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning 262 test: 263 suffix: 4 264 args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor 265 266 TEST*/ 267