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 { 21 SNES snes; 22 Vec x,r; 23 Mat J; 24 PetscInt its; 25 AppCtx user; 26 PetscMPIInt size; 27 28 PetscFunctionBeginUser; 29 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 30 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 31 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 32 33 /* Allocate vectors / matrix */ 34 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 35 PetscCall(VecSetSizes(x,PETSC_DECIDE,1)); 36 PetscCall(VecSetFromOptions(x)); 37 PetscCall(VecDuplicate(x,&r)); 38 39 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD,1,1,1,NULL,&J)); 40 41 /* Create / set-up SNES */ 42 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 43 PetscCall(SNESSetFunction(snes,r,UserFunction,&user)); 44 PetscCall(SNESSetJacobian(snes,J,J,UserJacobian,&user)); 45 PetscCall(SNESSetFromOptions(snes)); 46 47 /* Set initial guess (=1) and target value */ 48 user.value = 1e-4; 49 50 PetscCall(VecSet(x,1.0)); 51 52 /* Set initial guess / solve */ 53 PetscCall(SNESSolve(snes,NULL,x)); 54 PetscCall(SNESGetIterationNumber(snes,&its)); 55 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %" PetscInt_FMT "\n",its)); 56 PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 57 58 /* Done */ 59 PetscCall(VecDestroy(&x)); 60 PetscCall(VecDestroy(&r)); 61 PetscCall(MatDestroy(&J)); 62 PetscCall(SNESDestroy(&snes)); 63 PetscCall(PetscFinalize()); 64 return 0; 65 } 66 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; 77 78 half = 0.5; 79 80 PetscCall(VecGetSize(X,&N)); 81 PetscCall(VecGetArrayRead(X,&x)); 82 PetscCall(VecGetArray(F,&f)); 83 84 /* Calculate residual */ 85 for (i=0; i<N; ++i) { 86 /* 87 Test for domain error. 88 Artifical 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 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"DOMAIN ERROR: x=%g\n",(double)PetscRealPart(x[i]))); 94 PetscCall(SNESSetFunctionDomainError(snes)); 95 } 96 f[i] = x[i]*x[i] - user->value; 97 } 98 PetscCall(VecRestoreArrayRead(X,&x)); 99 PetscCall(VecRestoreArray(F,&f)); 100 return 0; 101 } 102 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; 111 112 PetscCall(VecGetSize(X,&N)); 113 PetscCall(VecGetArrayRead(X,&x)); 114 115 /* Calculate Jacobian */ 116 for (i=0; i<N; ++i) { 117 row = i; 118 col = i; 119 v = 2*x[i]; 120 PetscCall(MatSetValues(jac,1,&row,1,&col,&v,INSERT_VALUES)); 121 } 122 PetscCall(VecRestoreArrayRead(X,&x)); 123 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 124 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 125 126 if (jac != J) { 127 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 128 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 129 } 130 return 0; 131 } 132 133 /*TEST 134 135 build: 136 requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex 137 138 test: 139 args: -snes_monitor_solution -snes_linesearch_monitor 140 141 TEST*/ 142