1 static char help[] = "Newton methods to solve u'' = f in parallel with periodic boundary conditions.\n\n"; 2 3 /* 4 Compare this example to ex3.c that handles Dirichlet boundary conditions 5 6 Though this is a linear problem it is treated as a nonlinear problem in this example to demonstrate 7 handling periodic boundary conditions for nonlinear problems. 8 9 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 10 Include "petscsnes.h" so that we can use SNES solvers. Note that this 11 file automatically includes: 12 petscsys.h - base PETSc routines petscvec.h - vectors 13 petscmat.h - matrices 14 petscis.h - index sets petscksp.h - Krylov subspace methods 15 petscviewer.h - viewers petscpc.h - preconditioners 16 petscksp.h - linear solvers 17 */ 18 19 #include <petscdm.h> 20 #include <petscdmda.h> 21 #include <petscsnes.h> 22 23 PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 24 PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 25 26 int main(int argc, char **argv) 27 { 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 { 111 DM da = (DM)ctx; 112 PetscScalar *xx, *ff; 113 PetscReal h; 114 PetscInt i, M, xs, xm; 115 Vec xlocal; 116 117 PetscFunctionBeginUser; 118 /* Get local work vector */ 119 PetscCall(DMGetLocalVector(da, &xlocal)); 120 121 /* 122 Scatter ghost points to local vector, using the 2-step process 123 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 124 By placing code between these two statements, computations can 125 be done while messages are in transition. 126 */ 127 PetscCall(DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal)); 128 PetscCall(DMGlobalToLocalEnd(da, x, INSERT_VALUES, xlocal)); 129 130 /* 131 Get pointers to vector data. 132 - The vector xlocal includes ghost point; the vectors x and f do 133 NOT include ghost points. 134 - Using DMDAVecGetArray() allows accessing the values using global ordering 135 */ 136 PetscCall(DMDAVecGetArray(da, xlocal, &xx)); 137 PetscCall(DMDAVecGetArray(da, f, &ff)); 138 139 /* 140 Get local grid boundaries (for 1-dimensional DMDA): 141 xs, xm - starting grid index, width of local grid (no ghost points) 142 */ 143 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 144 PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 145 146 /* 147 Compute function over locally owned part of the grid 148 Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points. 149 */ 150 h = 1.0 / M; 151 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); 152 153 /* 154 Restore vectors 155 */ 156 PetscCall(DMDAVecRestoreArray(da, xlocal, &xx)); 157 PetscCall(DMDAVecRestoreArray(da, f, &ff)); 158 PetscCall(DMRestoreLocalVector(da, &xlocal)); 159 PetscFunctionReturn(PETSC_SUCCESS); 160 } 161 /* ------------------------------------------------------------------- */ 162 /* 163 FormJacobian - Evaluates Jacobian matrix. 164 165 Input Parameters: 166 . snes - the SNES context 167 . x - input vector 168 . dummy - optional user-defined context (not used here) 169 170 Output Parameters: 171 . jac - Jacobian matrix 172 . B - optionally different preconditioning matrix 173 . flag - flag indicating matrix structure 174 */ 175 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *ctx) 176 { 177 PetscScalar *xx, A[3]; 178 PetscInt i, M, xs, xm; 179 DM da = (DM)ctx; 180 MatStencil row, cols[3]; 181 PetscReal h; 182 183 PetscFunctionBeginUser; 184 /* 185 Get pointer to vector data 186 */ 187 PetscCall(DMDAVecGetArrayRead(da, x, &xx)); 188 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 189 190 /* 191 Get range of locally owned matrix 192 */ 193 PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 194 195 PetscCall(MatZeroEntries(jac)); 196 h = 1.0 / M; 197 /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */ 198 for (i = xs; i < xs + xm; i++) { 199 row.i = i; 200 cols[0].i = i - 1; 201 cols[1].i = i; 202 cols[2].i = i + 1; 203 A[0] = A[2] = 1.0 / (h * h); 204 A[1] = -2.0 / (h * h); 205 PetscCall(MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES)); 206 } 207 208 PetscCall(DMDAVecRestoreArrayRead(da, x, &xx)); 209 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 210 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 /*TEST 215 216 test: 217 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 218 requires: !single 219 220 test: 221 suffix: 2 222 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc 223 requires: !single 224 225 test: 226 suffix: 3 227 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc -snes_trdc_use_cauchy false 228 requires: !single 229 230 TEST*/ 231