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 PetscScalar pfive = .50; 171 PetscFunctionBeginUser; 172 PetscCall(VecSet(x, pfive)); 173 PetscFunctionReturn(PETSC_SUCCESS); 174 } 175 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 176 /* Evaluates a matrix that is used to precondition the matrix-free 177 jacobian. In this case, the explicit preconditioner matrix is 178 also EXACTLY the Jacobian. In general, it would be some lower 179 order, simplified apprioximation */ 180 181 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 182 { 183 const PetscScalar *xx; 184 PetscScalar A[3], d; 185 PetscInt i, n, j[3]; 186 AppCtx *user = (AppCtx *)dummy; 187 188 PetscFunctionBeginUser; 189 PetscCall(VecGetArrayRead(x, &xx)); 190 PetscCall(VecGetSize(x, &n)); 191 d = (PetscReal)(n - 1); 192 d = d * d; 193 194 i = 0; 195 A[0] = 1.0; 196 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 197 for (i = 1; i < n - 1; i++) { 198 j[0] = i - 1; 199 j[1] = i; 200 j[2] = i + 1; 201 A[0] = d; 202 A[1] = -2.0 * d + 2.0 * xx[i]; 203 A[2] = d; 204 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 205 } 206 i = n - 1; 207 A[0] = 1.0; 208 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 209 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 210 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 211 PetscCall(VecRestoreArrayRead(x, &xx)); 212 213 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 214 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 215 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 216 PetscFunctionReturn(PETSC_SUCCESS); 217 } 218 219 PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 220 { 221 AppCtx *user = (AppCtx *)dummy; 222 223 PetscFunctionBeginUser; 224 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 225 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 226 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 227 PetscFunctionReturn(PETSC_SUCCESS); 228 } 229 230 /* -------------------- User-defined monitor ----------------------- */ 231 232 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy) 233 { 234 MonitorCtx *monP = (MonitorCtx *)dummy; 235 Vec x; 236 MPI_Comm comm; 237 238 PetscFunctionBeginUser; 239 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 240 PetscCall(PetscFPrintf(comm, stdout, "iter = %" PetscInt_FMT ", SNES Function norm %g \n", its, (double)fnorm)); 241 PetscCall(SNESGetSolution(snes, &x)); 242 PetscCall(VecView(x, monP->viewer)); 243 PetscFunctionReturn(PETSC_SUCCESS); 244 } 245 246 /*TEST 247 248 test: 249 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 250 251 test: 252 suffix: 2 253 args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 254 output_file: output/ex7_1.out 255 256 # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning 257 test: 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