Actual source code: ex241.cxx

  1: static char help[] = "Artificial test to check that snes->domainerror is being reset appropriately";

  3: /* ------------------------------------------------------------------------

  5:     Artificial test to check that snes->domainerror is being reset appropriately

  7:   ------------------------------------------------------------------------- */

  9: #define PETSC_SKIP_COMPLEX
 10: #include <petscsnes.h>

 12: typedef struct {
 13:   PetscReal value;              /* parameter in nonlinear function */
 14: } AppCtx;

 16: PetscErrorCode UserFunction(SNES,Vec,Vec,void*);
 17: PetscErrorCode UserJacobian(SNES,Vec,Mat,Mat,void*);

 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;

 28:   PetscInitialize(&argc,&argv,(char*)0,help);
 29:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 32:   /* Allocate vectors / matrix */
 33:   VecCreate(PETSC_COMM_WORLD,&x);
 34:   VecSetSizes(x,PETSC_DECIDE,1);
 35:   VecSetFromOptions(x);
 36:   VecDuplicate(x,&r);

 38:   MatCreateSeqAIJ(PETSC_COMM_WORLD,1,1,1,NULL,&J);

 40:   /* Create / set-up SNES */
 41:   SNESCreate(PETSC_COMM_WORLD,&snes);
 42:   SNESSetFunction(snes,r,UserFunction,&user);
 43:   SNESSetJacobian(snes,J,J,UserJacobian,&user);
 44:   SNESSetFromOptions(snes);

 46:   /* Set initial guess (=1) and target value */
 47:   user.value = 1e-4;

 49:   VecSet(x,1.0);

 51:   /* Set initial guess / solve */
 52:   SNESSolve(snes,NULL,x);
 53:   SNESGetIterationNumber(snes,&its);
 54:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
 55:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);

 57:   /* Done */
 58:   VecDestroy(&x);
 59:   VecDestroy(&r);
 60:   MatDestroy(&J);
 61:   SNESDestroy(&snes);
 62:   PetscFinalize();
 63:   return 0;
 64: }

 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;

 77:   half = 0.5;

 79:   VecGetSize(X,&N);
 80:   VecGetArrayRead(X,&x);
 81:   VecGetArray(F,&f);

 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:       PetscPrintf(PETSC_COMM_WORLD,"DOMAIN ERROR: x=%g\n",(double)PetscRealPart(x[i]));
 93:       SNESSetFunctionDomainError(snes);
 94:     }
 95:     f[i] = x[i]*x[i] - user->value;
 96:   }
 97:   VecRestoreArrayRead(X,&x);
 98:   VecRestoreArray(F,&f);
 99:   return 0;
100: }

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;

111:   VecGetSize(X,&N);
112:   VecGetArrayRead(X,&x);

114:   /* Calculate Jacobian */
115:   for (i=0; i<N; ++i) {
116:     row = i;
117:     col = i;
118:     v = 2*x[i];
119:     MatSetValues(jac,1,&row,1,&col,&v,INSERT_VALUES);
120:   }
121:   VecRestoreArrayRead(X,&x);
122:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
123:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

125:   if (jac != J) {
126:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
127:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
128:   }
129:   return 0;
130: }

132: /*TEST

134:    build:
135:       requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex

137:    test:
138:       args:  -snes_monitor_solution -snes_linesearch_monitor

140: TEST*/