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 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 const PetscScalar *xx, *FF; 106 PetscScalar *ff, d; 107 PetscInt i, n; 108 109 PetscCall(VecGetArrayRead(x, &xx)); 110 PetscCall(VecGetArray(f, &ff)); 111 PetscCall(VecGetArrayRead((Vec)dummy, &FF)); 112 PetscCall(VecGetSize(x, &n)); 113 d = (PetscReal)(n - 1); 114 d = d * d; 115 ff[0] = xx[0]; 116 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]; 117 ff[n - 1] = xx[n - 1] - 1.0; 118 PetscCall(VecRestoreArrayRead(x, &xx)); 119 PetscCall(VecRestoreArray(f, &ff)); 120 PetscCall(VecRestoreArrayRead((Vec)dummy, &FF)); 121 return 0; 122 } 123 124 PetscErrorCode FormFunctioni(void *dummy, PetscInt i, Vec x, PetscScalar *s) { 125 const PetscScalar *xx, *FF; 126 PetscScalar d; 127 PetscInt n; 128 SNES snes = (SNES)dummy; 129 Vec F; 130 131 PetscCall(SNESGetFunction(snes, NULL, NULL, (void **)&F)); 132 PetscCall(VecGetArrayRead(x, &xx)); 133 PetscCall(VecGetArrayRead(F, &FF)); 134 PetscCall(VecGetSize(x, &n)); 135 d = (PetscReal)(n - 1); 136 d = d * d; 137 if (i == 0) { 138 *s = xx[0]; 139 } else if (i == n - 1) { 140 *s = xx[n - 1] - 1.0; 141 } else { 142 *s = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i]; 143 } 144 PetscCall(VecRestoreArrayRead(x, &xx)); 145 PetscCall(VecRestoreArrayRead(F, &FF)); 146 return 0; 147 } 148 149 /* 150 151 Example function that when differenced produces the same matrix free Jacobian as FormFunction() 152 this is provided to show how a user can provide a different function 153 */ 154 PetscErrorCode OtherFunctionForDifferencing(void *dummy, Vec x, Vec f) { 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 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 const PetscScalar *xx; 175 PetscScalar A[3], d; 176 PetscInt i, n, j[3]; 177 AppCtx *user = (AppCtx *)dummy; 178 179 PetscCall(VecGetArrayRead(x, &xx)); 180 PetscCall(VecGetSize(x, &n)); 181 d = (PetscReal)(n - 1); 182 d = d * d; 183 184 i = 0; 185 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; 189 j[1] = i; 190 j[2] = i + 1; 191 A[0] = d; 192 A[1] = -2.0 * d + 2.0 * xx[i]; 193 A[2] = d; 194 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 195 } 196 i = n - 1; 197 A[0] = 1.0; 198 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 199 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 200 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 201 PetscCall(VecRestoreArrayRead(x, &xx)); 202 203 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 204 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 205 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 206 return 0; 207 } 208 209 PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy) { 210 AppCtx *user = (AppCtx *)dummy; 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 return 0; 216 } 217 218 /* -------------------- User-defined monitor ----------------------- */ 219 220 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy) { 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