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