Actual source code: ex241.cxx
1: static char help[] = "Artificial test to check that snes->domainerror is being reset appropriately";
3: /* ------------------------------------------------------------------------
5: Artificial test to check that snes->domainerror is being reset appropriately
7: ------------------------------------------------------------------------- */
9: #define PETSC_SKIP_COMPLEX
10: #include <petscsnes.h>
12: typedef struct {
13: PetscReal value; /* parameter in nonlinear function */
14: } AppCtx;
16: PetscErrorCode UserFunction(SNES,Vec,Vec,void*);
17: PetscErrorCode UserJacobian(SNES,Vec,Mat,Mat,void*);
19: int main(int argc,char **argv)
20: {
21: SNES snes;
22: Vec x,r;
23: Mat J;
24: PetscInt its;
25: AppCtx user;
26: PetscMPIInt size;
28: PetscInitialize(&argc,&argv,(char*)0,help);
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
32: /* Allocate vectors / matrix */
33: VecCreate(PETSC_COMM_WORLD,&x);
34: VecSetSizes(x,PETSC_DECIDE,1);
35: VecSetFromOptions(x);
36: VecDuplicate(x,&r);
38: MatCreateSeqAIJ(PETSC_COMM_WORLD,1,1,1,NULL,&J);
40: /* Create / set-up SNES */
41: SNESCreate(PETSC_COMM_WORLD,&snes);
42: SNESSetFunction(snes,r,UserFunction,&user);
43: SNESSetJacobian(snes,J,J,UserJacobian,&user);
44: SNESSetFromOptions(snes);
46: /* Set initial guess (=1) and target value */
47: user.value = 1e-4;
49: VecSet(x,1.0);
51: /* Set initial guess / solve */
52: SNESSolve(snes,NULL,x);
53: SNESGetIterationNumber(snes,&its);
54: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
55: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
57: /* Done */
58: VecDestroy(&x);
59: VecDestroy(&r);
60: MatDestroy(&J);
61: SNESDestroy(&snes);
62: PetscFinalize();
63: return 0;
64: }
66: /*
67: UserFunction - for nonlinear function x^2 - value = 0
68: */
69: PetscErrorCode UserFunction(SNES snes,Vec X,Vec F,void *ptr)
70: {
71: AppCtx *user = (AppCtx*)ptr;
72: PetscInt N,i;
73: PetscScalar *f;
74: PetscReal half;
75: const PetscScalar *x;
77: half = 0.5;
79: VecGetSize(X,&N);
80: VecGetArrayRead(X,&x);
81: VecGetArray(F,&f);
83: /* Calculate residual */
84: for (i=0; i<N; ++i) {
85: /*
86: Test for domain error.
87: Artifical test is applied. With starting value 1.0, first iterate will be 0.5 + user->value/2.
88: Declare (0.5-value,0.5+value) to be infeasible.
89: In later iterations, snes->domainerror should be cleared, allowing iterations in the feasible region to be accepted.
90: */
91: if ((half-user->value) < PetscRealPart(x[i]) && PetscRealPart(x[i]) < (half+user->value)) {
92: PetscPrintf(PETSC_COMM_WORLD,"DOMAIN ERROR: x=%g\n",(double)PetscRealPart(x[i]));
93: SNESSetFunctionDomainError(snes);
94: }
95: f[i] = x[i]*x[i] - user->value;
96: }
97: VecRestoreArrayRead(X,&x);
98: VecRestoreArray(F,&f);
99: return 0;
100: }
102: /*
103: UserJacobian - for nonlinear function x^2 - value = 0
104: */
105: PetscErrorCode UserJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
106: {
107: PetscInt N,i,row,col;
108: const PetscScalar *x;
109: PetscScalar v;
111: VecGetSize(X,&N);
112: VecGetArrayRead(X,&x);
114: /* Calculate Jacobian */
115: for (i=0; i<N; ++i) {
116: row = i;
117: col = i;
118: v = 2*x[i];
119: MatSetValues(jac,1,&row,1,&col,&v,INSERT_VALUES);
120: }
121: VecRestoreArrayRead(X,&x);
122: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
123: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
125: if (jac != J) {
126: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
128: }
129: return 0;
130: }
132: /*TEST
134: build:
135: requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex
137: test:
138: args: -snes_monitor_solution -snes_linesearch_monitor
140: TEST*/