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