1 static char help[] = "Artificial test to check that snes->functiondomainerror is being reset appropriately";
2
3 /* ------------------------------------------------------------------------
4
5 Artificial test to check that snes->functiondomainerror is being reset appropriately
6
7 ------------------------------------------------------------------------- */
8
9 #if !defined(PETSC_SKIP_COMPLEX)
10 #define PETSC_SKIP_COMPLEX
11 #endif
12 #include <petscsnes.h>
13
14 typedef struct {
15 PetscReal value; /* parameter in nonlinear function */
16 } AppCtx;
17
18 PetscErrorCode UserFunction(SNES, Vec, Vec, void *);
19 PetscErrorCode UserJacobian(SNES, Vec, Mat, Mat, void *);
20
main(int argc,char ** argv)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;
29
30 PetscFunctionBeginUser;
31 PetscCall(PetscInitialize(&argc, &argv, nullptr, 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!");
34
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));
40
41 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 1, 1, 1, NULL, &J));
42
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));
48
49 /* Set initial guess (=1) and target value */
50 user.value = 1e-4;
51
52 PetscCall(VecSet(x, 1.0));
53
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));
59
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 }
68
69 /*
70 UserFunction - for nonlinear function x^2 - value = 0
71 */
UserFunction(SNES snes,Vec X,Vec F,void * ptr)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 = 0.5;
78 const PetscScalar *x;
79
80 PetscFunctionBeginUser;
81 PetscCall(VecGetSize(X, &N));
82 PetscCall(VecGetArrayRead(X, &x));
83 PetscCall(VecGetArray(F, &f));
84
85 /* Calculate residual */
86 for (i = 0; i < N; ++i) {
87 /*
88 Test for domain error.
89 Artificial test is applied. With starting value 1.0, first iterate will be 0.5 + user->value/2.
90 Declare (0.5-value,0.5+value) to be infeasible.
91 In later iterations, snes->functiondomainerror should be cleared, allowing iterations in the feasible region to be accepted.
92 */
93 if ((half - user->value) < PetscRealPart(x[i]) && PetscRealPart(x[i]) < (half + user->value)) {
94 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DOMAIN ERROR: x=%g\n", (double)PetscRealPart(x[i])));
95 PetscCall(SNESSetFunctionDomainError(snes));
96 }
97 f[i] = x[i] * x[i] - user->value;
98 }
99 PetscCall(VecRestoreArrayRead(X, &x));
100 PetscCall(VecRestoreArray(F, &f));
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
104 /*
105 UserJacobian - for nonlinear function x^2 - value = 0
106 */
UserJacobian(SNES snes,Vec X,Mat J,Mat jac,void * ptr)107 PetscErrorCode UserJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
108 {
109 PetscInt N, i, row, col;
110 const PetscScalar *x;
111 PetscScalar v;
112
113 PetscFunctionBeginUser;
114 PetscCall(VecGetSize(X, &N));
115 PetscCall(VecGetArrayRead(X, &x));
116
117 /* Calculate Jacobian */
118 for (i = 0; i < N; ++i) {
119 row = i;
120 col = i;
121 v = 2 * x[i];
122 PetscCall(MatSetValues(jac, 1, &row, 1, &col, &v, INSERT_VALUES));
123 }
124 PetscCall(VecRestoreArrayRead(X, &x));
125 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
126 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
127
128 if (jac != J) {
129 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
130 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
131 }
132 PetscFunctionReturn(PETSC_SUCCESS);
133 }
134
135 /*TEST
136
137 build:
138 requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex
139
140 test:
141 args: -snes_monitor_solution -snes_linesearch_monitor
142
143 TEST*/
144