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