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: #if !defined(PETSC_SKIP_COMPLEX)
 10:   #define PETSC_SKIP_COMPLEX
 11: #endif
 12: #include <petscsnes.h>

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

 18: PetscErrorCode UserFunction(SNES, Vec, Vec, void *);
 19: PetscErrorCode UserJacobian(SNES, Vec, Mat, Mat, void *);

 21: int main(int argc, char **argv)
 22: {
 23:   SNES        snes;
 24:   Vec         x, r;
 25:   Mat         J;
 26:   PetscInt    its;
 27:   AppCtx      user;
 28:   PetscMPIInt size;

 30:   PetscFunctionBeginUser;
 31:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 32:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 33:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

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

 41:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 1, 1, 1, NULL, &J));

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

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

 52:   PetscCall(VecSet(x, 1.0));

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

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

 69: /*
 70:     UserFunction - for nonlinear function x^2 - value = 0
 71: */
 72: PetscErrorCode UserFunction(SNES snes, Vec X, Vec F, void *ptr)
 73: {
 74:   AppCtx            *user = (AppCtx *)ptr;
 75:   PetscInt           N, i;
 76:   PetscScalar       *f;
 77:   PetscReal          half;
 78:   const PetscScalar *x;

 80:   half = 0.5;
 81:   PetscFunctionBeginUser;
 82:   PetscCall(VecGetSize(X, &N));
 83:   PetscCall(VecGetArrayRead(X, &x));
 84:   PetscCall(VecGetArray(F, &f));

 86:   /* Calculate residual */
 87:   for (i = 0; i < N; ++i) {
 88:     /*
 89:        Test for domain error.
 90:        Artificial test is applied.  With starting value 1.0, first iterate will be 0.5 + user->value/2.
 91:        Declare (0.5-value,0.5+value) to be infeasible.
 92:        In later iterations, snes->domainerror should be cleared, allowing iterations in the feasible region to be accepted.
 93:     */
 94:     if ((half - user->value) < PetscRealPart(x[i]) && PetscRealPart(x[i]) < (half + user->value)) {
 95:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DOMAIN ERROR: x=%g\n", (double)PetscRealPart(x[i])));
 96:       PetscCall(SNESSetFunctionDomainError(snes));
 97:     }
 98:     f[i] = x[i] * x[i] - user->value;
 99:   }
100:   PetscCall(VecRestoreArrayRead(X, &x));
101:   PetscCall(VecRestoreArray(F, &f));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: /*
106:     UserJacobian - for nonlinear function x^2 - value = 0
107: */
108: PetscErrorCode UserJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
109: {
110:   PetscInt           N, i, row, col;
111:   const PetscScalar *x;
112:   PetscScalar        v;

114:   PetscFunctionBeginUser;
115:   PetscCall(VecGetSize(X, &N));
116:   PetscCall(VecGetArrayRead(X, &x));

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

129:   if (jac != J) {
130:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
131:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
132:   }
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: /*TEST

138:    build:
139:       requires: !single !defined(PETSC_HAVE_SUN_CXX) !complex

141:    test:
142:       args:  -snes_monitor_solution -snes_linesearch_monitor

144: TEST*/