Actual source code: ex3.c
petsc-3.8.4 2018-03-24
2: static char help[] = "Tests MATNEST containing MATSHELL\n\n";
5: /* ------------------------------------------------------------------------
7: This tests if MATNEST properly detects it cannot perform MatMultTranspose()
8: when one of its shell matrices cannot perform MatMultTransposeAdd()
10: The Jacobian is intentionally wrong to force a DIVERGED_LINE_SEARCH which
11: forces a call to SNESNEWTONLSCheckLocalMin_Private() to confirm that
12: MatHasOperation(mat,MATOP_MULT_TRANSPOSE,&flg) returns that the nest matrix
13: cannot apply the transpose.
15: The original bug that was fixed by providing a MatHasOperation_Nest()
16: was reported by Manav Bhatia <bhatiamanav@gmail.com>.
18: ------------------------------------------------------------------------- */
21: #include <petscsnes.h>
23: /*
24: User-defined application context - contains data needed by the
25: application-provided call-back routines, FormJacobian() and
26: FormFunction().
27: */
28: typedef struct {
29: PetscReal param; /* test problem parameter */
30: PetscInt mx; /* Discretization in x-direction */
31: PetscInt my; /* Discretization in y-direction */
32: } AppCtx;
34: PetscErrorCode mymatmult(Mat shell,Vec x,Vec y)
35: {
36: Mat J;
40: MatShellGetContext(shell,&J);
41: MatMult(J,x,y);
42: return(0);
43: }
45: PetscErrorCode mymatdestroy(Mat shell)
46: {
47: Mat J;
51: MatShellGetContext(shell,&J);
52: MatDestroy(&J);
53: return(0);
54: }
56: /*
57: User-defined routines
58: */
59: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
60: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
61: extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
63: int main(int argc,char **argv)
64: {
65: SNES snes; /* nonlinear solver context */
66: Vec x,r; /* solution, residual vectors */
67: Mat J,nest,shell; /* Jacobian matrix */
68: AppCtx user; /* user-defined application context */
70: PetscInt its,N;
71: PetscMPIInt size;
72: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
73: IS is;
75: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
76: MPI_Comm_size(PETSC_COMM_WORLD,&size);
77: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
79: /*
80: Initialize problem parameters
81: */
82: user.mx = 4; user.my = 4; user.param = 6.0;
83: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
84: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
85: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
86: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
87: N = user.mx*user.my;
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Create nonlinear solver context
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: SNESCreate(PETSC_COMM_WORLD,&snes);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Create vector data structures; set function evaluation routine
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: VecCreate(PETSC_COMM_WORLD,&x);
100: VecSetSizes(x,PETSC_DECIDE,N);
101: VecSetFromOptions(x);
102: VecDuplicate(x,&r);
104: /*
105: Set function evaluation routine and vector. Whenever the nonlinear
106: solver needs to evaluate the nonlinear function, it will call this
107: routine.
108: - Note that the final routine argument is the user-defined
109: context that provides application-specific data for the
110: function evaluation routine.
111: */
112: SNESSetFunction(snes,r,FormFunction,(void*)&user);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Create matrix data structure; set Jacobian evaluation routine
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: /*
119: Create matrix. Here we only approximately preallocate storage space
120: for the Jacobian. See the users manual for a discussion of better
121: techniques for preallocating matrix memory.
122: */
123: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
124: MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,J,&shell);
125: MatShellSetOperation(shell,MATOP_MULT,(void(*)(void))mymatmult);
126: MatShellSetOperation(shell,MATOP_DESTROY,(void(*)(void))mymatdestroy);
127: ISCreateStride(PETSC_COMM_WORLD,user.mx*user.my,0,1,&is);
128: MatCreateNest(PETSC_COMM_WORLD,1,&is,1,&is,&shell,&nest);
129: MatDestroy(&shell);
130: ISDestroy(&is);
132: SNESSetJacobian(snes,nest,nest,FormJacobian,(void*)&user);
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Customize nonlinear solver; set runtime options
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: /*
139: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
140: */
141: SNESSetFromOptions(snes);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Evaluate initial guess; then solve nonlinear system
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: /*
148: Note: The user should initialize the vector, x, with the initial guess
149: for the nonlinear solver prior to calling SNESSolve(). In particular,
150: to employ an initial guess of zero, the user should explicitly set
151: this vector to zero by calling VecSet().
152: */
153: FormInitialGuess(&user,x);
154: SNESSolve(snes,NULL,x);
155: SNESGetIterationNumber(snes,&its);
156: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Free work space. All PETSc objects should be destroyed when they
160: are no longer needed.
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: MatDestroy(&nest);
164: VecDestroy(&x);
165: VecDestroy(&r);
166: SNESDestroy(&snes);
167: PetscFinalize();
168: return ierr;
169: }
171: /* ------------------------------------------------------------------- */
172: /*
173: FormInitialGuess - Forms initial approximation.
175: Input Parameters:
176: user - user-defined application context
177: X - vector
179: Output Parameter:
180: X - vector
181: */
182: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
183: {
184: PetscInt i,j,row,mx,my;
186: PetscReal lambda,temp1,temp,hx,hy;
187: PetscScalar *x;
189: mx = user->mx;
190: my = user->my;
191: lambda = user->param;
193: hx = 1.0 / (PetscReal)(mx-1);
194: hy = 1.0 / (PetscReal)(my-1);
196: /*
197: Get a pointer to vector data.
198: - For default PETSc vectors, VecGetArray() returns a pointer to
199: the data array. Otherwise, the routine is implementation dependent.
200: - You MUST call VecRestoreArray() when you no longer need access to
201: the array.
202: */
203: VecGetArray(X,&x);
204: temp1 = lambda/(lambda + 1.0);
205: for (j=0; j<my; j++) {
206: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
207: for (i=0; i<mx; i++) {
208: row = i + j*mx;
209: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
210: x[row] = 0.0;
211: continue;
212: }
213: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
214: }
215: }
217: /*
218: Restore vector
219: */
220: VecRestoreArray(X,&x);
221: return 0;
222: }
223: /* ------------------------------------------------------------------- */
224: /*
225: FormFunction - Evaluates nonlinear function, F(x).
227: Input Parameters:
228: . snes - the SNES context
229: . X - input vector
230: . ptr - optional user-defined context, as set by SNESSetFunction()
232: Output Parameter:
233: . F - function vector
234: */
235: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
236: {
237: AppCtx *user = (AppCtx*)ptr;
238: PetscInt i,j,row,mx,my;
239: PetscErrorCode ierr;
240: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
241: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
242: const PetscScalar *x;
244: mx = user->mx;
245: my = user->my;
246: lambda = user->param;
247: hx = one / (PetscReal)(mx-1);
248: hy = one / (PetscReal)(my-1);
249: sc = hx*hy;
250: hxdhy = hx/hy;
251: hydhx = hy/hx;
253: /*
254: Get pointers to vector data
255: */
256: VecGetArrayRead(X,&x);
257: VecGetArray(F,&f);
259: /*
260: Compute function
261: */
262: for (j=0; j<my; j++) {
263: for (i=0; i<mx; i++) {
264: row = i + j*mx;
265: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
266: f[row] = x[row];
267: continue;
268: }
269: u = x[row];
270: ub = x[row - mx];
271: ul = x[row - 1];
272: ut = x[row + mx];
273: ur = x[row + 1];
274: uxx = (-ur + two*u - ul)*hydhx;
275: uyy = (-ut + two*u - ub)*hxdhy;
276: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
277: }
278: }
280: /*
281: Restore vectors
282: */
283: VecRestoreArrayRead(X,&x);
284: VecRestoreArray(F,&f);
285: return 0;
286: }
287: /* ------------------------------------------------------------------- */
288: /*
289: FormJacobian - Evaluates Jacobian matrix.
291: Input Parameters:
292: . snes - the SNES context
293: . x - input vector
294: . ptr - optional user-defined context, as set by SNESSetJacobian()
296: Output Parameters:
297: . A - Jacobian matrix
298: . B - optionally different preconditioning matrix
299: . flag - flag indicating matrix structure
300: */
301: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat nest,Mat jac,void *ptr)
302: {
303: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
304: PetscInt i,j,row,mx,my,col[5];
305: PetscErrorCode ierr;
306: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
307: const PetscScalar *x;
308: PetscReal hx,hy,hxdhy,hydhx;
309: Mat J,shell;
311: MatNestGetSubMat(nest,0,0,&shell);
312: MatShellGetContext(shell,&J);
314: mx = user->mx;
315: my = user->my;
316: lambda = user->param;
317: hx = 1.0 / (PetscReal)(mx-1);
318: hy = 1.0 / (PetscReal)(my-1);
319: sc = hx*hy;
320: hxdhy = hx/hy;
321: hydhx = hy/hx;
323: /*
324: Get pointer to vector data
325: */
326: VecGetArrayRead(X,&x);
328: /*
329: Compute entries of the Jacobian
330: */
331: for (j=0; j<my; j++) {
332: for (i=0; i<mx; i++) {
333: row = i + j*mx;
334: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
335: MatSetValues(J,1,&row,1,&row,&one,INSERT_VALUES);
336: continue;
337: }
338: /* the 10 in the next line is intentionally an incorrect addition to the Jacobian */
339: v[0] = 10. + -hxdhy; col[0] = row - mx;
340: v[1] = -hydhx; col[1] = row - 1;
341: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
342: v[3] = -hydhx; col[3] = row + 1;
343: v[4] = -hxdhy; col[4] = row + mx;
344: MatSetValues(J,1,&row,5,col,v,INSERT_VALUES);
345: }
346: }
348: /*
349: Restore vector
350: */
351: VecRestoreArrayRead(X,&x);
353: /*
354: Assemble matrix
355: */
356: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
357: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
358: MatAssemblyBegin(shell,MAT_FINAL_ASSEMBLY);
359: MatAssemblyEnd(shell,MAT_FINAL_ASSEMBLY);
360: MatAssemblyBegin(nest,MAT_FINAL_ASSEMBLY);
361: MatAssemblyEnd(nest,MAT_FINAL_ASSEMBLY);
362: return 0;
363: }