1 static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\ 2 matrix-free techniques with user-provided explicit preconditioner matrix.\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, (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)); 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 preconditioner matrix 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