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: #if !defined(PETSC_SKIP_COMPLEX)
10: #define PETSC_SKIP_COMPLEX
11: #endif
12: #include <petscsnes.h>
14: typedef struct {
15: PetscReal value; /* parameter in nonlinear function */
16: } AppCtx;
18: PetscErrorCode UserFunction(SNES, Vec, Vec, void *);
19: PetscErrorCode UserJacobian(SNES, Vec, Mat, Mat, void *);
21: int main(int argc, char **argv)
22: {
23: SNES snes;
24: Vec x, r;
25: Mat J;
26: PetscInt its;
27: AppCtx user;
28: PetscMPIInt size;
30: PetscFunctionBeginUser;
31: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
32: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
33: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
35: /* Allocate vectors / matrix */
36: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
37: PetscCall(VecSetSizes(x, PETSC_DECIDE, 1));
38: PetscCall(VecSetFromOptions(x));
39: PetscCall(VecDuplicate(x, &r));
41: PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 1, 1, 1, NULL, &J));
43: /* Create / set-up SNES */
44: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
45: PetscCall(SNESSetFunction(snes, r, UserFunction, &user));
46: PetscCall(SNESSetJacobian(snes, J, J, UserJacobian, &user));
47: PetscCall(SNESSetFromOptions(snes));
49: /* Set initial guess (=1) and target value */
50: user.value = 1e-4;
52: PetscCall(VecSet(x, 1.0));
54: /* Set initial guess / solve */
55: PetscCall(SNESSolve(snes, NULL, x));
56: PetscCall(SNESGetIterationNumber(snes, &its));
57: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
58: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
60: /* Done */
61: PetscCall(VecDestroy(&x));
62: PetscCall(VecDestroy(&r));
63: PetscCall(MatDestroy(&J));
64: PetscCall(SNESDestroy(&snes));
65: PetscCall(PetscFinalize());
66: return 0;
67: }
69: /*
70: UserFunction - for nonlinear function x^2 - value = 0
71: */
72: PetscErrorCode UserFunction(SNES snes, Vec X, Vec F, void *ptr)
73: {
74: AppCtx *user = (AppCtx *)ptr;
75: PetscInt N, i;
76: PetscScalar *f;
77: PetscReal half;
78: const PetscScalar *x;
80: half = 0.5;
81: PetscFunctionBeginUser;
82: PetscCall(VecGetSize(X, &N));
83: PetscCall(VecGetArrayRead(X, &x));
84: PetscCall(VecGetArray(F, &f));
86: /* Calculate residual */
87: for (i = 0; i < N; ++i) {
88: /*
89: Test for domain error.
90: Artificial test is applied. With starting value 1.0, first iterate will be 0.5 + user->value/2.
91: Declare (0.5-value,0.5+value) to be infeasible.
92: In later iterations, snes->domainerror should be cleared, allowing iterations in the feasible region to be accepted.
93: */
94: if ((half - user->value) < PetscRealPart(x[i]) && PetscRealPart(x[i]) < (half + user->value)) {
95: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DOMAIN ERROR: x=%g\n", (double)PetscRealPart(x[i])));
96: PetscCall(SNESSetFunctionDomainError(snes));
97: }
98: f[i] = x[i] * x[i] - user->value;
99: }
100: PetscCall(VecRestoreArrayRead(X, &x));
101: PetscCall(VecRestoreArray(F, &f));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /*
106: UserJacobian - for nonlinear function x^2 - value = 0
107: */
108: PetscErrorCode UserJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
109: {
110: PetscInt N, i, row, col;
111: const PetscScalar *x;
112: PetscScalar v;
114: PetscFunctionBeginUser;
115: PetscCall(VecGetSize(X, &N));
116: PetscCall(VecGetArrayRead(X, &x));
118: /* Calculate Jacobian */
119: for (i = 0; i < N; ++i) {
120: row = i;
121: col = i;
122: v = 2 * x[i];
123: PetscCall(MatSetValues(jac, 1, &row, 1, &col, &v, INSERT_VALUES));
124: }
125: PetscCall(VecRestoreArrayRead(X, &x));
126: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
127: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
129: if (jac != J) {
130: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
132: }
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /*TEST
138: build:
139: requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex
141: test:
142: args: -snes_monitor_solution -snes_linesearch_monitor
144: TEST*/