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;
     12:   SNESNormSchedule normschedule;

     14:   PetscFunctionBegin;
     15:   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);

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

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

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

     30:   PetscCall(SNESGetNormSchedule(snes, &normschedule));
     31:   if (snes->numbermonitors && (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY)) {
     32:     PetscReal fnorm;
     33:     PetscCall(VecNorm(F, NORM_2, &fnorm));
     34:     SNESCheckFunctionDomainError(snes, fnorm);
     35:     PetscCall(SNESMonitor(snes, 0, fnorm));
     36:   }

     38:   /* Call general purpose update function */
     39:   PetscTryTypeMethod(snes, update, 0);

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

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

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

     58:   /* Take the computed step. */
     59:   PetscCall(VecAXPY(X, -1.0, Y));

     61:   if (snes->numbermonitors && (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_FINAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY)) {
     62:     PetscReal fnorm;
     63:     PetscCall(SNESComputeFunction(snes, X, F));
     64:     PetscCall(VecNorm(F, NORM_2, &fnorm));
     65:     SNESCheckFunctionDomainError(snes, fnorm);
     66:     PetscCall(SNESMonitor(snes, 1, fnorm));
     67:   }
     68:   PetscFunctionReturn(PETSC_SUCCESS);
     69: }

     71: static PetscErrorCode SNESSetUp_KSPONLY(SNES snes)
     72: {
     73:   PetscFunctionBegin;
     74:   PetscCall(SNESSetUpMatrices(snes));
     75:   PetscFunctionReturn(PETSC_SUCCESS);
     76: }

     78: static PetscErrorCode SNESDestroy_KSPONLY(SNES snes)
     79: {
     80:   PetscFunctionBegin;
     81:   PetscCall(PetscFree(snes->data));
     82:   PetscFunctionReturn(PETSC_SUCCESS);
     83: }

     85: /*MC
     86:    SNESKSPONLY - Nonlinear solver that performs one Newton step with `KSPSolve()` and does not compute any norms.

     88:    Level: beginner

     90:    Notes:
     91:    The main purpose of this solver is to solve linear problems using the `SNES` interface, without
     92:    any additional overhead in the form of vector norm operations.

     94:    Use `SNESKSPTRANSPOSEONLY` for solving transposed systems.

     96: .seealso: [](ch_snes), `SNES`, `SNESType`, `SNESKSPTRANSPOSEONLY`, `SNESCreate()`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`
     97: M*/
     98: PETSC_EXTERN PetscErrorCode SNESCreate_KSPONLY(SNES snes)
     99: {
    100:   SNES_KSPONLY *ksponly;

    102:   PetscFunctionBegin;
    103:   snes->ops->setup          = SNESSetUp_KSPONLY;
    104:   snes->ops->solve          = SNESSolve_KSPONLY;
    105:   snes->ops->destroy        = SNESDestroy_KSPONLY;
    106:   snes->ops->setfromoptions = NULL;
    107:   snes->ops->view           = NULL;
    108:   snes->ops->reset          = NULL;

    110:   snes->usesksp = PETSC_TRUE;
    111:   snes->usesnpc = PETSC_FALSE;

    113:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

    115:   PetscCall(SNESParametersInitialize(snes));

    117:   PetscCall(PetscNew(&ksponly));
    118:   snes->data = (void *)ksponly;
    119:   PetscFunctionReturn(PETSC_SUCCESS);
    120: }

    122: /*MC
    123:    SNESKSPTRANSPOSEONLY - Nonlinear solver that performs one linear solve with `KSPSolveTranspose()` and does not compute any norms.

    125:    Level: beginner

    127:    Notes:
    128:    The main purpose of this solver is to solve transposed linear problems using the `SNES` interface, without
    129:    any additional overhead in the form of vector operations within adjoint solvers.

    131:    Use `SNESKSPONLY` for solving non-transposed systems.

    133: .seealso: [](ch_snes), `SNES`, `SNESType`, `SNESKSPONLY`, `SNESCreate()`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`
    134: M*/
    135: PETSC_EXTERN PetscErrorCode SNESCreate_KSPTRANSPOSEONLY(SNES snes)
    136: {
    137:   SNES_KSPONLY *kspo;

    139:   PetscFunctionBegin;
    140:   PetscCall(SNESCreate_KSPONLY(snes));
    141:   PetscCall(PetscObjectChangeTypeName((PetscObject)snes, SNESKSPTRANSPOSEONLY));
    142:   kspo                  = (SNES_KSPONLY *)snes->data;
    143:   kspo->transpose_solve = PETSC_TRUE;
    144:   PetscFunctionReturn(PETSC_SUCCESS);
    145: }