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 PetscFunctionBeginUser; 112 PetscCall(VecGetArrayRead(x, &xx)); 113 PetscCall(VecGetArray(f, &ff)); 114 PetscCall(VecGetArrayRead((Vec)dummy, &FF)); 115 PetscCall(VecGetSize(x, &n)); 116 d = (PetscReal)(n - 1); 117 d = d * d; 118 ff[0] = xx[0]; 119 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]; 120 ff[n - 1] = xx[n - 1] - 1.0; 121 PetscCall(VecRestoreArrayRead(x, &xx)); 122 PetscCall(VecRestoreArray(f, &ff)); 123 PetscCall(VecRestoreArrayRead((Vec)dummy, &FF)); 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 PetscErrorCode FormFunctioni(void *dummy, PetscInt i, Vec x, PetscScalar *s) 128 { 129 const PetscScalar *xx, *FF; 130 PetscScalar d; 131 PetscInt n; 132 SNES snes = (SNES)dummy; 133 Vec F; 134 135 PetscFunctionBeginUser; 136 PetscCall(SNESGetFunction(snes, NULL, NULL, (void **)&F)); 137 PetscCall(VecGetArrayRead(x, &xx)); 138 PetscCall(VecGetArrayRead(F, &FF)); 139 PetscCall(VecGetSize(x, &n)); 140 d = (PetscReal)(n - 1); 141 d = d * d; 142 if (i == 0) { 143 *s = xx[0]; 144 } else if (i == n - 1) { 145 *s = xx[n - 1] - 1.0; 146 } else { 147 *s = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i]; 148 } 149 PetscCall(VecRestoreArrayRead(x, &xx)); 150 PetscCall(VecRestoreArrayRead(F, &FF)); 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 /* 155 156 Example function that when differenced produces the same matrix-free Jacobian as FormFunction() 157 this is provided to show how a user can provide a different function 158 */ 159 PetscErrorCode OtherFunctionForDifferencing(void *dummy, Vec x, Vec f) 160 { 161 PetscFunctionBeginUser; 162 PetscCall(FormFunction(NULL, x, f, dummy)); 163 PetscCall(VecShift(f, 1.0)); 164 PetscFunctionReturn(PETSC_SUCCESS); 165 } 166 167 /* -------------------- Form initial approximation ----------------- */ 168 169 PetscErrorCode FormInitialGuess(SNES snes, Vec x) 170 { 171 PetscScalar pfive = .50; 172 PetscFunctionBeginUser; 173 PetscCall(VecSet(x, pfive)); 174 PetscFunctionReturn(PETSC_SUCCESS); 175 } 176 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 177 /* Evaluates a matrix that is used to precondition the matrix-free 178 jacobian. In this case, the explicit preconditioner matrix is 179 also EXACTLY the Jacobian. In general, it would be some lower 180 order, simplified apprioximation */ 181 182 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 183 { 184 const PetscScalar *xx; 185 PetscScalar A[3], d; 186 PetscInt i, n, j[3]; 187 AppCtx *user = (AppCtx *)dummy; 188 189 PetscFunctionBeginUser; 190 PetscCall(VecGetArrayRead(x, &xx)); 191 PetscCall(VecGetSize(x, &n)); 192 d = (PetscReal)(n - 1); 193 d = d * d; 194 195 i = 0; 196 A[0] = 1.0; 197 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 198 for (i = 1; i < n - 1; i++) { 199 j[0] = i - 1; 200 j[1] = i; 201 j[2] = i + 1; 202 A[0] = d; 203 A[1] = -2.0 * d + 2.0 * xx[i]; 204 A[2] = d; 205 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 206 } 207 i = n - 1; 208 A[0] = 1.0; 209 PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES)); 210 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 211 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 212 PetscCall(VecRestoreArrayRead(x, &xx)); 213 214 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 215 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 216 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 217 PetscFunctionReturn(PETSC_SUCCESS); 218 } 219 220 PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 221 { 222 AppCtx *user = (AppCtx *)dummy; 223 224 PetscFunctionBeginUser; 225 if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL)); 226 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 227 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 228 PetscFunctionReturn(PETSC_SUCCESS); 229 } 230 231 /* -------------------- User-defined monitor ----------------------- */ 232 233 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy) 234 { 235 MonitorCtx *monP = (MonitorCtx *)dummy; 236 Vec x; 237 MPI_Comm comm; 238 239 PetscFunctionBeginUser; 240 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 241 PetscCall(PetscFPrintf(comm, stdout, "iter = %" PetscInt_FMT ", SNES Function norm %g \n", its, (double)fnorm)); 242 PetscCall(SNESGetSolution(snes, &x)); 243 PetscCall(VecView(x, monP->viewer)); 244 PetscFunctionReturn(PETSC_SUCCESS); 245 } 246 247 /*TEST 248 249 test: 250 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 251 252 test: 253 suffix: 2 254 args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 255 output_file: output/ex7_1.out 256 257 # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning 258 test: 259 suffix: 3 260 args: -variant -pc_type jacobi -snes_view -ksp_monitor 261 262 # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning 263 test: 264 suffix: 4 265 args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor 266 267 TEST*/ 268