Actual source code: ksponly.c

  1: #include <petsc/private/snesimpl.h>

  3: typedef struct {
  4:   PetscBool transpose_solve;
  5: } SNES_KSPONLY;

  7: static PetscErrorCode SNESSolve_KSPONLY(SNES snes)
  8: {
  9:   SNES_KSPONLY *ksponly = (SNES_KSPONLY *)snes->data;
 10:   PetscInt      lits;
 11:   Vec           Y, X, F;

 13:   PetscFunctionBegin;
 14:   PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

 16:   snes->numFailures            = 0;
 17:   snes->numLinearSolveFailures = 0;
 18:   snes->reason                 = SNES_CONVERGED_ITERATING;
 19:   snes->iter                   = 0;
 20:   snes->norm                   = 0.0;

 22:   X = snes->vec_sol;
 23:   F = snes->vec_func;
 24:   Y = snes->vec_sol_update;

 26:   if (!snes->vec_func_init_set) {
 27:     PetscCall(SNESComputeFunction(snes, X, F));
 28:   } else snes->vec_func_init_set = PETSC_FALSE;

 30:   if (snes->numbermonitors) {
 31:     PetscReal fnorm;
 32:     PetscCall(VecNorm(F, NORM_2, &fnorm));
 33:     PetscCall(SNESMonitor(snes, 0, fnorm));
 34:   }

 36:   /* Call general purpose update function */
 37:   PetscTryTypeMethod(snes, update, 0);

 39:   /* Solve J Y = F, where J is Jacobian matrix */
 40:   PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));

 42:   SNESCheckJacobianDomainerror(snes);

 44:   PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
 45:   if (ksponly->transpose_solve) {
 46:     PetscCall(KSPSolveTranspose(snes->ksp, F, Y));
 47:   } else {
 48:     PetscCall(KSPSolve(snes->ksp, F, Y));
 49:   }
 50:   snes->reason = SNES_CONVERGED_ITS;
 51:   SNESCheckKSPSolve(snes);

 53:   PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
 54:   PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
 55:   snes->iter++;

 57:   /* Take the computed step. */
 58:   PetscCall(VecAXPY(X, -1.0, Y));
 59:   if (snes->numbermonitors) {
 60:     PetscReal fnorm;
 61:     PetscCall(SNESComputeFunction(snes, X, F));
 62:     PetscCall(VecNorm(F, NORM_2, &fnorm));
 63:     PetscCall(SNESMonitor(snes, 1, fnorm));
 64:   }
 65:   PetscFunctionReturn(PETSC_SUCCESS);
 66: }

 68: static PetscErrorCode SNESSetUp_KSPONLY(SNES snes)
 69: {
 70:   PetscFunctionBegin;
 71:   PetscCall(SNESSetUpMatrices(snes));
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

 75: static PetscErrorCode SNESDestroy_KSPONLY(SNES snes)
 76: {
 77:   PetscFunctionBegin;
 78:   PetscCall(PetscFree(snes->data));
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: /*MC
 83:       SNESKSPONLY - Nonlinear solver that performs one Newton step and does not compute any norms.
 84:       The main purpose of this solver is to solve linear problems using the `SNES` interface, without
 85:       any additional overhead in the form of vector operations.

 87:    Level: beginner

 89: .seealso: `SNES`, `SNESType`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`
 90: M*/
 91: PETSC_EXTERN PetscErrorCode SNESCreate_KSPONLY(SNES snes)
 92: {
 93:   SNES_KSPONLY *ksponly;

 95:   PetscFunctionBegin;
 96:   snes->ops->setup          = SNESSetUp_KSPONLY;
 97:   snes->ops->solve          = SNESSolve_KSPONLY;
 98:   snes->ops->destroy        = SNESDestroy_KSPONLY;
 99:   snes->ops->setfromoptions = NULL;
100:   snes->ops->view           = NULL;
101:   snes->ops->reset          = NULL;

103:   snes->usesksp = PETSC_TRUE;
104:   snes->usesnpc = PETSC_FALSE;

106:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

108:   PetscCall(PetscNew(&ksponly));
109:   snes->data = (void *)ksponly;
110:   PetscFunctionReturn(PETSC_SUCCESS);
111: }

113: /*MC
114:       SNESKSPTRANSPOSEONLY - Nonlinear solver that performs one Newton step and does not compute any norms.
115:       The main purpose of this solver is to solve transposed linear problems using the `SNES` interface, without
116:       any additional overhead in the form of vector operations within adjoint solvers.

118:    Level: beginner

120: .seealso: `SNES`, `SNESType`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESKSPTRANSPOSEONLY`, `SNESNEWTONLS`, `SNESNEWTONTR`
121: M*/
122: PETSC_EXTERN PetscErrorCode SNESCreate_KSPTRANSPOSEONLY(SNES snes)
123: {
124:   SNES_KSPONLY *kspo;

126:   PetscFunctionBegin;
127:   PetscCall(SNESCreate_KSPONLY(snes));
128:   PetscCall(PetscObjectChangeTypeName((PetscObject)snes, SNESKSPTRANSPOSEONLY));
129:   kspo                  = (SNES_KSPONLY *)snes->data;
130:   kspo->transpose_solve = PETSC_TRUE;
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }