1 static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n"; 2 3 /* 4 Include "petscsnes.h" so that we can use SNES solvers. Note that this 5 file automatically includes: 6 petscsys.h - base PETSc routines petscvec.h - vectors 7 petscmat.h - matrices 8 petscis.h - index sets petscksp.h - Krylov subspace methods 9 petscviewer.h - viewers petscpc.h - preconditioners 10 petscksp.h - linear solvers 11 */ 12 #include <petscsnes.h> 13 14 extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *); 15 extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *); 16 17 int main(int argc, char **argv) 18 { 19 SNES snes; /* nonlinear solver context */ 20 Vec x; /* solution vector */ 21 Mat J; /* Jacobian matrix */ 22 PetscInt its; 23 PetscScalar *xx; 24 SNESConvergedReason reason; 25 PetscBool test_ghost = PETSC_FALSE; 26 27 PetscFunctionBeginUser; 28 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 29 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_ghost", &test_ghost, NULL)); 30 31 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 32 Create nonlinear solver context 33 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 34 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 35 36 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 37 Create matrix and vector data structures; set corresponding routines 38 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 39 /* 40 Create vectors for solution and nonlinear function 41 */ 42 if (test_ghost) { 43 PetscInt gIdx[] = {0, 1}; 44 PetscInt nghost = 2; 45 46 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nghost", &nghost, NULL)); 47 PetscCall(VecCreateGhost(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PetscMin(nghost, 2), gIdx, &x)); 48 } else { 49 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 50 PetscCall(VecSetSizes(x, 2, PETSC_DECIDE)); 51 PetscCall(VecSetFromOptions(x)); 52 } 53 54 /* 55 Create Jacobian matrix data structure 56 */ 57 PetscCall(MatCreate(PETSC_COMM_WORLD, &J)); 58 PetscCall(MatSetSizes(J, 2, 2, PETSC_DECIDE, PETSC_DECIDE)); 59 PetscCall(MatSetFromOptions(J)); 60 PetscCall(MatSetUp(J)); 61 62 /* 63 Set function evaluation routine and vector. 64 */ 65 PetscCall(SNESSetFunction(snes, NULL, FormFunction1, &test_ghost)); 66 67 /* 68 Set Jacobian matrix data structure and Jacobian evaluation routine 69 */ 70 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, &test_ghost)); 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Customize nonlinear solver; set runtime options 74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 75 PetscCall(SNESSetFromOptions(snes)); 76 77 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 78 Evaluate initial guess; then solve nonlinear system 79 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 80 PetscCall(VecGetArray(x, &xx)); 81 xx[0] = -1.2; 82 xx[1] = 1.0; 83 PetscCall(VecRestoreArray(x, &xx)); 84 85 /* 86 Note: The user should initialize the vector, x, with the initial guess 87 for the nonlinear solver prior to calling SNESSolve(). In particular, 88 to employ an initial guess of zero, the user should explicitly set 89 this vector to zero by calling VecSet(). 90 */ 91 92 PetscCall(SNESSolve(snes, NULL, x)); 93 PetscCall(VecViewFromOptions(x, NULL, "-sol_view")); 94 PetscCall(SNESGetIterationNumber(snes, &its)); 95 PetscCall(SNESGetConvergedReason(snes, &reason)); 96 /* 97 Some systems computes a residual that is identically zero, thus converging 98 due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only 99 report the reason if the iteration did not converge so that the tests are 100 reproducible. 101 */ 102 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s number of SNES iterations = %" PetscInt_FMT "\n", reason > 0 ? "CONVERGED" : SNESConvergedReasons[reason], its)); 103 104 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 105 Free work space. All PETSc objects should be destroyed when they 106 are no longer needed. 107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 108 109 PetscCall(VecDestroy(&x)); 110 PetscCall(MatDestroy(&J)); 111 PetscCall(SNESDestroy(&snes)); 112 PetscCall(PetscFinalize()); 113 return 0; 114 } 115 116 PetscErrorCode VecCheckGhosted(Vec X, PetscBool test_rev) 117 { 118 PetscFunctionBeginUser; 119 PetscCall(VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD)); 120 PetscCall(VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD)); 121 if (test_rev) { 122 PetscCall(VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_REVERSE)); 123 PetscCall(VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_REVERSE)); 124 } 125 PetscFunctionReturn(PETSC_SUCCESS); 126 } 127 128 PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx) 129 { 130 PetscScalar *ff; 131 const PetscScalar *xx; 132 133 PetscFunctionBeginUser; 134 if (*(PetscBool *)ctx) { 135 PetscCall(VecCheckGhosted(x, PETSC_FALSE)); 136 PetscCall(VecCheckGhosted(f, PETSC_TRUE)); 137 } 138 139 /* 140 Get pointers to vector data. 141 - For default PETSc vectors, VecGetArray() returns a pointer to 142 the data array. Otherwise, the routine is implementation dependent. 143 - You MUST call VecRestoreArray() when you no longer need access to 144 the array. 145 */ 146 PetscCall(VecGetArrayRead(x, &xx)); 147 PetscCall(VecGetArray(f, &ff)); 148 149 /* Compute function */ 150 ff[0] = -2.0 + 2.0 * xx[0] + 400.0 * xx[0] * xx[0] * xx[0] - 400.0 * xx[0] * xx[1]; 151 ff[1] = -200.0 * xx[0] * xx[0] + 200.0 * xx[1]; 152 153 /* Restore vectors */ 154 PetscCall(VecRestoreArrayRead(x, &xx)); 155 PetscCall(VecRestoreArray(f, &ff)); 156 PetscFunctionReturn(PETSC_SUCCESS); 157 } 158 159 PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *ctx) 160 { 161 const PetscScalar *xx; 162 PetscScalar A[4]; 163 PetscInt idx[2]; 164 PetscMPIInt rank; 165 166 PetscFunctionBeginUser; 167 if (*(PetscBool *)ctx) PetscCall(VecCheckGhosted(x, PETSC_FALSE)); 168 /* 169 Get pointer to vector data 170 */ 171 PetscCall(VecGetArrayRead(x, &xx)); 172 173 /* 174 Compute Jacobian entries and insert into matrix. 175 - Since this is such a small problem, we set all entries for 176 the matrix at once. 177 */ 178 A[0] = 2.0 + 1200.0 * xx[0] * xx[0] - 400.0 * xx[1]; 179 A[1] = -400.0 * xx[0]; 180 A[2] = -400.0 * xx[0]; 181 A[3] = 200; 182 183 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)x), &rank)); 184 idx[0] = 2 * rank; 185 idx[1] = 2 * rank + 1; 186 187 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 188 189 /* 190 Restore vector 191 */ 192 PetscCall(VecRestoreArrayRead(x, &xx)); 193 194 /* 195 Assemble matrix 196 */ 197 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 198 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 199 if (jac != B) { 200 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 201 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 202 } 203 PetscFunctionReturn(PETSC_SUCCESS); 204 } 205 206 /*TEST 207 208 test: 209 suffix: 1 210 args: -snes_monitor_short -snes_max_it 1000 -sol_view 211 requires: !single 212 213 test: 214 suffix: 2 215 args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false -sol_view 216 requires: !single 217 218 test: 219 suffix: ghosts 220 nsize: {{1 2}} 221 args: -snes_max_it 4 -snes_type {{newtontr newtonls}} -nghost {{0 1 2}} -test_ghost 222 requires: !single 223 224 TEST*/ 225