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