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