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 #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 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, (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!"); 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 */ 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; 79 80 half = 0.5; 81 PetscFunctionBeginUser; 82 PetscCall(VecGetSize(X, &N)); 83 PetscCall(VecGetArrayRead(X, &x)); 84 PetscCall(VecGetArray(F, &f)); 85 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 } 104 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; 113 114 PetscFunctionBeginUser; 115 PetscCall(VecGetSize(X, &N)); 116 PetscCall(VecGetArrayRead(X, &x)); 117 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)); 128 129 if (jac != J) { 130 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 131 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 132 } 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 136 /*TEST 137 138 build: 139 requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex 140 141 test: 142 args: -snes_monitor_solution -snes_linesearch_monitor 143 144 TEST*/ 145