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;
29: PetscInitialize(&argc, &argv, (char *)0, help);
30: MPI_Comm_size(PETSC_COMM_WORLD, &size);
33: /* Allocate vectors / matrix */
34: VecCreate(PETSC_COMM_WORLD, &x);
35: VecSetSizes(x, PETSC_DECIDE, 1);
36: VecSetFromOptions(x);
37: VecDuplicate(x, &r);
39: MatCreateSeqAIJ(PETSC_COMM_WORLD, 1, 1, 1, NULL, &J);
41: /* Create / set-up SNES */
42: SNESCreate(PETSC_COMM_WORLD, &snes);
43: SNESSetFunction(snes, r, UserFunction, &user);
44: SNESSetJacobian(snes, J, J, UserJacobian, &user);
45: SNESSetFromOptions(snes);
47: /* Set initial guess (=1) and target value */
48: user.value = 1e-4;
50: VecSet(x, 1.0);
52: /* Set initial guess / solve */
53: SNESSolve(snes, NULL, x);
54: SNESGetIterationNumber(snes, &its);
55: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its);
56: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
58: /* Done */
59: VecDestroy(&x);
60: VecDestroy(&r);
61: MatDestroy(&J);
62: SNESDestroy(&snes);
63: PetscFinalize();
64: return 0;
65: }
67: /*
68: UserFunction - for nonlinear function x^2 - value = 0
69: */
70: PetscErrorCode UserFunction(SNES snes, Vec X, Vec F, void *ptr)
71: {
72: AppCtx *user = (AppCtx *)ptr;
73: PetscInt N, i;
74: PetscScalar *f;
75: PetscReal half;
76: const PetscScalar *x;
78: half = 0.5;
80: VecGetSize(X, &N);
81: VecGetArrayRead(X, &x);
82: VecGetArray(F, &f);
84: /* Calculate residual */
85: for (i = 0; i < N; ++i) {
86: /*
87: Test for domain error.
88: Artificial test is applied. With starting value 1.0, first iterate will be 0.5 + user->value/2.
89: Declare (0.5-value,0.5+value) to be infeasible.
90: In later iterations, snes->domainerror should be cleared, allowing iterations in the feasible region to be accepted.
91: */
92: if ((half - user->value) < PetscRealPart(x[i]) && PetscRealPart(x[i]) < (half + user->value)) {
93: PetscPrintf(PETSC_COMM_WORLD, "DOMAIN ERROR: x=%g\n", (double)PetscRealPart(x[i]));
94: SNESSetFunctionDomainError(snes);
95: }
96: f[i] = x[i] * x[i] - user->value;
97: }
98: VecRestoreArrayRead(X, &x);
99: VecRestoreArray(F, &f);
100: return 0;
101: }
103: /*
104: UserJacobian - for nonlinear function x^2 - value = 0
105: */
106: PetscErrorCode UserJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
107: {
108: PetscInt N, i, row, col;
109: const PetscScalar *x;
110: PetscScalar v;
112: VecGetSize(X, &N);
113: VecGetArrayRead(X, &x);
115: /* Calculate Jacobian */
116: for (i = 0; i < N; ++i) {
117: row = i;
118: col = i;
119: v = 2 * x[i];
120: MatSetValues(jac, 1, &row, 1, &col, &v, INSERT_VALUES);
121: }
122: VecRestoreArrayRead(X, &x);
123: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
124: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
126: if (jac != J) {
127: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
129: }
130: return 0;
131: }
133: /*TEST
135: build:
136: requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex
138: test:
139: args: -snes_monitor_solution -snes_linesearch_monitor
141: TEST*/