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(0); 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); 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(0); 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