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 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL)); 45 46 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47 Create nonlinear solver context 48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 49 50 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 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 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,N,1,1,NULL,&da)); 60 CHKERRQ(DMSetFromOptions(da)); 61 CHKERRQ(DMSetUp(da)); 62 63 /* 64 Extract global and local vectors from DMDA; then duplicate for remaining 65 vectors that are the same types 66 */ 67 CHKERRQ(DMCreateGlobalVector(da,&x)); 68 CHKERRQ(VecDuplicate(x,&r)); 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 CHKERRQ(SNESSetFunction(snes,r,FormFunction,da)); 79 80 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 81 Create matrix data structure; set Jacobian evaluation routine 82 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 83 CHKERRQ(DMCreateMatrix(da,&J)); 84 CHKERRQ(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants)); 85 CHKERRQ(MatSetNullSpace(J,constants)); 86 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian,da)); 87 88 CHKERRQ(SNESSetFromOptions(snes)); 89 CHKERRQ(SNESSolve(snes,NULL,x)); 90 91 CHKERRQ(VecDestroy(&x)); 92 CHKERRQ(VecDestroy(&r)); 93 CHKERRQ(MatDestroy(&J)); 94 CHKERRQ(MatNullSpaceDestroy(&constants)); 95 CHKERRQ(SNESDestroy(&snes)); 96 CHKERRQ(DMDestroy(&da)); 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 PetscInt i,M,xs,xm; 122 Vec xlocal; 123 124 PetscFunctionBeginUser; 125 /* Get local work vector */ 126 CHKERRQ(DMGetLocalVector(da,&xlocal)); 127 128 /* 129 Scatter ghost points to local vector, using the 2-step process 130 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 131 By placing code between these two statements, computations can 132 be done while messages are in transition. 133 */ 134 CHKERRQ(DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal)); 135 CHKERRQ(DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal)); 136 137 /* 138 Get pointers to vector data. 139 - The vector xlocal includes ghost point; the vectors x and f do 140 NOT include ghost points. 141 - Using DMDAVecGetArray() allows accessing the values using global ordering 142 */ 143 CHKERRQ(DMDAVecGetArray(da,xlocal,&xx)); 144 CHKERRQ(DMDAVecGetArray(da,f,&ff)); 145 146 /* 147 Get local grid boundaries (for 1-dimensional DMDA): 148 xs, xm - starting grid index, width of local grid (no ghost points) 149 */ 150 CHKERRQ(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 151 CHKERRQ(DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 152 153 /* 154 Compute function over locally owned part of the grid 155 Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points. 156 */ 157 h = 1.0/M; 158 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); 159 160 /* 161 Restore vectors 162 */ 163 CHKERRQ(DMDAVecRestoreArray(da,xlocal,&xx)); 164 CHKERRQ(DMDAVecRestoreArray(da,f,&ff)); 165 CHKERRQ(DMRestoreLocalVector(da,&xlocal)); 166 PetscFunctionReturn(0); 167 } 168 /* ------------------------------------------------------------------- */ 169 /* 170 FormJacobian - Evaluates Jacobian matrix. 171 172 Input Parameters: 173 . snes - the SNES context 174 . x - input vector 175 . dummy - optional user-defined context (not used here) 176 177 Output Parameters: 178 . jac - Jacobian matrix 179 . B - optionally different preconditioning matrix 180 . flag - flag indicating matrix structure 181 */ 182 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx) 183 { 184 PetscScalar *xx,A[3]; 185 PetscInt i,M,xs,xm; 186 DM da = (DM) ctx; 187 MatStencil row,cols[3]; 188 PetscReal h; 189 190 PetscFunctionBeginUser; 191 /* 192 Get pointer to vector data 193 */ 194 CHKERRQ(DMDAVecGetArrayRead(da,x,&xx)); 195 CHKERRQ(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 196 197 /* 198 Get range of locally owned matrix 199 */ 200 CHKERRQ(DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 201 202 CHKERRQ(MatZeroEntries(jac)); 203 h = 1.0/M; 204 /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */ 205 for (i=xs; i<xs+xm; i++) { 206 row.i = i; 207 cols[0].i = i - 1; 208 cols[1].i = i; 209 cols[2].i = i + 1; 210 A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h); 211 CHKERRQ(MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES)); 212 } 213 214 CHKERRQ(DMDAVecRestoreArrayRead(da,x,&xx)); 215 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 216 CHKERRQ(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 217 PetscFunctionReturn(0); 218 } 219 220 /*TEST 221 222 test: 223 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 224 requires: !single 225 226 test: 227 suffix: 2 228 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc 229 requires: !single 230 231 test: 232 suffix: 3 233 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc -snes_trdc_use_cauchy false 234 requires: !single 235 236 TEST*/ 237