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