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;

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

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

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

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

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

 50:   VecSet(x, 1.0);

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

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

 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;

 78:   half = 0.5;

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

 84:   /* Calculate residual */
 85:   for (i = 0; i < N; ++i) {
 86:     /*
 87:        Test for domain error.
 88:        Artificial 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:       PetscPrintf(PETSC_COMM_WORLD, "DOMAIN ERROR: x=%g\n", (double)PetscRealPart(x[i]));
 94:       SNESSetFunctionDomainError(snes);
 95:     }
 96:     f[i] = x[i] * x[i] - user->value;
 97:   }
 98:   VecRestoreArrayRead(X, &x);
 99:   VecRestoreArray(F, &f);
100:   return 0;
101: }

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;

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

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

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

133: /*TEST

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

138:    test:
139:       args:  -snes_monitor_solution -snes_linesearch_monitor

141: TEST*/