1 2 static char help[] = "Newton methods to solve u'' = f in parallel with periodic boundary conditions.\n\n"; 3 4 /* 5 Compare this example to ex3.c that handles Dirichlet boundary conditions 6 7 Though this is a linear problem it is treated as a nonlinear problem in this example to demonstrate 8 handling periodic boundary conditions for nonlinear problems. 9 10 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 11 Include "petscsnes.h" so that we can use SNES solvers. Note that this 12 file automatically includes: 13 petscsys.h - base PETSc routines petscvec.h - vectors 14 petscmat.h - matrices 15 petscis.h - index sets petscksp.h - Krylov subspace methods 16 petscviewer.h - viewers petscpc.h - preconditioners 17 petscksp.h - linear solvers 18 */ 19 20 #include <petscdm.h> 21 #include <petscdmda.h> 22 #include <petscsnes.h> 23 24 PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 25 PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 26 27 int main(int argc, char **argv) { 28 SNES snes; /* SNES context */ 29 Mat J; /* Jacobian matrix */ 30 DM da; 31 Vec x, r; /* vectors */ 32 PetscInt N = 5; 33 MatNullSpace constants; 34 35 PetscFunctionBeginUser; 36 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 37 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL)); 38 39 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 40 Create nonlinear solver context 41 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 42 43 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 44 45 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 46 Create vector data structures; set function evaluation routine 47 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 48 49 /* 50 Create distributed array (DMDA) to manage parallel grid and vectors 51 */ 52 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, N, 1, 1, NULL, &da)); 53 PetscCall(DMSetFromOptions(da)); 54 PetscCall(DMSetUp(da)); 55 56 /* 57 Extract global and local vectors from DMDA; then duplicate for remaining 58 vectors that are the same types 59 */ 60 PetscCall(DMCreateGlobalVector(da, &x)); 61 PetscCall(VecDuplicate(x, &r)); 62 63 /* 64 Set function evaluation routine and vector. Whenever the nonlinear 65 solver needs to compute the nonlinear function, it will call this 66 routine. 67 - Note that the final routine argument is the user-defined 68 context that provides application-specific data for the 69 function evaluation routine. 70 */ 71 PetscCall(SNESSetFunction(snes, r, FormFunction, da)); 72 73 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 74 Create matrix data structure; set Jacobian evaluation routine 75 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 76 PetscCall(DMCreateMatrix(da, &J)); 77 PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &constants)); 78 PetscCall(MatSetNullSpace(J, constants)); 79 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, da)); 80 81 PetscCall(SNESSetFromOptions(snes)); 82 PetscCall(SNESSolve(snes, NULL, x)); 83 84 PetscCall(VecDestroy(&x)); 85 PetscCall(VecDestroy(&r)); 86 PetscCall(MatDestroy(&J)); 87 PetscCall(MatNullSpaceDestroy(&constants)); 88 PetscCall(SNESDestroy(&snes)); 89 PetscCall(DMDestroy(&da)); 90 PetscCall(PetscFinalize()); 91 return 0; 92 } 93 94 /* 95 FormFunction - Evaluates nonlinear function, F(x). 96 97 Input Parameters: 98 . snes - the SNES context 99 . x - input vector 100 . ctx - optional user-defined context, as set by SNESSetFunction() 101 102 Output Parameter: 103 . f - function vector 104 105 Note: 106 The user-defined context can contain any application-specific 107 data needed for the function evaluation. 108 */ 109 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) { 110 DM da = (DM)ctx; 111 PetscScalar *xx, *ff; 112 PetscReal h; 113 PetscInt i, M, xs, xm; 114 Vec xlocal; 115 116 PetscFunctionBeginUser; 117 /* Get local work vector */ 118 PetscCall(DMGetLocalVector(da, &xlocal)); 119 120 /* 121 Scatter ghost points to local vector, using the 2-step process 122 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 123 By placing code between these two statements, computations can 124 be done while messages are in transition. 125 */ 126 PetscCall(DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal)); 127 PetscCall(DMGlobalToLocalEnd(da, x, INSERT_VALUES, xlocal)); 128 129 /* 130 Get pointers to vector data. 131 - The vector xlocal includes ghost point; the vectors x and f do 132 NOT include ghost points. 133 - Using DMDAVecGetArray() allows accessing the values using global ordering 134 */ 135 PetscCall(DMDAVecGetArray(da, xlocal, &xx)); 136 PetscCall(DMDAVecGetArray(da, f, &ff)); 137 138 /* 139 Get local grid boundaries (for 1-dimensional DMDA): 140 xs, xm - starting grid index, width of local grid (no ghost points) 141 */ 142 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 143 PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 144 145 /* 146 Compute function over locally owned part of the grid 147 Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points. 148 */ 149 h = 1.0 / M; 150 for (i = xs; i < xs + xm; i++) ff[i] = (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) / (h * h) - PetscSinReal(2.0 * PETSC_PI * i * h); 151 152 /* 153 Restore vectors 154 */ 155 PetscCall(DMDAVecRestoreArray(da, xlocal, &xx)); 156 PetscCall(DMDAVecRestoreArray(da, f, &ff)); 157 PetscCall(DMRestoreLocalVector(da, &xlocal)); 158 PetscFunctionReturn(0); 159 } 160 /* ------------------------------------------------------------------- */ 161 /* 162 FormJacobian - Evaluates Jacobian matrix. 163 164 Input Parameters: 165 . snes - the SNES context 166 . x - input vector 167 . dummy - optional user-defined context (not used here) 168 169 Output Parameters: 170 . jac - Jacobian matrix 171 . B - optionally different preconditioning matrix 172 . flag - flag indicating matrix structure 173 */ 174 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *ctx) { 175 PetscScalar *xx, A[3]; 176 PetscInt i, M, xs, xm; 177 DM da = (DM)ctx; 178 MatStencil row, cols[3]; 179 PetscReal h; 180 181 PetscFunctionBeginUser; 182 /* 183 Get pointer to vector data 184 */ 185 PetscCall(DMDAVecGetArrayRead(da, x, &xx)); 186 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 187 188 /* 189 Get range of locally owned matrix 190 */ 191 PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 192 193 PetscCall(MatZeroEntries(jac)); 194 h = 1.0 / M; 195 /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */ 196 for (i = xs; i < xs + xm; i++) { 197 row.i = i; 198 cols[0].i = i - 1; 199 cols[1].i = i; 200 cols[2].i = i + 1; 201 A[0] = A[2] = 1.0 / (h * h); 202 A[1] = -2.0 / (h * h); 203 PetscCall(MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES)); 204 } 205 206 PetscCall(DMDAVecRestoreArrayRead(da, x, &xx)); 207 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 208 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 209 PetscFunctionReturn(0); 210 } 211 212 /*TEST 213 214 test: 215 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 216 requires: !single 217 218 test: 219 suffix: 2 220 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc 221 requires: !single 222 223 test: 224 suffix: 3 225 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc -snes_trdc_use_cauchy false 226 requires: !single 227 228 TEST*/ 229