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