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