#include /*I "petscsnes.h" I*/ #include #include #include #include #include PetscBool SNESRegisterAllCalled = PETSC_FALSE; PetscFunctionList SNESList = NULL; /* Logging support */ PetscClassId SNES_CLASSID, DMSNES_CLASSID; PetscLogEvent SNES_Solve, SNES_SetUp, SNES_FunctionEval, SNES_JacobianEval, SNES_NGSEval, SNES_NGSFuncEval, SNES_NewtonALEval, SNES_NPCSolve, SNES_ObjectiveEval; /*@ SNESSetErrorIfNotConverged - Causes `SNESSolve()` to generate an error immediately if the solver has not converged. Logically Collective Input Parameters: + snes - iterative context obtained from `SNESCreate()` - flg - `PETSC_TRUE` indicates you want the error generated Options Database Key: . -snes_error_if_not_converged - cause an immediate error condition and stop the program if the solver does not converge Level: intermediate Note: Normally PETSc continues if a solver fails to converge, you can call `SNESGetConvergedReason()` after a `SNESSolve()` to determine if it has converged. Otherwise the solution may be inaccurate or wrong .seealso: [](ch_snes), `SNES`, `SNESGetErrorIfNotConverged()`, `KSPGetErrorIfNotConverged()`, `KSPSetErrorIfNotConverged()` @*/ PetscErrorCode SNESSetErrorIfNotConverged(SNES snes, PetscBool flg) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveBool(snes, flg, 2); snes->errorifnotconverged = flg; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetErrorIfNotConverged - Indicates if `SNESSolve()` will generate an error if the solver does not converge? Not Collective Input Parameter: . snes - iterative context obtained from `SNESCreate()` Output Parameter: . flag - `PETSC_TRUE` if it will generate an error, else `PETSC_FALSE` Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESSetErrorIfNotConverged()`, `KSPGetErrorIfNotConverged()`, `KSPSetErrorIfNotConverged()` @*/ PetscErrorCode SNESGetErrorIfNotConverged(SNES snes, PetscBool *flag) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(flag, 2); *flag = snes->errorifnotconverged; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetAlwaysComputesFinalResidual - tells the `SNES` to always compute the residual (nonlinear function value) at the final solution Logically Collective Input Parameters: + snes - the shell `SNES` - flg - `PETSC_TRUE` to always compute the residual Level: advanced Note: Some solvers (such as smoothers in a `SNESFAS`) do not need the residual computed at the final solution so skip computing it to save time. .seealso: [](ch_snes), `SNES`, `SNESFAS`, `SNESSolve()`, `SNESGetAlwaysComputesFinalResidual()` @*/ PetscErrorCode SNESSetAlwaysComputesFinalResidual(SNES snes, PetscBool flg) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->alwayscomputesfinalresidual = flg; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetAlwaysComputesFinalResidual - checks if the `SNES` always computes the residual at the final solution Logically Collective Input Parameter: . snes - the `SNES` context Output Parameter: . flg - `PETSC_TRUE` if the residual is computed Level: advanced .seealso: [](ch_snes), `SNES`, `SNESFAS`, `SNESSolve()`, `SNESSetAlwaysComputesFinalResidual()` @*/ PetscErrorCode SNESGetAlwaysComputesFinalResidual(SNES snes, PetscBool *flg) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *flg = snes->alwayscomputesfinalresidual; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetFunctionDomainError - tells `SNES` that the input vector, a proposed new solution, to your function you provided to `SNESSetFunction()` is not in the functions domain. For example, a step with negative pressure. Not Collective Input Parameter: . snes - the `SNES` context Level: advanced Notes: This does not need to be called by all processes in the `SNES` MPI communicator. If this is called the `SNESSolve()` stops iterating and returns with a `SNESConvergedReason` of `SNES_DIVERGED_FUNCTION_DOMAIN` You should always call `SNESGetConvergedReason()` after each `SNESSolve()` and verify if the iteration converged (positive result) or diverged (negative result). You can direct `SNES` to avoid certain steps by using `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()` or `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()` You can call `SNESSetJacobianDomainError()` during a Jacobian computation to indicate the proposed solution is not in the domain. Developer Note: This value is used by `SNESCheckFunctionNorm()` to determine if the `SNESConvergedReason` is set to `SNES_DIVERGED_FUNCTION_DOMAIN` .seealso: [](ch_snes), `SNESCreate()`, `SNESSetFunction()`, `SNESFunctionFn`, `SNESSetJacobianDomainError()`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`, `SNESConvergedReason`, `SNESGetConvergedReason()`, `SNES_DIVERGED_FUNCTION_DOMAIN` @*/ PetscErrorCode SNESSetFunctionDomainError(SNES snes) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->domainerror = PETSC_TRUE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetJacobianDomainError - tells `SNES` that the function you provided to `SNESSetJacobian()` at the proposed step. For example there is a negative element transformation. Logically Collective Input Parameter: . snes - the `SNES` context Level: advanced Notes: If this is called the `SNESSolve()` stops iterating and returns with a `SNESConvergedReason` of `SNES_DIVERGED_FUNCTION_DOMAIN` You should always call `SNESGetConvergedReason()` after each `SNESSolve()` and verify if the iteration converged (positive result) or diverged (negative result). You can direct `SNES` to avoid certain steps by using `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()` or `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()` .seealso: [](ch_snes), `SNESCreate()`, `SNESSetFunction()`, `SNESFunctionFn`, `SNESSetFunctionDomainError()`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`, `SNESConvergedReason`, `SNESGetConvergedReason()` @*/ PetscErrorCode SNESSetJacobianDomainError(SNES snes) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCheck(!snes->errorifnotconverged, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "User code indicates computeJacobian does not make sense"); snes->jacobiandomainerror = PETSC_TRUE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetCheckJacobianDomainError - tells `SNESSolve()` whether to check if the user called `SNESSetJacobianDomainError()` Jacobian domain error after each Jacobian evaluation. By default, it checks for the Jacobian domain error in the debug mode, and does not check it in the optimized mode. Logically Collective Input Parameters: + snes - the `SNES` context - flg - indicates if or not to check Jacobian domain error after each Jacobian evaluation Level: advanced Note: Checks require one extra parallel synchronization for each Jacobian evaluation .seealso: [](ch_snes), `SNES`, `SNESConvergedReason`, `SNESCreate()`, `SNESSetFunction()`, `SNESFunctionFn`, `SNESSetFunctionDomainError()`, `SNESGetCheckJacobianDomainError()` @*/ PetscErrorCode SNESSetCheckJacobianDomainError(SNES snes, PetscBool flg) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->checkjacdomainerror = flg; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetCheckJacobianDomainError - Get an indicator whether or not `SNES` is checking Jacobian domain errors after each Jacobian evaluation. Logically Collective Input Parameter: . snes - the `SNES` context Output Parameter: . flg - `PETSC_FALSE` indicates that it is not checking Jacobian domain errors after each Jacobian evaluation Level: advanced .seealso: [](ch_snes), `SNES`, `SNESCreate()`, `SNESSetFunction()`, `SNESFunctionFn`, `SNESSetFunctionDomainError()`, `SNESSetCheckJacobianDomainError()` @*/ PetscErrorCode SNESGetCheckJacobianDomainError(SNES snes, PetscBool *flg) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(flg, 2); *flg = snes->checkjacdomainerror; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetFunctionDomainError - Gets the status of the domain error after a call to `SNESComputeFunction()` Not Collective, different MPI processes may return different values Input Parameter: . snes - the `SNES` context Output Parameter: . domainerror - Set to `PETSC_TRUE` if there's a domain error; `PETSC_FALSE` otherwise. Level: developer Notes: The value will only be true on those MPI processes that called `SNESSetFunctionDomainError()` The value is reset to `PETSC_FALSE` when `SNESCheckFunctionNorm()` is called. .seealso: [](ch_snes), `SNES`, `SNESSetFunctionDomainError()`, `SNESComputeFunction()` @*/ PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(domainerror, 2); *domainerror = snes->domainerror; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetJacobianDomainError - Gets the status of the Jacobian domain error after a call to `SNESComputeJacobian()` Not Collective, different MPI processes may return different values Input Parameter: . snes - the `SNES` context Output Parameter: . domainerror - Set to `PETSC_TRUE` if there's a Jacobian domain error; `PETSC_FALSE` otherwise. Level: advanced Notes: The value will only be true on those MPI processes that called `SNESSetJacobianDomainError()` The value is reset to `PETSC_FALSE` when `SNESCheckJacobianDomainerror()` is called but only `SNESSetCheckJacobianDomainError()` was called .seealso: [](ch_snes), `SNES`, `SNESSetFunctionDomainError()`, `SNESComputeFunction()`, `SNESGetFunctionDomainError()` @*/ PetscErrorCode SNESGetJacobianDomainError(SNES snes, PetscBool *domainerror) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(domainerror, 2); *domainerror = snes->jacobiandomainerror; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESLoad - Loads a `SNES` that has been stored in `PETSCVIEWERBINARY` with `SNESView()`. Collective Input Parameters: + snes - the newly loaded `SNES`, this needs to have been created with `SNESCreate()` or some related function before a call to `SNESLoad()`. - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` Level: intermediate Note: The `SNESType` is determined by the data in the file, any type set into the `SNES` before this call is ignored. .seealso: [](ch_snes), `SNES`, `PetscViewer`, `SNESCreate()`, `SNESType`, `PetscViewerBinaryOpen()`, `SNESView()`, `MatLoad()`, `VecLoad()` @*/ PetscErrorCode SNESLoad(SNES snes, PetscViewer viewer) { PetscBool isbinary; PetscInt classid; char type[256]; KSP ksp; DM dm; DMSNES dmsnes; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()"); PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT)); PetscCheck(classid == SNES_FILE_CLASSID, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Not SNES next in file"); PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR)); PetscCall(SNESSetType(snes, type)); PetscTryTypeMethod(snes, load, viewer); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &dmsnes)); PetscCall(DMSNESLoad(dmsnes, viewer)); PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(KSPLoad(ksp, viewer)); PetscFunctionReturn(PETSC_SUCCESS); } #include #if defined(PETSC_HAVE_SAWS) #include #endif /*@ SNESViewFromOptions - View a `SNES` based on values in the options database Collective Input Parameters: + A - the `SNES` context . obj - Optional object that provides the options prefix for the checks - name - command line option Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESView`, `PetscObjectViewFromOptions()`, `SNESCreate()` @*/ PetscErrorCode SNESViewFromOptions(SNES A, PetscObject obj, const char name[]) { PetscFunctionBegin; PetscValidHeaderSpecific(A, SNES_CLASSID, 1); PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name)); PetscFunctionReturn(PETSC_SUCCESS); } PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES, Vec, Mat, Mat, void *); /*@ SNESView - Prints or visualizes the `SNES` data structure. Collective Input Parameters: + snes - the `SNES` context - viewer - the `PetscViewer` Options Database Key: . -snes_view - Calls `SNESView()` at end of `SNESSolve()` Level: beginner Notes: The available visualization contexts include + `PETSC_VIEWER_STDOUT_SELF` - standard output (default) - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard output where only the first processor opens the file. All other processors send their data to the first processor to print. The available formats include + `PETSC_VIEWER_DEFAULT` - standard output (default) - `PETSC_VIEWER_ASCII_INFO_DETAIL` - more verbose output for `SNESNASM` The user can open an alternative visualization context with `PetscViewerASCIIOpen()` - output to a specified file. In the debugger you can do "call `SNESView`(snes,0)" to display the `SNES` solver. (The same holds for any PETSc object viewer). .seealso: [](ch_snes), `SNES`, `SNESLoad()`, `SNESCreate()`, `PetscViewerASCIIOpen()` @*/ PetscErrorCode SNESView(SNES snes, PetscViewer viewer) { SNESKSPEW *kctx; KSP ksp; SNESLineSearch linesearch; PetscBool isascii, isstring, isbinary, isdraw; DMSNES dmsnes; #if defined(PETSC_HAVE_SAWS) PetscBool issaws; #endif PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &viewer)); PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); PetscCheckSameComm(snes, 1, viewer, 2); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); #if defined(PETSC_HAVE_SAWS) PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws)); #endif if (isascii) { SNESNormSchedule normschedule; DM dm; SNESJacobianFn *cJ; void *ctx; const char *pre = ""; PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)snes, viewer)); if (!snes->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, " SNES has not been set up so information may be incomplete\n")); if (snes->ops->view) { PetscCall(PetscViewerASCIIPushTab(viewer)); PetscUseTypeMethod(snes, view, viewer); PetscCall(PetscViewerASCIIPopTab(viewer)); } if (snes->max_funcs == PETSC_UNLIMITED) { PetscCall(PetscViewerASCIIPrintf(viewer, " maximum iterations=%" PetscInt_FMT ", maximum function evaluations=unlimited\n", snes->max_its)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " maximum iterations=%" PetscInt_FMT ", maximum function evaluations=%" PetscInt_FMT "\n", snes->max_its, snes->max_funcs)); } PetscCall(PetscViewerASCIIPrintf(viewer, " tolerances: relative=%g, absolute=%g, solution=%g\n", (double)snes->rtol, (double)snes->abstol, (double)snes->stol)); if (snes->usesksp) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of linear solver iterations=%" PetscInt_FMT "\n", snes->linear_its)); PetscCall(PetscViewerASCIIPrintf(viewer, " total number of function evaluations=%" PetscInt_FMT "\n", snes->nfuncs)); PetscCall(SNESGetNormSchedule(snes, &normschedule)); if (normschedule > 0) PetscCall(PetscViewerASCIIPrintf(viewer, " norm schedule %s\n", SNESNormSchedules[normschedule])); if (snes->gridsequence) PetscCall(PetscViewerASCIIPrintf(viewer, " total number of grid sequence refinements=%" PetscInt_FMT "\n", snes->gridsequence)); if (snes->ksp_ewconv) { kctx = (SNESKSPEW *)snes->kspconvctx; if (kctx) { PetscCall(PetscViewerASCIIPrintf(viewer, " Eisenstat-Walker computation of KSP relative tolerance (version %" PetscInt_FMT ")\n", kctx->version)); PetscCall(PetscViewerASCIIPrintf(viewer, " rtol_0=%g, rtol_max=%g, threshold=%g\n", (double)kctx->rtol_0, (double)kctx->rtol_max, (double)kctx->threshold)); PetscCall(PetscViewerASCIIPrintf(viewer, " gamma=%g, alpha=%g, alpha2=%g\n", (double)kctx->gamma, (double)kctx->alpha, (double)kctx->alpha2)); } } if (snes->lagpreconditioner == -1) { PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioned is never rebuilt\n")); } else if (snes->lagpreconditioner > 1) { PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioned is rebuilt every %" PetscInt_FMT " new Jacobians\n", snes->lagpreconditioner)); } if (snes->lagjacobian == -1) { PetscCall(PetscViewerASCIIPrintf(viewer, " Jacobian is never rebuilt\n")); } else if (snes->lagjacobian > 1) { PetscCall(PetscViewerASCIIPrintf(viewer, " Jacobian is rebuilt every %" PetscInt_FMT " SNES iterations\n", snes->lagjacobian)); } PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESGetJacobian(dm, &cJ, &ctx)); if (snes->mf_operator) { PetscCall(PetscViewerASCIIPrintf(viewer, " Jacobian is applied matrix-free with differencing\n")); pre = "Preconditioning "; } if (cJ == SNESComputeJacobianDefault) { PetscCall(PetscViewerASCIIPrintf(viewer, " %sJacobian is built using finite differences one column at a time\n", pre)); } else if (cJ == SNESComputeJacobianDefaultColor) { PetscCall(PetscViewerASCIIPrintf(viewer, " %sJacobian is built using finite differences with coloring\n", pre)); /* it slightly breaks data encapsulation for access the DMDA information directly */ } else if (cJ == SNESComputeJacobian_DMDA) { MatFDColoring fdcoloring; PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring)); if (fdcoloring) { PetscCall(PetscViewerASCIIPrintf(viewer, " %sJacobian is built using colored finite differences on a DMDA\n", pre)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " %sJacobian is built using a DMDA local Jacobian\n", pre)); } } else if (snes->mf && !snes->mf_operator) { PetscCall(PetscViewerASCIIPrintf(viewer, " Jacobian is applied matrix-free with differencing, no explicit Jacobian\n")); } } else if (isstring) { const char *type; PetscCall(SNESGetType(snes, &type)); PetscCall(PetscViewerStringSPrintf(viewer, " SNESType: %-7.7s", type)); PetscTryTypeMethod(snes, view, viewer); } else if (isbinary) { PetscInt classid = SNES_FILE_CLASSID; MPI_Comm comm; PetscMPIInt rank; char type[256]; PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); PetscCallMPI(MPI_Comm_rank(comm, &rank)); if (rank == 0) { PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT)); PetscCall(PetscStrncpy(type, ((PetscObject)snes)->type_name, sizeof(type))); PetscCall(PetscViewerBinaryWrite(viewer, type, sizeof(type), PETSC_CHAR)); } PetscTryTypeMethod(snes, view, viewer); } else if (isdraw) { PetscDraw draw; char str[36]; PetscReal x, y, bottom, h; PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); PetscCall(PetscStrncpy(str, "SNES: ", sizeof(str))); PetscCall(PetscStrlcat(str, ((PetscObject)snes)->type_name, sizeof(str))); PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_BLUE, PETSC_DRAW_BLACK, str, NULL, &h)); bottom = y - h; PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom)); PetscTryTypeMethod(snes, view, viewer); #if defined(PETSC_HAVE_SAWS) } else if (issaws) { PetscMPIInt rank; const char *name; PetscCall(PetscObjectGetName((PetscObject)snes, &name)); PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); if (!((PetscObject)snes)->amsmem && rank == 0) { char dir[1024]; PetscCall(PetscObjectViewSAWs((PetscObject)snes, viewer)); PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/its", name)); PetscCallSAWs(SAWs_Register, (dir, &snes->iter, 1, SAWs_READ, SAWs_INT)); if (!snes->conv_hist) PetscCall(SNESSetConvergenceHistory(snes, NULL, NULL, PETSC_DECIDE, PETSC_TRUE)); PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/conv_hist", name)); PetscCallSAWs(SAWs_Register, (dir, snes->conv_hist, 10, SAWs_READ, SAWs_DOUBLE)); } #endif } if (snes->linesearch) { PetscCall(SNESGetLineSearch(snes, &linesearch)); PetscCall(PetscViewerASCIIPushTab(viewer)); PetscCall(SNESLineSearchView(linesearch, viewer)); PetscCall(PetscViewerASCIIPopTab(viewer)); } if (snes->npc && snes->usesnpc) { PetscCall(PetscViewerASCIIPushTab(viewer)); PetscCall(SNESView(snes->npc, viewer)); PetscCall(PetscViewerASCIIPopTab(viewer)); } PetscCall(PetscViewerASCIIPushTab(viewer)); PetscCall(DMGetDMSNES(snes->dm, &dmsnes)); PetscCall(DMSNESView(dmsnes, viewer)); PetscCall(PetscViewerASCIIPopTab(viewer)); if (snes->usesksp) { PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(PetscViewerASCIIPushTab(viewer)); PetscCall(KSPView(ksp, viewer)); PetscCall(PetscViewerASCIIPopTab(viewer)); } if (isdraw) { PetscDraw draw; PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); PetscCall(PetscDrawPopCurrentPoint(draw)); } PetscFunctionReturn(PETSC_SUCCESS); } /* We retain a list of functions that also take SNES command line options. These are called at the end SNESSetFromOptions() */ #define MAXSETFROMOPTIONS 5 static PetscInt numberofsetfromoptions; static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES); /*@C SNESAddOptionsChecker - Adds an additional function to check for `SNES` options. Not Collective Input Parameter: . snescheck - function that checks for options Calling sequence of `snescheck`: . snes - the `SNES` object for which it is checking options Level: developer .seealso: [](ch_snes), `SNES`, `SNESSetFromOptions()` @*/ PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES snes)) { PetscFunctionBegin; PetscCheck(numberofsetfromoptions < MAXSETFROMOPTIONS, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %d allowed", MAXSETFROMOPTIONS); othersetfromoptions[numberofsetfromoptions++] = snescheck; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version) { Mat J; MatNullSpace nullsp; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) { Mat A = snes->jacobian, B = snes->jacobian_pre; PetscCall(MatCreateVecs(A ? A : B, NULL, &snes->vec_func)); } PetscCheck(version == 1 || version == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator routines, only version 1 and 2"); if (version == 1) { PetscCall(MatCreateSNESMF(snes, &J)); PetscCall(MatMFFDSetOptionsPrefix(J, ((PetscObject)snes)->prefix)); PetscCall(MatSetFromOptions(J)); /* TODO: the version 2 code should be merged into the MatCreateSNESMF() and MatCreateMFFD() infrastructure and then removed */ } else /* if (version == 2) */ { PetscCheck(snes->vec_func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "SNESSetFunction() must be called first"); #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) && !defined(PETSC_USE_REAL___FP16) PetscCall(MatCreateSNESMFMore(snes, snes->vec_func, &J)); #else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "matrix-free operator routines (version 2)"); #endif } /* attach any user provided null space that was on Amat to the newly created matrix-free matrix */ if (snes->jacobian) { PetscCall(MatGetNullSpace(snes->jacobian, &nullsp)); if (nullsp) PetscCall(MatSetNullSpace(J, nullsp)); } PetscCall(PetscInfo(snes, "Setting default matrix-free operator routines (version %" PetscInt_FMT ")\n", version)); if (hasOperator) { /* This version replaces the user provided Jacobian matrix with a matrix-free version but still employs the user-provided matrix used for computing the preconditioner. */ PetscCall(SNESSetJacobian(snes, J, NULL, NULL, NULL)); } else { /* This version replaces both the user-provided Jacobian and the user- provided preconditioner Jacobian with the default matrix-free version. */ if (snes->npcside == PC_LEFT && snes->npc) { if (!snes->jacobian) PetscCall(SNESSetJacobian(snes, J, NULL, NULL, NULL)); } else { KSP ksp; PC pc; PetscBool match; PetscCall(SNESSetJacobian(snes, J, J, MatMFFDComputeJacobian, NULL)); /* Force no preconditioner */ PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(KSPGetPC(ksp, &pc)); PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &match, PCSHELL, PCH2OPUS, "")); if (!match) { PetscCall(PetscInfo(snes, "Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n")); PetscCall(PCSetType(pc, PCNONE)); } } } PetscCall(MatDestroy(&J)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMRestrictHook_SNESVecSol(DM dmfine, Mat Restrict, Vec Rscale, Mat Inject, DM dmcoarse, void *ctx) { SNES snes = (SNES)ctx; Vec Xfine, Xfine_named = NULL, Xcoarse; PetscFunctionBegin; if (PetscLogPrintInfo) { PetscInt finelevel, coarselevel, fineclevel, coarseclevel; PetscCall(DMGetRefineLevel(dmfine, &finelevel)); PetscCall(DMGetCoarsenLevel(dmfine, &fineclevel)); PetscCall(DMGetRefineLevel(dmcoarse, &coarselevel)); PetscCall(DMGetCoarsenLevel(dmcoarse, &coarseclevel)); PetscCall(PetscInfo(dmfine, "Restricting SNES solution vector from level %" PetscInt_FMT "-%" PetscInt_FMT " to level %" PetscInt_FMT "-%" PetscInt_FMT "\n", finelevel, fineclevel, coarselevel, coarseclevel)); } if (dmfine == snes->dm) Xfine = snes->vec_sol; else { PetscCall(DMGetNamedGlobalVector(dmfine, "SNESVecSol", &Xfine_named)); Xfine = Xfine_named; } PetscCall(DMGetNamedGlobalVector(dmcoarse, "SNESVecSol", &Xcoarse)); if (Inject) { PetscCall(MatRestrict(Inject, Xfine, Xcoarse)); } else { PetscCall(MatRestrict(Restrict, Xfine, Xcoarse)); PetscCall(VecPointwiseMult(Xcoarse, Xcoarse, Rscale)); } PetscCall(DMRestoreNamedGlobalVector(dmcoarse, "SNESVecSol", &Xcoarse)); if (Xfine_named) PetscCall(DMRestoreNamedGlobalVector(dmfine, "SNESVecSol", &Xfine_named)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMCoarsenHook_SNESVecSol(DM dm, DM dmc, void *ctx) { PetscFunctionBegin; PetscCall(DMCoarsenHookAdd(dmc, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } /* This may be called to rediscretize the operator on levels of linear multigrid. The DM shuffle is so the user can * safely call SNESGetDM() in their residual evaluation routine. */ static PetscErrorCode KSPComputeOperators_SNES(KSP ksp, Mat A, Mat B, void *ctx) { SNES snes = (SNES)ctx; DMSNES sdm; Vec X, Xnamed = NULL; DM dmsave; void *ctxsave; SNESJacobianFn *jac = NULL; PetscFunctionBegin; dmsave = snes->dm; PetscCall(KSPGetDM(ksp, &snes->dm)); if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */ else { PetscBool has; /* We are on a coarser level, this vec was initialized using a DM restrict hook */ PetscCall(DMHasNamedGlobalVector(snes->dm, "SNESVecSol", &has)); PetscCheck(has, PetscObjectComm((PetscObject)snes->dm), PETSC_ERR_PLIB, "Missing SNESVecSol"); PetscCall(DMGetNamedGlobalVector(snes->dm, "SNESVecSol", &Xnamed)); X = Xnamed; PetscCall(SNESGetJacobian(snes, NULL, NULL, &jac, &ctxsave)); /* If the DM's don't match up, the MatFDColoring context needed for the jacobian won't match up either -- fixit. */ if (jac == SNESComputeJacobianDefaultColor) PetscCall(SNESSetJacobian(snes, NULL, NULL, SNESComputeJacobianDefaultColor, NULL)); } /* Compute the operators */ PetscCall(DMGetDMSNES(snes->dm, &sdm)); if (Xnamed && sdm->ops->computefunction) { /* The SNES contract with the user is that ComputeFunction is always called before ComputeJacobian. We make sure of this here. Disable affine shift since it is for the finest level */ Vec F, saverhs = snes->vec_rhs; snes->vec_rhs = NULL; PetscCall(DMGetGlobalVector(snes->dm, &F)); PetscCall(SNESComputeFunction(snes, X, F)); PetscCall(DMRestoreGlobalVector(snes->dm, &F)); snes->vec_rhs = saverhs; snes->nfuncs--; /* Do not log coarser level evaluations */ } /* Make sure KSP DM has the Jacobian computation routine */ if (!sdm->ops->computejacobian) PetscCall(DMCopyDMSNES(dmsave, snes->dm)); PetscCall(SNESComputeJacobian(snes, X, A, B)); /* Put the previous context back */ if (snes->dm != dmsave && jac == SNESComputeJacobianDefaultColor) PetscCall(SNESSetJacobian(snes, NULL, NULL, jac, ctxsave)); if (Xnamed) PetscCall(DMRestoreNamedGlobalVector(snes->dm, "SNESVecSol", &Xnamed)); snes->dm = dmsave; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetUpMatrices - ensures that matrices are available for `SNES` Newton-like methods, this is called by `SNESSetUp_XXX()` Collective Input Parameter: . snes - `SNES` object to configure Level: developer Note: If the matrices do not yet exist it attempts to create them based on options previously set for the `SNES` such as `-snes_mf` Developer Note: The functionality of this routine overlaps in a confusing way with the functionality of `SNESSetUpMatrixFree_Private()` which is called by `SNESSetUp()` but sometimes `SNESSetUpMatrices()` is called without `SNESSetUp()` being called. A refactorization to simplify the logic that handles the matrix-free case is desirable. .seealso: [](ch_snes), `SNES`, `SNESSetUp()` @*/ PetscErrorCode SNESSetUpMatrices(SNES snes) { DM dm; DMSNES sdm; PetscFunctionBegin; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); if (!snes->jacobian && snes->mf && !snes->mf_operator && !snes->jacobian_pre) { Mat J; void *functx; PetscCall(MatCreateSNESMF(snes, &J)); PetscCall(MatMFFDSetOptionsPrefix(J, ((PetscObject)snes)->prefix)); PetscCall(MatSetFromOptions(J)); PetscCall(SNESGetFunction(snes, NULL, NULL, &functx)); PetscCall(SNESSetJacobian(snes, J, J, NULL, NULL)); PetscCall(MatDestroy(&J)); } else if (snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) { Mat J, B; PetscCall(MatCreateSNESMF(snes, &J)); PetscCall(MatMFFDSetOptionsPrefix(J, ((PetscObject)snes)->prefix)); PetscCall(MatSetFromOptions(J)); PetscCall(DMCreateMatrix(snes->dm, &B)); /* sdm->computejacobian was already set to reach here */ PetscCall(SNESSetJacobian(snes, J, B, NULL, NULL)); PetscCall(MatDestroy(&J)); PetscCall(MatDestroy(&B)); } else if (!snes->jacobian_pre) { PetscDS prob; Mat J, B; PetscBool hasPrec = PETSC_FALSE; J = snes->jacobian; PetscCall(DMGetDS(dm, &prob)); if (prob) PetscCall(PetscDSHasJacobianPreconditioner(prob, &hasPrec)); if (J) PetscCall(PetscObjectReference((PetscObject)J)); else if (hasPrec) PetscCall(DMCreateMatrix(snes->dm, &J)); PetscCall(DMCreateMatrix(snes->dm, &B)); PetscCall(SNESSetJacobian(snes, J ? J : B, B, NULL, NULL)); PetscCall(MatDestroy(&J)); PetscCall(MatDestroy(&B)); } { KSP ksp; PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(KSPSetComputeOperators(ksp, KSPComputeOperators_SNES, snes)); PetscCall(DMCoarsenHookAdd(snes->dm, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, snes)); } PetscFunctionReturn(PETSC_SUCCESS); } PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode PetscMonitorPauseFinal_Internal(PetscInt, void *); static PetscErrorCode SNESMonitorPauseFinal_Internal(SNES snes) { PetscFunctionBegin; if (!snes->pauseFinal) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(PetscMonitorPauseFinal_Internal(snes->numbermonitors, snes->monitorcontext)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user Collective Input Parameters: + snes - `SNES` object you wish to monitor . name - the monitor type one is seeking . help - message indicating what monitoring is done . manual - manual page for the monitor . monitor - the monitor function, this must use a `PetscViewerFormat` as its context - monitorsetup - a function that is called once ONLY if the user selected this monitor that may set additional features of the `SNES` or `PetscViewer` objects Calling sequence of `monitor`: + snes - the nonlinear solver context . it - the current iteration . r - the current function norm - vf - a `PetscViewerAndFormat` struct that contains the `PetscViewer` and `PetscViewerFormat` to use Calling sequence of `monitorsetup`: + snes - the nonlinear solver context - vf - a `PetscViewerAndFormat` struct that contains the `PetscViewer` and `PetscViewerFormat` to use Options Database Key: . -name - trigger the use of this monitor in `SNESSetFromOptions()` Level: advanced .seealso: [](ch_snes), `PetscOptionsCreateViewer()`, `PetscOptionsGetReal()`, `PetscOptionsHasName()`, `PetscOptionsGetString()`, `PetscOptionsGetIntArray()`, `PetscOptionsGetRealArray()`, `PetscOptionsBool()` `PetscOptionsInt()`, `PetscOptionsString()`, `PetscOptionsReal()`, `PetscOptionsName()`, `PetscOptionsBegin()`, `PetscOptionsEnd()`, `PetscOptionsHeadBegin()`, `PetscOptionsStringArray()`, `PetscOptionsRealArray()`, `PetscOptionsScalar()`, `PetscOptionsBoolGroupBegin()`, `PetscOptionsBoolGroup()`, `PetscOptionsBoolGroupEnd()`, `PetscOptionsFList()`, `PetscOptionsEList()` @*/ PetscErrorCode SNESMonitorSetFromOptions(SNES snes, const char name[], const char help[], const char manual[], PetscErrorCode (*monitor)(SNES snes, PetscInt it, PetscReal r, PetscViewerAndFormat *vf), PetscErrorCode (*monitorsetup)(SNES snes, PetscViewerAndFormat *vf)) { PetscViewer viewer; PetscViewerFormat format; PetscBool flg; PetscFunctionBegin; PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, name, &viewer, &format, &flg)); if (flg) { PetscViewerAndFormat *vf; PetscCall(PetscViewerAndFormatCreate(viewer, format, &vf)); PetscCall(PetscViewerDestroy(&viewer)); if (monitorsetup) PetscCall((*monitorsetup)(snes, vf)); PetscCall(SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *))monitor, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy)); } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode SNESEWSetFromOptions_Private(SNESKSPEW *kctx, PetscBool print_api, MPI_Comm comm, const char *prefix) { const char *api = print_api ? "SNESKSPSetParametersEW" : NULL; PetscFunctionBegin; PetscOptionsBegin(comm, prefix, "Eisenstat and Walker type forcing options", "KSP"); PetscCall(PetscOptionsInt("-ksp_ew_version", "Version 1, 2 or 3", api, kctx->version, &kctx->version, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_rtol0", "0 <= rtol0 < 1", api, kctx->rtol_0, &kctx->rtol_0, NULL)); kctx->rtol_max = PetscMax(kctx->rtol_0, kctx->rtol_max); PetscCall(PetscOptionsReal("-ksp_ew_rtolmax", "0 <= rtolmax < 1", api, kctx->rtol_max, &kctx->rtol_max, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_gamma", "0 <= gamma <= 1", api, kctx->gamma, &kctx->gamma, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_alpha", "1 < alpha <= 2", api, kctx->alpha, &kctx->alpha, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_alpha2", "alpha2", NULL, kctx->alpha2, &kctx->alpha2, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_threshold", "0 < threshold < 1", api, kctx->threshold, &kctx->threshold, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_v4_p1", "p1", NULL, kctx->v4_p1, &kctx->v4_p1, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_v4_p2", "p2", NULL, kctx->v4_p2, &kctx->v4_p2, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_v4_p3", "p3", NULL, kctx->v4_p3, &kctx->v4_p3, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_v4_m1", "Scaling when rk-1 in [p2,p3)", NULL, kctx->v4_m1, &kctx->v4_m1, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_v4_m2", "Scaling when rk-1 in [p3,+infty)", NULL, kctx->v4_m2, &kctx->v4_m2, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_v4_m3", "Threshold for successive rtol (0.1 in Eq.7)", NULL, kctx->v4_m3, &kctx->v4_m3, NULL)); PetscCall(PetscOptionsReal("-ksp_ew_v4_m4", "Adaptation scaling (0.5 in Eq.7)", NULL, kctx->v4_m4, &kctx->v4_m4, NULL)); PetscOptionsEnd(); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetFromOptions - Sets various `SNES` and `KSP` parameters from user options. Collective Input Parameter: . snes - the `SNES` context Options Database Keys: + -snes_type - newtonls, newtontr, ngmres, ncg, nrichardson, qn, vi, fas, `SNESType` for complete list . -snes_rtol - relative decrease in tolerance norm from initial . -snes_atol - absolute tolerance of residual norm . -snes_stol - convergence tolerance in terms of the norm of the change in the solution between steps . -snes_divergence_tolerance - if the residual goes above divtol*rnorm0, exit with divergence . -snes_max_it - maximum number of iterations . -snes_max_funcs - maximum number of function evaluations . -snes_force_iteration - force `SNESSolve()` to take at least one iteration . -snes_max_fail - maximum number of line search failures allowed before stopping, default is none . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops . -snes_lag_preconditioner - how often preconditioner is rebuilt (use -1 to never rebuild) . -snes_lag_preconditioner_persists - retains the -snes_lag_preconditioner information across multiple SNESSolve() . -snes_lag_jacobian - how often Jacobian is rebuilt (use -1 to never rebuild) . -snes_lag_jacobian_persists - retains the -snes_lag_jacobian information across multiple SNESSolve() . -snes_convergence_test - convergence test in nonlinear solver. default `SNESConvergedDefault()`. skip `SNESConvergedSkip()` means continue iterating until max_it or some other criterion is reached, saving expense of convergence test. correct_pressure `SNESConvergedCorrectPressure()` has special handling of a pressure null space. . -snes_monitor [ascii][:filename][:viewer format] - prints residual norm at each iteration. if no filename given prints to stdout . -snes_monitor_solution [ascii binary draw][:filename][:viewer format] - plots solution at each iteration . -snes_monitor_residual [ascii binary draw][:filename][:viewer format] - plots residual (not its norm) at each iteration . -snes_monitor_solution_update [ascii binary draw][:filename][:viewer format] - plots update to solution at each iteration . -snes_monitor_lg_residualnorm - plots residual norm at each iteration . -snes_monitor_lg_range - plots residual norm at each iteration . -snes_monitor_pause_final - Pauses all monitor drawing after the solver ends . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing . -snes_fd_color - use finite differences with coloring to compute Jacobian . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each `KSP` iteration . -snes_converged_reason - print the reason for convergence/divergence after each solve . -npc_snes_type - the `SNES` type to use as a nonlinear preconditioner . -snes_test_jacobian - compare the user provided Jacobian with one computed via finite differences to check for errors. If a threshold is given, display only those entries whose difference is greater than the threshold. - -snes_test_jacobian_view - display the user provided Jacobian, the finite difference Jacobian and the difference between them to help users detect the location of errors in the user provided Jacobian. Options Database Keys for Eisenstat-Walker method: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence . -snes_ksp_ew_version ver - version of Eisenstat-Walker method . -snes_ksp_ew_rtol0 - Sets rtol0 . -snes_ksp_ew_rtolmax - Sets rtolmax . -snes_ksp_ew_gamma - Sets gamma . -snes_ksp_ew_alpha - Sets alpha . -snes_ksp_ew_alpha2 - Sets alpha2 - -snes_ksp_ew_threshold - Sets threshold Level: beginner Notes: To see all options, run your program with the -help option or consult the users manual `SNES` supports three approaches for computing (approximate) Jacobians: user provided via `SNESSetJacobian()`, matrix-free using `MatCreateSNESMF()`, and computing explicitly with finite differences and coloring using `MatFDColoring`. It is also possible to use automatic differentiation and the `MatFDColoring` object. .seealso: [](ch_snes), `SNESType`, `SNESSetOptionsPrefix()`, `SNESResetFromOptions()`, `SNES`, `SNESCreate()`, `MatCreateSNESMF()`, `MatFDColoring` @*/ PetscErrorCode SNESSetFromOptions(SNES snes) { PetscBool flg, pcset, persist, set; PetscInt i, indx, lag, grids, max_its, max_funcs; const char *deft = SNESNEWTONLS; const char *convtests[] = {"default", "skip", "correct_pressure"}; SNESKSPEW *kctx = NULL; char type[256], monfilename[PETSC_MAX_PATH_LEN], ewprefix[256]; PCSide pcside; const char *optionsprefix; PetscReal rtol, abstol, stol; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(SNESRegisterAll()); PetscObjectOptionsBegin((PetscObject)snes); if (((PetscObject)snes)->type_name) deft = ((PetscObject)snes)->type_name; PetscCall(PetscOptionsFList("-snes_type", "Nonlinear solver method", "SNESSetType", SNESList, deft, type, 256, &flg)); if (flg) { PetscCall(SNESSetType(snes, type)); } else if (!((PetscObject)snes)->type_name) { PetscCall(SNESSetType(snes, deft)); } abstol = snes->abstol; rtol = snes->rtol; stol = snes->stol; max_its = snes->max_its; max_funcs = snes->max_funcs; PetscCall(PetscOptionsReal("-snes_rtol", "Stop if decrease in function norm less than", "SNESSetTolerances", snes->rtol, &rtol, NULL)); PetscCall(PetscOptionsReal("-snes_atol", "Stop if function norm less than", "SNESSetTolerances", snes->abstol, &abstol, NULL)); PetscCall(PetscOptionsReal("-snes_stol", "Stop if step length less than", "SNESSetTolerances", snes->stol, &stol, NULL)); PetscCall(PetscOptionsInt("-snes_max_it", "Maximum iterations", "SNESSetTolerances", snes->max_its, &max_its, NULL)); PetscCall(PetscOptionsInt("-snes_max_funcs", "Maximum function evaluations", "SNESSetTolerances", snes->max_funcs, &max_funcs, NULL)); PetscCall(SNESSetTolerances(snes, abstol, rtol, stol, max_its, max_funcs)); PetscCall(PetscOptionsReal("-snes_divergence_tolerance", "Stop if residual norm increases by this factor", "SNESSetDivergenceTolerance", snes->divtol, &snes->divtol, &flg)); if (flg) PetscCall(SNESSetDivergenceTolerance(snes, snes->divtol)); PetscCall(PetscOptionsInt("-snes_max_fail", "Maximum nonlinear step failures", "SNESSetMaxNonlinearStepFailures", snes->maxFailures, &snes->maxFailures, &flg)); if (flg) PetscCall(SNESSetMaxNonlinearStepFailures(snes, snes->maxFailures)); PetscCall(PetscOptionsInt("-snes_max_linear_solve_fail", "Maximum failures in linear solves allowed", "SNESSetMaxLinearSolveFailures", snes->maxLinearSolveFailures, &snes->maxLinearSolveFailures, &flg)); if (flg) PetscCall(SNESSetMaxLinearSolveFailures(snes, snes->maxLinearSolveFailures)); PetscCall(PetscOptionsBool("-snes_error_if_not_converged", "Generate error if solver does not converge", "SNESSetErrorIfNotConverged", snes->errorifnotconverged, &snes->errorifnotconverged, NULL)); PetscCall(PetscOptionsBool("-snes_force_iteration", "Force SNESSolve() to take at least one iteration", "SNESSetForceIteration", snes->forceiteration, &snes->forceiteration, NULL)); PetscCall(PetscOptionsBool("-snes_check_jacobian_domain_error", "Check Jacobian domain error after Jacobian evaluation", "SNESCheckJacobianDomainError", snes->checkjacdomainerror, &snes->checkjacdomainerror, NULL)); PetscCall(PetscOptionsInt("-snes_lag_preconditioner", "How often to rebuild preconditioner", "SNESSetLagPreconditioner", snes->lagpreconditioner, &lag, &flg)); if (flg) { PetscCheck(lag != -1, PetscObjectComm((PetscObject)snes), PETSC_ERR_USER, "Cannot set the lag to -1 from the command line since the preconditioner must be built as least once, perhaps you mean -2"); PetscCall(SNESSetLagPreconditioner(snes, lag)); } PetscCall(PetscOptionsBool("-snes_lag_preconditioner_persists", "Preconditioner lagging through multiple SNES solves", "SNESSetLagPreconditionerPersists", snes->lagjac_persist, &persist, &flg)); if (flg) PetscCall(SNESSetLagPreconditionerPersists(snes, persist)); PetscCall(PetscOptionsInt("-snes_lag_jacobian", "How often to rebuild Jacobian", "SNESSetLagJacobian", snes->lagjacobian, &lag, &flg)); if (flg) { PetscCheck(lag != -1, PetscObjectComm((PetscObject)snes), PETSC_ERR_USER, "Cannot set the lag to -1 from the command line since the Jacobian must be built as least once, perhaps you mean -2"); PetscCall(SNESSetLagJacobian(snes, lag)); } PetscCall(PetscOptionsBool("-snes_lag_jacobian_persists", "Jacobian lagging through multiple SNES solves", "SNESSetLagJacobianPersists", snes->lagjac_persist, &persist, &flg)); if (flg) PetscCall(SNESSetLagJacobianPersists(snes, persist)); PetscCall(PetscOptionsInt("-snes_grid_sequence", "Use grid sequencing to generate initial guess", "SNESSetGridSequence", snes->gridsequence, &grids, &flg)); if (flg) PetscCall(SNESSetGridSequence(snes, grids)); PetscCall(PetscOptionsEList("-snes_convergence_test", "Convergence test", "SNESSetConvergenceTest", convtests, PETSC_STATIC_ARRAY_LENGTH(convtests), "default", &indx, &flg)); if (flg) { switch (indx) { case 0: PetscCall(SNESSetConvergenceTest(snes, SNESConvergedDefault, NULL, NULL)); break; case 1: PetscCall(SNESSetConvergenceTest(snes, SNESConvergedSkip, NULL, NULL)); break; case 2: PetscCall(SNESSetConvergenceTest(snes, SNESConvergedCorrectPressure, NULL, NULL)); break; } } PetscCall(PetscOptionsEList("-snes_norm_schedule", "SNES Norm schedule", "SNESSetNormSchedule", SNESNormSchedules, 5, "function", &indx, &flg)); if (flg) PetscCall(SNESSetNormSchedule(snes, (SNESNormSchedule)indx)); PetscCall(PetscOptionsEList("-snes_function_type", "SNES Norm schedule", "SNESSetFunctionType", SNESFunctionTypes, 2, "unpreconditioned", &indx, &flg)); if (flg) PetscCall(SNESSetFunctionType(snes, (SNESFunctionType)indx)); kctx = (SNESKSPEW *)snes->kspconvctx; PetscCall(PetscOptionsBool("-snes_ksp_ew", "Use Eisentat-Walker linear system convergence test", "SNESKSPSetUseEW", snes->ksp_ewconv, &snes->ksp_ewconv, NULL)); PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix)); PetscCall(PetscSNPrintf(ewprefix, sizeof(ewprefix), "%s%s", optionsprefix ? optionsprefix : "", "snes_")); PetscCall(SNESEWSetFromOptions_Private(kctx, PETSC_TRUE, PetscObjectComm((PetscObject)snes), ewprefix)); flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_monitor_cancel", "Remove all monitors", "SNESMonitorCancel", flg, &flg, &set)); if (set && flg) PetscCall(SNESMonitorCancel(snes)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor", "Monitor norm of function", "SNESMonitorDefault", SNESMonitorDefault, SNESMonitorDefaultSetUp)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_short", "Monitor norm of function with fewer digits", "SNESMonitorDefaultShort", SNESMonitorDefaultShort, NULL)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_range", "Monitor range of elements of function", "SNESMonitorRange", SNESMonitorRange, NULL)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_ratio", "Monitor ratios of the norm of function for consecutive steps", "SNESMonitorRatio", SNESMonitorRatio, SNESMonitorRatioSetUp)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_field", "Monitor norm of function (split into fields)", "SNESMonitorDefaultField", SNESMonitorDefaultField, NULL)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_solution", "View solution at each iteration", "SNESMonitorSolution", SNESMonitorSolution, NULL)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_solution_update", "View correction at each iteration", "SNESMonitorSolutionUpdate", SNESMonitorSolutionUpdate, NULL)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_residual", "View residual at each iteration", "SNESMonitorResidual", SNESMonitorResidual, NULL)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_jacupdate_spectrum", "Print the change in the spectrum of the Jacobian", "SNESMonitorJacUpdateSpectrum", SNESMonitorJacUpdateSpectrum, NULL)); PetscCall(SNESMonitorSetFromOptions(snes, "-snes_monitor_fields", "Monitor norm of function per field", "SNESMonitorSet", SNESMonitorFields, NULL)); PetscCall(PetscOptionsBool("-snes_monitor_pause_final", "Pauses all draw monitors at the final iterate", "SNESMonitorPauseFinal_Internal", PETSC_FALSE, &snes->pauseFinal, NULL)); PetscCall(PetscOptionsString("-snes_monitor_python", "Use Python function", "SNESMonitorSet", NULL, monfilename, sizeof(monfilename), &flg)); if (flg) PetscCall(PetscPythonMonitorSet((PetscObject)snes, monfilename)); flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_monitor_lg_range", "Plot function range at each iteration", "SNESMonitorLGRange", flg, &flg, NULL)); if (flg) { PetscViewer ctx; PetscCall(PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes), NULL, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &ctx)); PetscCall(SNESMonitorSet(snes, SNESMonitorLGRange, ctx, (PetscCtxDestroyFn *)PetscViewerDestroy)); } PetscCall(PetscViewerDestroy(&snes->convergedreasonviewer)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_converged_reason", &snes->convergedreasonviewer, &snes->convergedreasonformat, NULL)); flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_converged_reason_view_cancel", "Remove all converged reason viewers", "SNESConvergedReasonViewCancel", flg, &flg, &set)); if (set && flg) PetscCall(SNESConvergedReasonViewCancel(snes)); flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_fd", "Use finite differences (slow) to compute Jacobian", "SNESComputeJacobianDefault", flg, &flg, NULL)); if (flg) { void *functx; DM dm; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESUnsetJacobianContext_Internal(dm)); PetscCall(SNESGetFunction(snes, NULL, NULL, &functx)); PetscCall(SNESSetJacobian(snes, snes->jacobian, snes->jacobian_pre, SNESComputeJacobianDefault, functx)); PetscCall(PetscInfo(snes, "Setting default finite difference Jacobian matrix\n")); } flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_fd_function", "Use finite differences (slow) to compute function from user objective", "SNESObjectiveComputeFunctionDefaultFD", flg, &flg, NULL)); if (flg) PetscCall(SNESSetFunction(snes, NULL, SNESObjectiveComputeFunctionDefaultFD, NULL)); flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_fd_color", "Use finite differences with coloring to compute Jacobian", "SNESComputeJacobianDefaultColor", flg, &flg, NULL)); if (flg) { DM dm; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESUnsetJacobianContext_Internal(dm)); PetscCall(SNESSetJacobian(snes, snes->jacobian, snes->jacobian_pre, SNESComputeJacobianDefaultColor, NULL)); PetscCall(PetscInfo(snes, "Setting default finite difference coloring Jacobian matrix\n")); } flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_mf_operator", "Use a Matrix-Free Jacobian with user-provided matrix for computing the preconditioner", "SNESSetUseMatrixFree", PETSC_FALSE, &snes->mf_operator, &flg)); if (flg && snes->mf_operator) { snes->mf_operator = PETSC_TRUE; snes->mf = PETSC_TRUE; } flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_mf", "Use a Matrix-Free Jacobian with no matrix for computing the preconditioner", "SNESSetUseMatrixFree", PETSC_FALSE, &snes->mf, &flg)); if (!flg && snes->mf_operator) snes->mf = PETSC_TRUE; PetscCall(PetscOptionsInt("-snes_mf_version", "Matrix-Free routines version 1 or 2", "None", snes->mf_version, &snes->mf_version, NULL)); PetscCall(PetscOptionsName("-snes_test_function", "Compare hand-coded and finite difference functions", "None", &snes->testFunc)); PetscCall(PetscOptionsName("-snes_test_jacobian", "Compare hand-coded and finite difference Jacobians", "None", &snes->testJac)); flg = PETSC_FALSE; PetscCall(SNESGetNPCSide(snes, &pcside)); PetscCall(PetscOptionsEnum("-snes_npc_side", "SNES nonlinear preconditioner side", "SNESSetNPCSide", PCSides, (PetscEnum)pcside, (PetscEnum *)&pcside, &flg)); if (flg) PetscCall(SNESSetNPCSide(snes, pcside)); #if defined(PETSC_HAVE_SAWS) /* Publish convergence information using SAWs */ flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_monitor_saws", "Publish SNES progress using SAWs", "SNESMonitorSet", flg, &flg, NULL)); if (flg) { void *ctx; PetscCall(SNESMonitorSAWsCreate(snes, &ctx)); PetscCall(SNESMonitorSet(snes, SNESMonitorSAWs, ctx, SNESMonitorSAWsDestroy)); } #endif #if defined(PETSC_HAVE_SAWS) { PetscBool set; flg = PETSC_FALSE; PetscCall(PetscOptionsBool("-snes_saws_block", "Block for SAWs at end of SNESSolve", "PetscObjectSAWsBlock", ((PetscObject)snes)->amspublishblock, &flg, &set)); if (set) PetscCall(PetscObjectSAWsSetBlock((PetscObject)snes, flg)); } #endif for (i = 0; i < numberofsetfromoptions; i++) PetscCall((*othersetfromoptions[i])(snes)); PetscTryTypeMethod(snes, setfromoptions, PetscOptionsObject); /* process any options handlers added with PetscObjectAddOptionsHandler() */ PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)snes, PetscOptionsObject)); PetscOptionsEnd(); if (snes->linesearch) { PetscCall(SNESGetLineSearch(snes, &snes->linesearch)); PetscCall(SNESLineSearchSetFromOptions(snes->linesearch)); } if (snes->usesksp) { if (!snes->ksp) PetscCall(SNESGetKSP(snes, &snes->ksp)); PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre)); PetscCall(KSPSetFromOptions(snes->ksp)); } /* if user has set the SNES NPC type via options database, create it. */ PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix)); PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, optionsprefix, "-npc_snes_type", &pcset)); if (pcset && (!snes->npc)) PetscCall(SNESGetNPC(snes, &snes->npc)); if (snes->npc) PetscCall(SNESSetFromOptions(snes->npc)); snes->setfromoptionscalled++; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESResetFromOptions - Sets various `SNES` and `KSP` parameters from user options ONLY if the `SNESSetFromOptions()` was previously called Collective Input Parameter: . snes - the `SNES` context Level: advanced .seealso: [](ch_snes), `SNES`, `SNESSetFromOptions()`, `SNESSetOptionsPrefix()` @*/ PetscErrorCode SNESResetFromOptions(SNES snes) { PetscFunctionBegin; if (snes->setfromoptionscalled) PetscCall(SNESSetFromOptions(snes)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for the nonlinear solvers. Logically Collective; No Fortran Support Input Parameters: + snes - the `SNES` context . compute - function to compute the context - destroy - function to destroy the context, see `PetscCtxDestroyFn` for the calling sequence Calling sequence of `compute`: + snes - the `SNES` context - ctx - context to be computed Level: intermediate Note: This routine is useful if you are performing grid sequencing or using `SNESFAS` and need the appropriate context generated for each level. Use `SNESSetApplicationContext()` to see the context immediately .seealso: [](ch_snes), `SNESGetApplicationContext()`, `SNESSetApplicationContext()`, `PetscCtxDestroyFn` @*/ PetscErrorCode SNESSetComputeApplicationContext(SNES snes, PetscErrorCode (*compute)(SNES snes, void **ctx), PetscCtxDestroyFn *destroy) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->ops->usercompute = compute; snes->ops->ctxdestroy = destroy; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetApplicationContext - Sets the optional user-defined context for the nonlinear solvers. Logically Collective Input Parameters: + snes - the `SNES` context - ctx - the user context Level: intermediate Notes: Users can provide a context when constructing the `SNES` options and then access it inside their function, Jacobian computation, or other evaluation function with `SNESGetApplicationContext()` To provide a function that computes the context for you use `SNESSetComputeApplicationContext()` Fortran Note: This only works when `ctx` is a Fortran derived type (it cannot be a `PetscObject`), we recommend writing a Fortran interface definition for this function that tells the Fortran compiler the derived data type that is passed in as the `ctx` argument. See `SNESGetApplicationContext()` for an example. .seealso: [](ch_snes), `SNES`, `SNESSetComputeApplicationContext()`, `SNESGetApplicationContext()` @*/ PetscErrorCode SNESSetApplicationContext(SNES snes, void *ctx) { KSP ksp; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(KSPSetApplicationContext(ksp, ctx)); snes->ctx = ctx; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetApplicationContext - Gets the user-defined context for the nonlinear solvers set with `SNESGetApplicationContext()` or `SNESSetComputeApplicationContext()` Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . ctx - user context Level: intermediate Fortran Notes: This only works when the context is a Fortran derived type (it cannot be a `PetscObject`) and you **must** write a Fortran interface definition for this function that tells the Fortran compiler the derived data type that is returned as the `ctx` argument. For example, .vb Interface SNESGetApplicationContext Subroutine SNESGetApplicationContext(snes,ctx,ierr) #include use petscsnes SNES snes type(tUsertype), pointer :: ctx PetscErrorCode ierr End Subroutine End Interface SNESGetApplicationContext .ve The prototype for `ctx` must be .vb type(tUsertype), pointer :: ctx .ve .seealso: [](ch_snes), `SNESSetApplicationContext()`, `SNESSetComputeApplicationContext()` @*/ PetscErrorCode SNESGetApplicationContext(SNES snes, PeCtx ctx) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *(void **)ctx = snes->ctx; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetUseMatrixFree - indicates that `SNES` should use matrix-free finite difference matrix-vector products to apply the Jacobian. Logically Collective Input Parameters: + snes - `SNES` context . mf_operator - use matrix-free only for the Amat used by `SNESSetJacobian()`, this means the user provided Pmat will continue to be used - mf - use matrix-free for both the Amat and Pmat used by `SNESSetJacobian()`, both the Amat and Pmat set in `SNESSetJacobian()` will be ignored. With this option no matrix-element based preconditioners can be used in the linear solve since the matrix won't be explicitly available Options Database Keys: + -snes_mf_operator - use matrix-free only for the mat operator . -snes_mf - use matrix-free for both the mat and pmat operator . -snes_fd_color - compute the Jacobian via coloring and finite differences. - -snes_fd - compute the Jacobian via finite differences (slow) Level: intermediate Note: `SNES` supports three approaches for computing (approximate) Jacobians: user provided via `SNESSetJacobian()`, matrix-free using `MatCreateSNESMF()`, and computing explicitly with finite differences and coloring using `MatFDColoring`. It is also possible to use automatic differentiation and the `MatFDColoring` object. .seealso: [](ch_snes), `SNES`, `SNESGetUseMatrixFree()`, `MatCreateSNESMF()`, `SNESComputeJacobianDefaultColor()`, `MatFDColoring` @*/ PetscErrorCode SNESSetUseMatrixFree(SNES snes, PetscBool mf_operator, PetscBool mf) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveBool(snes, mf_operator, 2); PetscValidLogicalCollectiveBool(snes, mf, 3); snes->mf = mf_operator ? PETSC_TRUE : mf; snes->mf_operator = mf_operator; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetUseMatrixFree - indicates if the `SNES` uses matrix-free finite difference matrix vector products to apply the Jacobian. Not Collective, but the resulting flags will be the same on all MPI processes Input Parameter: . snes - `SNES` context Output Parameters: + mf_operator - use matrix-free only for the Amat used by `SNESSetJacobian()`, this means the user provided Pmat will continue to be used - mf - use matrix-free for both the Amat and Pmat used by `SNESSetJacobian()`, both the Amat and Pmat set in `SNESSetJacobian()` will be ignored Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESSetUseMatrixFree()`, `MatCreateSNESMF()` @*/ PetscErrorCode SNESGetUseMatrixFree(SNES snes, PetscBool *mf_operator, PetscBool *mf) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (mf) *mf = snes->mf; if (mf_operator) *mf_operator = snes->mf_operator; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetIterationNumber - Gets the number of nonlinear iterations completed in the current or most recent `SNESSolve()` Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . iter - iteration number Level: intermediate Notes: For example, during the computation of iteration 2 this would return 1. This is useful for using lagged Jacobians (where one does not recompute the Jacobian at each `SNES` iteration). For example, the code .vb ierr = SNESGetIterationNumber(snes,&it); if (!(it % 2)) { [compute Jacobian here] } .ve can be used in your function that computes the Jacobian to cause the Jacobian to be recomputed every second `SNES` iteration. See also `SNESSetLagJacobian()` After the `SNES` solve is complete this will return the number of nonlinear iterations used. .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESSetLagJacobian()`, `SNESGetLinearSolveIterations()`, `SNESSetMonitor()` @*/ PetscErrorCode SNESGetIterationNumber(SNES snes, PetscInt *iter) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(iter, 2); *iter = snes->iter; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetIterationNumber - Sets the current iteration number. Not Collective Input Parameters: + snes - `SNES` context - iter - iteration number Level: developer Note: This should only be called inside a `SNES` nonlinear solver. .seealso: [](ch_snes), `SNESGetLinearSolveIterations()` @*/ PetscErrorCode SNESSetIterationNumber(SNES snes, PetscInt iter) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); snes->iter = iter; PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps attempted by the nonlinear solver in the current or most recent `SNESSolve()` . Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . nfails - number of unsuccessful steps attempted Level: intermediate Note: This counter is reset to zero for each successive call to `SNESSolve()`. .seealso: [](ch_snes), `SNES`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`, `SNESSetMaxNonlinearStepFailures()`, `SNESGetMaxNonlinearStepFailures()` @*/ PetscErrorCode SNESGetNonlinearStepFailures(SNES snes, PetscInt *nfails) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(nfails, 2); *nfails = snes->numFailures; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps attempted by the nonlinear solver before it gives up and returns unconverged or generates an error Not Collective Input Parameters: + snes - `SNES` context - maxFails - maximum of unsuccessful steps allowed, use `PETSC_UNLIMITED` to have no limit on the number of failures Options Database Key: . -snes_max_fail - maximum number of unsuccessful steps allowed Level: intermediate Developer Note: The options database key is wrong for this function name .seealso: [](ch_snes), `SNESSetErrorIfNotConverged()`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`, `SNESGetMaxNonlinearStepFailures()`, `SNESGetNonlinearStepFailures()` @*/ PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (maxFails == PETSC_UNLIMITED) { snes->maxFailures = PETSC_INT_MAX; } else { PetscCheck(maxFails >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures"); snes->maxFailures = maxFails; } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps attempted by the nonlinear solver before it gives up and returns unconverged or generates an error Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . maxFails - maximum of unsuccessful steps Level: intermediate .seealso: [](ch_snes), `SNESSetErrorIfNotConverged()`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`, `SNESSetMaxNonlinearStepFailures()`, `SNESGetNonlinearStepFailures()` @*/ PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(maxFails, 2); *maxFails = snes->maxFailures; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations done by the `SNES` object in the current or most recent `SNESSolve()` Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . nfuncs - number of evaluations Level: intermediate Note: Reset every time `SNESSolve()` is called unless `SNESSetCountersReset()` is used. .seealso: [](ch_snes), `SNES`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, `SNESGetLinearSolveFailures()`, `SNESSetCountersReset()` @*/ PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(nfuncs, 2); *nfuncs = snes->nfuncs; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetLinearSolveFailures - Gets the number of failed (non-converged) linear solvers in the current or most recent `SNESSolve()` Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . nfails - number of failed solves Options Database Key: . -snes_max_linear_solve_fail - The number of failures before the solve is terminated Level: intermediate Note: This counter is reset to zero for each successive call to `SNESSolve()`. .seealso: [](ch_snes), `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()` @*/ PetscErrorCode SNESGetLinearSolveFailures(SNES snes, PetscInt *nfails) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(nfails, 2); *nfails = snes->numLinearSolveFailures; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts allowed before `SNES` returns with a diverged reason of `SNES_DIVERGED_LINEAR_SOLVE` Logically Collective Input Parameters: + snes - `SNES` context - maxFails - maximum allowed linear solve failures, use `PETSC_UNLIMITED` to have no limit on the number of failures Options Database Key: . -snes_max_linear_solve_fail - The number of failures before the solve is terminated Level: intermediate Note: By default this is 0; that is `SNES` returns on the first failed linear solve Developer Note: The options database key is wrong for this function name .seealso: [](ch_snes), `SNESSetErrorIfNotConverged()`, `SNESGetLinearSolveFailures()`, `SNESGetMaxLinearSolveFailures()`, `SNESGetLinearSolveIterations()` @*/ PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveInt(snes, maxFails, 2); if (maxFails == PETSC_UNLIMITED) { snes->maxLinearSolveFailures = PETSC_INT_MAX; } else { PetscCheck(maxFails >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a negative maximum number of failures"); snes->maxLinearSolveFailures = maxFails; } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that are allowed before `SNES` returns as unsuccessful Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . maxFails - maximum of unsuccessful solves allowed Level: intermediate Note: By default this is 1; that is `SNES` returns on the first failed linear solve .seealso: [](ch_snes), `SNESSetErrorIfNotConverged()`, `SNESGetLinearSolveFailures()`, `SNESGetLinearSolveIterations()`, `SNESSetMaxLinearSolveFailures()`, @*/ PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(maxFails, 2); *maxFails = snes->maxLinearSolveFailures; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetLinearSolveIterations - Gets the total number of linear iterations used by the nonlinear solver in the most recent `SNESSolve()` Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . lits - number of linear iterations Level: intermediate Notes: This counter is reset to zero for each successive call to `SNESSolve()` unless `SNESSetCountersReset()` is used. If the linear solver fails inside the `SNESSolve()` the iterations for that call to the linear solver are not included. If you wish to count them then call `KSPGetIterationNumber()` after the failed solve. .seealso: [](ch_snes), `SNES`, `SNESGetIterationNumber()`, `SNESGetLinearSolveFailures()`, `SNESGetMaxLinearSolveFailures()`, `SNESSetCountersReset()` @*/ PetscErrorCode SNESGetLinearSolveIterations(SNES snes, PetscInt *lits) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(lits, 2); *lits = snes->linear_its; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetCountersReset - Sets whether or not the counters for linear iterations and function evaluations are reset every time `SNESSolve()` is called. Logically Collective Input Parameters: + snes - `SNES` context - reset - whether to reset the counters or not, defaults to `PETSC_TRUE` Level: developer .seealso: [](ch_snes), `SNESGetNumberFunctionEvals()`, `SNESGetLinearSolveIterations()`, `SNESGetNPC()` @*/ PetscErrorCode SNESSetCountersReset(SNES snes, PetscBool reset) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveBool(snes, reset, 2); snes->counters_reset = reset; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESResetCounters - Reset counters for linear iterations and function evaluations. Logically Collective Input Parameters: . snes - `SNES` context Level: developer Note: It honors the flag set with `SNESSetCountersReset()` .seealso: [](ch_snes), `SNESGetNumberFunctionEvals()`, `SNESGetLinearSolveIterations()`, `SNESGetNPC()` @*/ PetscErrorCode SNESResetCounters(SNES snes) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (snes->counters_reset) { snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetKSP - Sets a `KSP` context for the `SNES` object to use Not Collective, but the `SNES` and `KSP` objects must live on the same `MPI_Comm` Input Parameters: + snes - the `SNES` context - ksp - the `KSP` context Level: developer Notes: The `SNES` object already has its `KSP` object, you can obtain with `SNESGetKSP()` so this routine is rarely needed. The `KSP` object that is already in the `SNES` object has its reference count decreased by one when this is called. .seealso: [](ch_snes), `SNES`, `KSP`, `KSPGetPC()`, `SNESCreate()`, `KSPCreate()` @*/ PetscErrorCode SNESSetKSP(SNES snes, KSP ksp) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(ksp, KSP_CLASSID, 2); PetscCheckSameComm(snes, 1, ksp, 2); PetscCall(PetscObjectReference((PetscObject)ksp)); if (snes->ksp) PetscCall(PetscObjectDereference((PetscObject)snes->ksp)); snes->ksp = ksp; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESParametersInitialize - Sets all the parameters in `snes` to their default value (when `SNESCreate()` was called) if they currently contain default values Collective Input Parameter: . snes - the `SNES` object Level: developer Developer Note: This is called by all the `SNESCreate_XXX()` routines. .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESDestroy()`, `SNESSetLagPreconditioner()`, `SNESSetLagJacobian()`, `PetscObjectParameterSetDefault()` @*/ PetscErrorCode SNESParametersInitialize(SNES snes) { PetscObjectParameterSetDefault(snes, max_its, 50); PetscObjectParameterSetDefault(snes, max_funcs, 10000); PetscObjectParameterSetDefault(snes, rtol, PetscDefined(USE_REAL_SINGLE) ? 1.e-5 : 1.e-8); PetscObjectParameterSetDefault(snes, abstol, PetscDefined(USE_REAL_SINGLE) ? 1.e-25 : 1.e-50); PetscObjectParameterSetDefault(snes, stol, PetscDefined(USE_REAL_SINGLE) ? 1.e-5 : 1.e-8); PetscObjectParameterSetDefault(snes, divtol, 1.e4); return PETSC_SUCCESS; } /*@ SNESCreate - Creates a nonlinear solver context used to manage a set of nonlinear solves Collective Input Parameter: . comm - MPI communicator Output Parameter: . outsnes - the new `SNES` context Options Database Keys: + -snes_mf - Activates default matrix-free Jacobian-vector products, and no matrix to construct a preconditioner . -snes_mf_operator - Activates default matrix-free Jacobian-vector products, and a user-provided matrix as set by `SNESSetJacobian()` . -snes_fd_coloring - uses a relative fast computation of the Jacobian using finite differences and a graph coloring - -snes_fd - Uses (slow!) finite differences to compute Jacobian Level: beginner Developer Notes: `SNES` always creates a `KSP` object even though many `SNES` methods do not use it. This is unfortunate and should be fixed at some point. The flag snes->usesksp indicates if the particular method does use `KSP` and regulates if the information about the `KSP` is printed in `SNESView()`. `TSSetFromOptions()` does call `SNESSetFromOptions()` which can lead to users being confused by help messages about meaningless `SNES` options. `SNES` always creates the `snes->kspconvctx` even though it is used by only one type. This should be fixed. .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESDestroy()`, `SNESSetLagPreconditioner()`, `SNESSetLagJacobian()` @*/ PetscErrorCode SNESCreate(MPI_Comm comm, SNES *outsnes) { SNES snes; SNESKSPEW *kctx; PetscFunctionBegin; PetscAssertPointer(outsnes, 2); PetscCall(SNESInitializePackage()); PetscCall(PetscHeaderCreate(snes, SNES_CLASSID, "SNES", "Nonlinear solver", "SNES", comm, SNESDestroy, SNESView)); snes->ops->converged = SNESConvergedDefault; snes->usesksp = PETSC_TRUE; snes->norm = 0.0; snes->xnorm = 0.0; snes->ynorm = 0.0; snes->normschedule = SNES_NORM_ALWAYS; snes->functype = SNES_FUNCTION_DEFAULT; snes->ttol = 0.0; snes->rnorm0 = 0; snes->nfuncs = 0; snes->numFailures = 0; snes->maxFailures = 1; snes->linear_its = 0; snes->lagjacobian = 1; snes->jac_iter = 0; snes->lagjac_persist = PETSC_FALSE; snes->lagpreconditioner = 1; snes->pre_iter = 0; snes->lagpre_persist = PETSC_FALSE; snes->numbermonitors = 0; snes->numberreasonviews = 0; snes->data = NULL; snes->setupcalled = PETSC_FALSE; snes->ksp_ewconv = PETSC_FALSE; snes->nwork = 0; snes->work = NULL; snes->nvwork = 0; snes->vwork = NULL; snes->conv_hist_len = 0; snes->conv_hist_max = 0; snes->conv_hist = NULL; snes->conv_hist_its = NULL; snes->conv_hist_reset = PETSC_TRUE; snes->counters_reset = PETSC_TRUE; snes->vec_func_init_set = PETSC_FALSE; snes->reason = SNES_CONVERGED_ITERATING; snes->npcside = PC_RIGHT; snes->setfromoptionscalled = 0; snes->mf = PETSC_FALSE; snes->mf_operator = PETSC_FALSE; snes->mf_version = 1; snes->numLinearSolveFailures = 0; snes->maxLinearSolveFailures = 1; snes->vizerotolerance = 1.e-8; snes->checkjacdomainerror = PetscDefined(USE_DEBUG) ? PETSC_TRUE : PETSC_FALSE; /* Set this to true if the implementation of SNESSolve_XXX does compute the residual at the final solution. */ snes->alwayscomputesfinalresidual = PETSC_FALSE; /* Create context to compute Eisenstat-Walker relative tolerance for KSP */ PetscCall(PetscNew(&kctx)); snes->kspconvctx = kctx; kctx->version = 2; kctx->rtol_0 = 0.3; /* Eisenstat and Walker suggest rtol_0=.5, but this was too large for some test cases */ kctx->rtol_last = 0.0; kctx->rtol_max = 0.9; kctx->gamma = 1.0; kctx->alpha = 0.5 * (1.0 + PetscSqrtReal(5.0)); kctx->alpha2 = kctx->alpha; kctx->threshold = 0.1; kctx->lresid_last = 0.0; kctx->norm_last = 0.0; kctx->rk_last = 0.0; kctx->rk_last_2 = 0.0; kctx->rtol_last_2 = 0.0; kctx->v4_p1 = 0.1; kctx->v4_p2 = 0.4; kctx->v4_p3 = 0.7; kctx->v4_m1 = 0.8; kctx->v4_m2 = 0.5; kctx->v4_m3 = 0.1; kctx->v4_m4 = 0.5; PetscCall(SNESParametersInitialize(snes)); *outsnes = snes; PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetFunction - Sets the function evaluation routine and function vector for use by the `SNES` routines in solving systems of nonlinear equations. Logically Collective Input Parameters: + snes - the `SNES` context . r - vector to store function values, may be `NULL` . f - function evaluation routine; for calling sequence see `SNESFunctionFn` - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`) Level: beginner .seealso: [](ch_snes), `SNES`, `SNESGetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`, `SNESSetPicard()`, `SNESFunctionFn` @*/ PetscErrorCode SNESSetFunction(SNES snes, Vec r, SNESFunctionFn *f, void *ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (r) { PetscValidHeaderSpecific(r, VEC_CLASSID, 2); PetscCheckSameComm(snes, 1, r, 2); PetscCall(PetscObjectReference((PetscObject)r)); PetscCall(VecDestroy(&snes->vec_func)); snes->vec_func = r; } PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESSetFunction(dm, f, ctx)); if (f == SNESPicardComputeFunction) PetscCall(DMSNESSetMFFunction(dm, SNESPicardComputeMFFunction, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetInitialFunction - Set an already computed function evaluation at the initial guess to be reused by `SNESSolve()`. Logically Collective Input Parameters: + snes - the `SNES` context - f - vector to store function value Level: developer Notes: This should not be modified during the solution procedure. This is used extensively in the `SNESFAS` hierarchy and in nonlinear preconditioning. .seealso: [](ch_snes), `SNES`, `SNESFAS`, `SNESSetFunction()`, `SNESComputeFunction()`, `SNESSetInitialFunctionNorm()` @*/ PetscErrorCode SNESSetInitialFunction(SNES snes, Vec f) { Vec vec_func; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(f, VEC_CLASSID, 2); PetscCheckSameComm(snes, 1, f, 2); if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { snes->vec_func_init_set = PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(SNESGetFunction(snes, &vec_func, NULL, NULL)); PetscCall(VecCopy(f, vec_func)); snes->vec_func_init_set = PETSC_TRUE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetNormSchedule - Sets the `SNESNormSchedule` used in convergence and monitoring of the `SNES` method, when norms are computed in the solving process Logically Collective Input Parameters: + snes - the `SNES` context - normschedule - the frequency of norm computation Options Database Key: . -snes_norm_schedule - set the schedule Level: advanced Notes: Only certain `SNES` methods support certain `SNESNormSchedules`. Most require evaluation of the nonlinear function and the taking of its norm at every iteration to even ensure convergence at all. However, methods such as custom Gauss-Seidel methods `SNESNGS` and the like do not require the norm of the function to be computed, and therefore may either be monitored for convergence or not. As these are often used as nonlinear preconditioners, monitoring the norm of their error is not a useful enterprise within their solution. .seealso: [](ch_snes), `SNESNormSchedule`, `SNESGetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()` @*/ PetscErrorCode SNESSetNormSchedule(SNES snes, SNESNormSchedule normschedule) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->normschedule = normschedule; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetNormSchedule - Gets the `SNESNormSchedule` used in convergence and monitoring of the `SNES` method. Logically Collective Input Parameters: + snes - the `SNES` context - normschedule - the type of the norm used Level: advanced .seealso: [](ch_snes), `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule` @*/ PetscErrorCode SNESGetNormSchedule(SNES snes, SNESNormSchedule *normschedule) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *normschedule = snes->normschedule; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetFunctionNorm - Sets the last computed residual norm. Logically Collective Input Parameters: + snes - the `SNES` context - norm - the value of the norm Level: developer .seealso: [](ch_snes), `SNES`, `SNESGetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule` @*/ PetscErrorCode SNESSetFunctionNorm(SNES snes, PetscReal norm) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->norm = norm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetFunctionNorm - Gets the last computed norm of the residual Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . norm - the last computed residual norm Level: developer .seealso: [](ch_snes), `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule` @*/ PetscErrorCode SNESGetFunctionNorm(SNES snes, PetscReal *norm) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(norm, 2); *norm = snes->norm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetUpdateNorm - Gets the last computed norm of the solution update Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . ynorm - the last computed update norm Level: developer Note: The new solution is the current solution plus the update, so this norm is an indication of the size of the update .seealso: [](ch_snes), `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `SNESGetFunctionNorm()` @*/ PetscErrorCode SNESGetUpdateNorm(SNES snes, PetscReal *ynorm) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(ynorm, 2); *ynorm = snes->ynorm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetSolutionNorm - Gets the last computed norm of the solution Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . xnorm - the last computed solution norm Level: developer .seealso: [](ch_snes), `SNES`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `SNESGetFunctionNorm()`, `SNESGetUpdateNorm()` @*/ PetscErrorCode SNESGetSolutionNorm(SNES snes, PetscReal *xnorm) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(xnorm, 2); *xnorm = snes->xnorm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetFunctionType - Sets the `SNESFunctionType` of the `SNES` method. Logically Collective Input Parameters: + snes - the `SNES` context - type - the function type Level: developer Values of the function type\: + `SNES_FUNCTION_DEFAULT` - the default for the given `SNESType` . `SNES_FUNCTION_UNPRECONDITIONED` - an unpreconditioned function evaluation (this is the function provided with `SNESSetFunction()` - `SNES_FUNCTION_PRECONDITIONED` - a transformation of the function provided with `SNESSetFunction()` Note: Different `SNESType`s use this value in different ways .seealso: [](ch_snes), `SNES`, `SNESFunctionType`, `SNESGetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule` @*/ PetscErrorCode SNESSetFunctionType(SNES snes, SNESFunctionType type) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->functype = type; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetFunctionType - Gets the `SNESFunctionType` used in convergence and monitoring set with `SNESSetFunctionType()` of the SNES method. Logically Collective Input Parameters: + snes - the `SNES` context - type - the type of the function evaluation, see `SNESSetFunctionType()` Level: advanced .seealso: [](ch_snes), `SNESSetFunctionType()`, `SNESFunctionType`, `SNESSetNormSchedule()`, `SNESComputeFunction()`, `VecNorm()`, `SNESSetFunction()`, `SNESSetInitialFunction()`, `SNESNormSchedule` @*/ PetscErrorCode SNESGetFunctionType(SNES snes, SNESFunctionType *type) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *type = snes->functype; PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetNGS - Sets the user nonlinear Gauss-Seidel routine for use with composed nonlinear solvers. Input Parameters: + snes - the `SNES` context, usually of the `SNESType` `SNESNGS` . f - function evaluation routine to apply Gauss-Seidel, see `SNESNGSFn` for calling sequence - ctx - [optional] user-defined context for private data for the smoother evaluation routine (may be `NULL`) Level: intermediate Note: The `SNESNGS` routines are used by the composed nonlinear solver to generate a problem appropriate update to the solution, particularly `SNESFAS`. .seealso: [](ch_snes), `SNESNGS`, `SNESGetNGS()`, `SNESNCG`, `SNESGetFunction()`, `SNESComputeNGS()`, `SNESNGSFn` @*/ PetscErrorCode SNESSetNGS(SNES snes, SNESNGSFn *f, void *ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESSetNGS(dm, f, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } /* This is used for -snes_mf_operator; it uses a duplicate of snes->jacobian_pre because snes->jacobian_pre cannot be changed during the KSPSolve() */ PetscErrorCode SNESPicardComputeMFFunction(SNES snes, Vec x, Vec f, void *ctx) { DM dm; DMSNES sdm; PetscFunctionBegin; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); /* A(x)*x - b(x) */ if (sdm->ops->computepfunction) { PetscCallBack("SNES Picard callback function", (*sdm->ops->computepfunction)(snes, x, f, sdm->pctx)); PetscCall(VecScale(f, -1.0)); /* Cannot share nonzero pattern because of the possible use of SNESComputeJacobianDefault() */ if (!snes->picard) PetscCall(MatDuplicate(snes->jacobian_pre, MAT_DO_NOT_COPY_VALUES, &snes->picard)); PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->picard, snes->picard, sdm->pctx)); PetscCall(MatMultAdd(snes->picard, x, f, f)); } else { PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->picard, snes->picard, sdm->pctx)); PetscCall(MatMult(snes->picard, x, f)); } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode SNESPicardComputeFunction(SNES snes, Vec x, Vec f, void *ctx) { DM dm; DMSNES sdm; PetscFunctionBegin; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); /* A(x)*x - b(x) */ if (sdm->ops->computepfunction) { PetscCallBack("SNES Picard callback function", (*sdm->ops->computepfunction)(snes, x, f, sdm->pctx)); PetscCall(VecScale(f, -1.0)); PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->jacobian, snes->jacobian_pre, sdm->pctx)); PetscCall(MatMultAdd(snes->jacobian_pre, x, f, f)); } else { PetscCallBack("SNES Picard callback Jacobian", (*sdm->ops->computepjacobian)(snes, x, snes->jacobian, snes->jacobian_pre, sdm->pctx)); PetscCall(MatMult(snes->jacobian_pre, x, f)); } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode SNESPicardComputeJacobian(SNES snes, Vec x1, Mat J, Mat B, void *ctx) { PetscFunctionBegin; /* the jacobian matrix should be pre-filled in SNESPicardComputeFunction */ /* must assembly if matrix-free to get the last SNES solution */ PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetPicard - Use `SNES` to solve the system $A(x) x = bp(x) + b $ via a Picard type iteration (Picard linearization) Logically Collective Input Parameters: + snes - the `SNES` context . r - vector to store function values, may be `NULL` . bp - function evaluation routine, may be `NULL`, for the calling sequence see `SNESFunctionFn` . Amat - matrix with which $A(x) x - bp(x) - b$ is to be computed . Pmat - matrix from which preconditioner is computed (usually the same as `Amat`) . J - function to compute matrix values, for the calling sequence see `SNESJacobianFn` - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`) Level: intermediate Notes: It is often better to provide the nonlinear function $F()$ and some approximation to its Jacobian directly and use an approximate Newton solver. This interface is provided to allow porting/testing a previous Picard based code in PETSc before converting it to approximate Newton. One can call `SNESSetPicard()` or `SNESSetFunction()` (and possibly `SNESSetJacobian()`) but cannot call both Solves the equation $A(x) x = bp(x) - b$ via the defect correction algorithm $A(x^{n}) (x^{n+1} - x^{n}) = bp(x^{n}) + b - A(x^{n})x^{n}$. When an exact solver is used this corresponds to the "classic" Picard $A(x^{n}) x^{n+1} = bp(x^{n}) + b$ iteration. Run with `-snes_mf_operator` to solve the system with Newton's method using $A(x^{n})$ to construct the preconditioner. We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then the direct Picard iteration $A(x^n) x^{n+1} = bp(x^n) + b$ There is some controversity over the definition of a Picard iteration for nonlinear systems but almost everyone agrees that it involves a linear solve and some believe it is the iteration $A(x^{n}) x^{n+1} = b(x^{n})$ hence we use the name Picard. If anyone has an authoritative reference that defines the Picard iteration different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument \:-). When used with `-snes_mf_operator` this will run matrix-free Newton's method where the matrix-vector product is of the true Jacobian of $A(x)x - bp(x) - b$ and $A(x^{n})$ is used to build the preconditioner When used with `-snes_fd` this will compute the true Jacobian (very slowly one column at a time) and thus represent Newton's method. When used with `-snes_fd_coloring` this will compute the Jacobian via coloring and thus represent a faster implementation of Newton's method. But the the nonzero structure of the Jacobian is, in general larger than that of the Picard matrix $A$ so you must provide in $A$ the needed nonzero structure for the correct coloring. When using `DMDA` this may mean creating the matrix $A$ with `DMCreateMatrix()` using a wider stencil than strictly needed for $A$ or with a `DMDA_STENCIL_BOX`. See the comment in src/snes/tutorials/ex15.c. .seealso: [](ch_snes), `SNES`, `SNESGetFunction()`, `SNESSetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`, `SNESGetPicard()`, `SNESLineSearchPreCheckPicard()`, `SNESFunctionFn`, `SNESJacobianFn` @*/ PetscErrorCode SNESSetPicard(SNES snes, Vec r, SNESFunctionFn *bp, Mat Amat, Mat Pmat, SNESJacobianFn *J, void *ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESSetPicard(dm, bp, J, ctx)); PetscCall(DMSNESSetMFFunction(dm, SNESPicardComputeMFFunction, ctx)); PetscCall(SNESSetFunction(snes, r, SNESPicardComputeFunction, ctx)); PetscCall(SNESSetJacobian(snes, Amat, Pmat, SNESPicardComputeJacobian, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESGetPicard - Returns the context for the Picard iteration Not Collective, but `Vec` is parallel if `SNES` is parallel. Collective if `Vec` is requested, but has not been created yet. Input Parameter: . snes - the `SNES` context Output Parameters: + r - the function (or `NULL`) . f - the function (or `NULL`); for calling sequence see `SNESFunctionFn` . Amat - the matrix used to defined the operation A(x) x - b(x) (or `NULL`) . Pmat - the matrix from which the preconditioner will be constructed (or `NULL`) . J - the function for matrix evaluation (or `NULL`); for calling sequence see `SNESJacobianFn` - ctx - the function context (or `NULL`) Level: advanced .seealso: [](ch_snes), `SNESSetFunction()`, `SNESSetPicard()`, `SNESGetFunction()`, `SNESGetJacobian()`, `SNESGetDM()`, `SNESFunctionFn`, `SNESJacobianFn` @*/ PetscErrorCode SNESGetPicard(SNES snes, Vec *r, SNESFunctionFn **f, Mat *Amat, Mat *Pmat, SNESJacobianFn **J, void **ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(SNESGetFunction(snes, r, NULL, NULL)); PetscCall(SNESGetJacobian(snes, Amat, Pmat, NULL, NULL)); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESGetPicard(dm, f, J, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the nonlinear problem Logically Collective Input Parameters: + snes - the `SNES` context . func - function evaluation routine, see `SNESInitialGuessFn` for the calling sequence - ctx - [optional] user-defined context for private data for the function evaluation routine (may be `NULL`) Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`, `SNESInitialGuessFn` @*/ PetscErrorCode SNESSetComputeInitialGuess(SNES snes, SNESInitialGuessFn *func, void *ctx) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (func) snes->ops->computeinitialguess = func; if (ctx) snes->initialguessP = ctx; PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESGetRhs - Gets the vector for solving F(x) = `rhs`. If `rhs` is not set it assumes a zero right-hand side. Logically Collective Input Parameter: . snes - the `SNES` context Output Parameter: . rhs - the right-hand side vector or `NULL` if there is no right-hand side vector Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESGetSolution()`, `SNESGetFunction()`, `SNESComputeFunction()`, `SNESSetJacobian()`, `SNESSetFunction()` @*/ PetscErrorCode SNESGetRhs(SNES snes, Vec *rhs) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(rhs, 2); *rhs = snes->vec_rhs; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESComputeFunction - Calls the function that has been set with `SNESSetFunction()`. Collective Input Parameters: + snes - the `SNES` context - x - input vector Output Parameter: . y - function vector, as set by `SNESSetFunction()` Level: developer Notes: `SNESComputeFunction()` is typically used within nonlinear solvers implementations, so users would not generally call this routine themselves. When solving for $F(x) = b$, this routine computes $y = F(x) - b$. .seealso: [](ch_snes), `SNES`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESComputeMFFunction()` @*/ PetscErrorCode SNESComputeFunction(SNES snes, Vec x, Vec y) { DM dm; DMSNES sdm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(x, VEC_CLASSID, 2); PetscValidHeaderSpecific(y, VEC_CLASSID, 3); PetscCheckSameComm(snes, 1, x, 2); PetscCheckSameComm(snes, 1, y, 3); PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE)); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); PetscCheck(sdm->ops->computefunction || snes->vec_rhs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); if (sdm->ops->computefunction) { if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) PetscCall(PetscLogEventBegin(SNES_FunctionEval, snes, x, y, 0)); PetscCall(VecLockReadPush(x)); /* ensure domainerror is false prior to computefunction evaluation (may not have been reset) */ snes->domainerror = PETSC_FALSE; { void *ctx; SNESFunctionFn *computefunction; PetscCall(DMSNESGetFunction(dm, &computefunction, &ctx)); PetscCallBack("SNES callback function", (*computefunction)(snes, x, y, ctx)); } PetscCall(VecLockReadPop(x)); if (sdm->ops->computefunction != SNESObjectiveComputeFunctionDefaultFD) PetscCall(PetscLogEventEnd(SNES_FunctionEval, snes, x, y, 0)); } else /* if (snes->vec_rhs) */ { PetscCall(MatMult(snes->jacobian, x, y)); } if (snes->vec_rhs) PetscCall(VecAXPY(y, -1.0, snes->vec_rhs)); snes->nfuncs++; /* domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will propagate the value to all processes */ PetscCall(VecFlag(y, snes->domainerror)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESComputeMFFunction - Calls the function that has been set with `DMSNESSetMFFunction()`. Collective Input Parameters: + snes - the `SNES` context - x - input vector Output Parameter: . y - output vector Level: developer Notes: `SNESComputeMFFunction()` is used within the matrix-vector products called by the matrix created with `MatCreateSNESMF()` so users would not generally call this routine themselves. Since this function is intended for use with finite differencing it does not subtract the right-hand side vector provided with `SNESSolve()` while `SNESComputeFunction()` does. As such, this routine cannot be used with `MatMFFDSetBase()` with a provided F function value even if it applies the same function as `SNESComputeFunction()` if a `SNESSolve()` right-hand side vector is use because the two functions difference would include this right hand side function. .seealso: [](ch_snes), `SNES`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESComputeFunction()`, `MatCreateSNESMF()`, `DMSNESSetMFFunction()` @*/ PetscErrorCode SNESComputeMFFunction(SNES snes, Vec x, Vec y) { DM dm; DMSNES sdm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(x, VEC_CLASSID, 2); PetscValidHeaderSpecific(y, VEC_CLASSID, 3); PetscCheckSameComm(snes, 1, x, 2); PetscCheckSameComm(snes, 1, y, 3); PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE)); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); PetscCall(PetscLogEventBegin(SNES_FunctionEval, snes, x, y, 0)); PetscCall(VecLockReadPush(x)); /* ensure domainerror is false prior to computefunction evaluation (may not have been reset) */ snes->domainerror = PETSC_FALSE; PetscCallBack("SNES callback function", (*sdm->ops->computemffunction)(snes, x, y, sdm->mffunctionctx)); PetscCall(VecLockReadPop(x)); PetscCall(PetscLogEventEnd(SNES_FunctionEval, snes, x, y, 0)); snes->nfuncs++; /* domainerror might not be set on all processes; so we tag vector locally with Inf and the next inner product or norm will propagate the value to all processes */ PetscCall(VecFlag(y, snes->domainerror)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESComputeNGS - Calls the Gauss-Seidel function that has been set with `SNESSetNGS()`. Collective Input Parameters: + snes - the `SNES` context . x - input vector - b - rhs vector Output Parameter: . x - new solution vector Level: developer Note: `SNESComputeNGS()` is typically used within composed nonlinear solver implementations, so most users would not generally call this routine themselves. .seealso: [](ch_snes), `SNESNGSFn`, `SNESSetNGS()`, `SNESComputeFunction()`, `SNESNGS` @*/ PetscErrorCode SNESComputeNGS(SNES snes, Vec b, Vec x) { DM dm; DMSNES sdm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(x, VEC_CLASSID, 3); if (b) PetscValidHeaderSpecific(b, VEC_CLASSID, 2); PetscCheckSameComm(snes, 1, x, 3); if (b) PetscCheckSameComm(snes, 1, b, 2); if (b) PetscCall(VecValidValues_Internal(b, 2, PETSC_TRUE)); PetscCall(PetscLogEventBegin(SNES_NGSEval, snes, x, b, 0)); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); PetscCheck(sdm->ops->computegs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetNGS() before SNESComputeNGS(), likely called from SNESSolve()."); if (b) PetscCall(VecLockReadPush(b)); PetscCallBack("SNES callback NGS", (*sdm->ops->computegs)(snes, x, b, sdm->gsctx)); if (b) PetscCall(VecLockReadPop(b)); PetscCall(PetscLogEventEnd(SNES_NGSEval, snes, x, b, 0)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESComputeFunction_FD(SNES snes, Vec Xin, Vec G) { Vec X; PetscScalar *g; PetscReal f, f2; PetscInt low, high, N, i; PetscBool flg; PetscReal h = .5 * PETSC_SQRT_MACHINE_EPSILON; PetscFunctionBegin; PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fd_delta", &h, &flg)); PetscCall(VecDuplicate(Xin, &X)); PetscCall(VecCopy(Xin, X)); PetscCall(VecGetSize(X, &N)); PetscCall(VecGetOwnershipRange(X, &low, &high)); PetscCall(VecSetOption(X, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE)); PetscCall(VecGetArray(G, &g)); for (i = 0; i < N; i++) { PetscCall(VecSetValue(X, i, -h, ADD_VALUES)); PetscCall(VecAssemblyBegin(X)); PetscCall(VecAssemblyEnd(X)); PetscCall(SNESComputeObjective(snes, X, &f)); PetscCall(VecSetValue(X, i, 2.0 * h, ADD_VALUES)); PetscCall(VecAssemblyBegin(X)); PetscCall(VecAssemblyEnd(X)); PetscCall(SNESComputeObjective(snes, X, &f2)); PetscCall(VecSetValue(X, i, -h, ADD_VALUES)); PetscCall(VecAssemblyBegin(X)); PetscCall(VecAssemblyEnd(X)); if (i >= low && i < high) g[i - low] = (f2 - f) / (2.0 * h); } PetscCall(VecRestoreArray(G, &g)); PetscCall(VecDestroy(&X)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESTestFunction - Computes the difference between the computed and finite-difference functions Collective Input Parameter: . snes - the `SNES` context Options Database Keys: + -snes_test_function - compare the user provided function with one compute via finite differences to check for errors. - -snes_test_function_view - display the user provided function, the finite difference function and the difference Level: developer .seealso: [](ch_snes), `SNESTestJacobian()`, `SNESSetFunction()`, `SNESComputeFunction()` @*/ PetscErrorCode SNESTestFunction(SNES snes) { Vec x, g1, g2, g3; PetscBool complete_print = PETSC_FALSE; PetscReal hcnorm, fdnorm, hcmax, fdmax, diffmax, diffnorm; PetscScalar dot; MPI_Comm comm; PetscViewer viewer, mviewer; PetscViewerFormat format; PetscInt tabs; static PetscBool directionsprinted = PETSC_FALSE; SNESObjectiveFn *objective; PetscFunctionBegin; PetscCall(SNESGetObjective(snes, &objective, NULL)); if (!objective) PetscFunctionReturn(PETSC_SUCCESS); PetscObjectOptionsBegin((PetscObject)snes); PetscCall(PetscOptionsViewer("-snes_test_function_view", "View difference between hand-coded and finite difference function element entries", "None", &mviewer, &format, &complete_print)); PetscOptionsEnd(); PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); PetscCall(PetscViewerASCIIGetStdout(comm, &viewer)); PetscCall(PetscViewerASCIIGetTab(viewer, &tabs)); PetscCall(PetscViewerASCIISetTab(viewer, ((PetscObject)snes)->tablevel)); PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Function -------------\n")); if (!complete_print && !directionsprinted) { PetscCall(PetscViewerASCIIPrintf(viewer, " Run with -snes_test_function_view and optionally -snes_test_function to show difference\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " of hand-coded and finite difference function entries greater than .\n")); } if (!directionsprinted) { PetscCall(PetscViewerASCIIPrintf(viewer, " Testing hand-coded Function, if (for double precision runs) ||F - Ffd||/||F|| is\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " O(1.e-8), the hand-coded Function is probably correct.\n")); directionsprinted = PETSC_TRUE; } if (complete_print) PetscCall(PetscViewerPushFormat(mviewer, format)); PetscCall(SNESGetSolution(snes, &x)); PetscCall(VecDuplicate(x, &g1)); PetscCall(VecDuplicate(x, &g2)); PetscCall(VecDuplicate(x, &g3)); PetscCall(SNESComputeFunction(snes, x, g1)); PetscCall(SNESComputeFunction_FD(snes, x, g2)); PetscCall(VecNorm(g2, NORM_2, &fdnorm)); PetscCall(VecNorm(g1, NORM_2, &hcnorm)); PetscCall(VecNorm(g2, NORM_INFINITY, &fdmax)); PetscCall(VecNorm(g1, NORM_INFINITY, &hcmax)); PetscCall(VecDot(g1, g2, &dot)); PetscCall(VecCopy(g1, g3)); PetscCall(VecAXPY(g3, -1.0, g2)); PetscCall(VecNorm(g3, NORM_2, &diffnorm)); PetscCall(VecNorm(g3, NORM_INFINITY, &diffmax)); PetscCall(PetscViewerASCIIPrintf(viewer, " ||Ffd|| %g, ||F|| = %g, angle cosine = (Ffd'F)/||Ffd||||F|| = %g\n", (double)fdnorm, (double)hcnorm, (double)(PetscRealPart(dot) / (fdnorm * hcnorm)))); PetscCall(PetscViewerASCIIPrintf(viewer, " 2-norm ||F - Ffd||/||F|| = %g, ||F - Ffd|| = %g\n", (double)(diffnorm / PetscMax(hcnorm, fdnorm)), (double)diffnorm)); PetscCall(PetscViewerASCIIPrintf(viewer, " max-norm ||F - Ffd||/||F|| = %g, ||F - Ffd|| = %g\n", (double)(diffmax / PetscMax(hcmax, fdmax)), (double)diffmax)); if (complete_print) { PetscCall(PetscViewerASCIIPrintf(viewer, " Hand-coded function ----------\n")); PetscCall(VecView(g1, mviewer)); PetscCall(PetscViewerASCIIPrintf(viewer, " Finite difference function ----------\n")); PetscCall(VecView(g2, mviewer)); PetscCall(PetscViewerASCIIPrintf(viewer, " Hand-coded minus finite-difference function ----------\n")); PetscCall(VecView(g3, mviewer)); } PetscCall(VecDestroy(&g1)); PetscCall(VecDestroy(&g2)); PetscCall(VecDestroy(&g3)); if (complete_print) { PetscCall(PetscViewerPopFormat(mviewer)); PetscCall(PetscViewerDestroy(&mviewer)); } PetscCall(PetscViewerASCIISetTab(viewer, tabs)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESTestJacobian - Computes the difference between the computed and finite-difference Jacobians Collective Input Parameter: . snes - the `SNES` context Output Parameters: + Jnorm - the Frobenius norm of the computed Jacobian, or `NULL` - diffNorm - the Frobenius norm of the difference of the computed and finite-difference Jacobians, or `NULL` Options Database Keys: + -snes_test_jacobian - compare the user provided Jacobian with one compute via finite differences to check for errors. If a threshold is given, display only those entries whose difference is greater than the threshold. - -snes_test_jacobian_view - display the user provided Jacobian, the finite difference Jacobian and the difference Level: developer Note: Directions and norms are printed to stdout if `diffNorm` is `NULL`. .seealso: [](ch_snes), `SNESTestFunction()`, `SNESSetJacobian()`, `SNESComputeJacobian()` @*/ PetscErrorCode SNESTestJacobian(SNES snes, PetscReal *Jnorm, PetscReal *diffNorm) { Mat A, B, C, D, jacobian; Vec x = snes->vec_sol, f; PetscReal nrm, gnorm; PetscReal threshold = 1.e-5; MatType mattype; PetscInt m, n, M, N; void *functx; PetscBool complete_print = PETSC_FALSE, threshold_print = PETSC_FALSE, flg, istranspose; PetscBool silent = diffNorm != PETSC_NULLPTR ? PETSC_TRUE : PETSC_FALSE; PetscViewer viewer, mviewer; MPI_Comm comm; PetscInt tabs; static PetscBool directionsprinted = PETSC_FALSE; PetscViewerFormat format; PetscFunctionBegin; PetscObjectOptionsBegin((PetscObject)snes); PetscCall(PetscOptionsReal("-snes_test_jacobian", "Threshold for element difference between hand-coded and finite difference being meaningful", "None", threshold, &threshold, NULL)); PetscCall(PetscOptionsDeprecated("-snes_test_jacobian_display", "-snes_test_jacobian_view", "3.13", NULL)); PetscCall(PetscOptionsViewer("-snes_test_jacobian_view", "View difference between hand-coded and finite difference Jacobians element entries", "None", &mviewer, &format, &complete_print)); PetscCall(PetscOptionsDeprecated("-snes_test_jacobian_display_threshold", "-snes_test_jacobian", "3.13", "-snes_test_jacobian accepts an optional threshold (since v3.10)")); PetscCall(PetscOptionsReal("-snes_test_jacobian_display_threshold", "Display difference between hand-coded and finite difference Jacobians which exceed input threshold", "None", threshold, &threshold, &threshold_print)); PetscOptionsEnd(); PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); PetscCall(PetscViewerASCIIGetStdout(comm, &viewer)); PetscCall(PetscViewerASCIIGetTab(viewer, &tabs)); PetscCall(PetscViewerASCIISetTab(viewer, ((PetscObject)snes)->tablevel)); if (!silent) PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Jacobian -------------\n")); if (!complete_print && !silent && !directionsprinted) { PetscCall(PetscViewerASCIIPrintf(viewer, " Run with -snes_test_jacobian_view and optionally -snes_test_jacobian to show difference\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " of hand-coded and finite difference Jacobian entries greater than .\n")); } if (!directionsprinted && !silent) { PetscCall(PetscViewerASCIIPrintf(viewer, " Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " O(1.e-8), the hand-coded Jacobian is probably correct.\n")); directionsprinted = PETSC_TRUE; } if (complete_print) PetscCall(PetscViewerPushFormat(mviewer, format)); PetscCall(PetscObjectTypeCompare((PetscObject)snes->jacobian, MATMFFD, &flg)); if (!flg) jacobian = snes->jacobian; else jacobian = snes->jacobian_pre; if (!x) PetscCall(MatCreateVecs(jacobian, &x, NULL)); else PetscCall(PetscObjectReference((PetscObject)x)); PetscCall(VecDuplicate(x, &f)); /* evaluate the function at this point because SNESComputeJacobianDefault() assumes that the function has been evaluated and put into snes->vec_func */ PetscCall(SNESComputeFunction(snes, x, f)); PetscCall(VecDestroy(&f)); PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESKSPTRANSPOSEONLY, &istranspose)); while (jacobian) { Mat JT = NULL, Jsave = NULL; if (istranspose) { PetscCall(MatCreateTranspose(jacobian, &JT)); Jsave = jacobian; jacobian = JT; } PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)jacobian, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, "")); if (flg) { A = jacobian; PetscCall(PetscObjectReference((PetscObject)A)); } else { PetscCall(MatComputeOperator(jacobian, MATAIJ, &A)); } PetscCall(MatGetType(A, &mattype)); PetscCall(MatGetSize(A, &M, &N)); PetscCall(MatGetLocalSize(A, &m, &n)); PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); PetscCall(MatSetType(B, mattype)); PetscCall(MatSetSizes(B, m, n, M, N)); PetscCall(MatSetBlockSizesFromMats(B, A, A)); PetscCall(MatSetUp(B)); PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); PetscCall(SNESGetFunction(snes, NULL, NULL, &functx)); PetscCall(SNESComputeJacobianDefault(snes, x, B, B, functx)); PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &D)); PetscCall(MatAYPX(D, -1.0, A, DIFFERENT_NONZERO_PATTERN)); PetscCall(MatNorm(D, NORM_FROBENIUS, &nrm)); PetscCall(MatNorm(A, NORM_FROBENIUS, &gnorm)); PetscCall(MatDestroy(&D)); if (!gnorm) gnorm = 1; /* just in case */ if (!silent) PetscCall(PetscViewerASCIIPrintf(viewer, " ||J - Jfd||_F/||J||_F = %g, ||J - Jfd||_F = %g\n", (double)(nrm / gnorm), (double)nrm)); if (complete_print) { PetscCall(PetscViewerASCIIPrintf(viewer, " Hand-coded Jacobian ----------\n")); PetscCall(MatView(A, mviewer)); PetscCall(PetscViewerASCIIPrintf(viewer, " Finite difference Jacobian ----------\n")); PetscCall(MatView(B, mviewer)); } if (threshold_print || complete_print) { PetscInt Istart, Iend, *ccols, bncols, cncols, j, row; PetscScalar *cvals; const PetscInt *bcols; const PetscScalar *bvals; PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); PetscCall(MatSetType(C, mattype)); PetscCall(MatSetSizes(C, m, n, M, N)); PetscCall(MatSetBlockSizesFromMats(C, A, A)); PetscCall(MatSetUp(C)); PetscCall(MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); PetscCall(MatAYPX(B, -1.0, A, DIFFERENT_NONZERO_PATTERN)); PetscCall(MatGetOwnershipRange(B, &Istart, &Iend)); for (row = Istart; row < Iend; row++) { PetscCall(MatGetRow(B, row, &bncols, &bcols, &bvals)); PetscCall(PetscMalloc2(bncols, &ccols, bncols, &cvals)); for (j = 0, cncols = 0; j < bncols; j++) { if (PetscAbsScalar(bvals[j]) > threshold) { ccols[cncols] = bcols[j]; cvals[cncols] = bvals[j]; cncols += 1; } } if (cncols) PetscCall(MatSetValues(C, 1, &row, cncols, ccols, cvals, INSERT_VALUES)); PetscCall(MatRestoreRow(B, row, &bncols, &bcols, &bvals)); PetscCall(PetscFree2(ccols, cvals)); } PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); PetscCall(PetscViewerASCIIPrintf(viewer, " Hand-coded minus finite-difference Jacobian with tolerance %g ----------\n", (double)threshold)); PetscCall(MatView(C, complete_print ? mviewer : viewer)); PetscCall(MatDestroy(&C)); } PetscCall(MatDestroy(&A)); PetscCall(MatDestroy(&B)); PetscCall(MatDestroy(&JT)); if (Jsave) jacobian = Jsave; if (jacobian != snes->jacobian_pre) { jacobian = snes->jacobian_pre; if (!silent) PetscCall(PetscViewerASCIIPrintf(viewer, " ---------- Testing Jacobian for preconditioner -------------\n")); } else jacobian = NULL; } PetscCall(VecDestroy(&x)); if (complete_print) PetscCall(PetscViewerPopFormat(mviewer)); if (mviewer) PetscCall(PetscViewerDestroy(&mviewer)); PetscCall(PetscViewerASCIISetTab(viewer, tabs)); if (Jnorm) *Jnorm = gnorm; if (diffNorm) *diffNorm = nrm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESComputeJacobian - Computes the Jacobian matrix that has been set with `SNESSetJacobian()`. Collective Input Parameters: + snes - the `SNES` context - X - input vector Output Parameters: + A - Jacobian matrix - B - optional matrix for building the preconditioner, usually the same as `A` Options Database Keys: + -snes_lag_preconditioner - how often to rebuild preconditioner . -snes_lag_jacobian - how often to rebuild Jacobian . -snes_test_jacobian - compare the user provided Jacobian with one compute via finite differences to check for errors. If a threshold is given, display only those entries whose difference is greater than the threshold. . -snes_test_jacobian_view - display the user provided Jacobian, the finite difference Jacobian and the difference between them to help users detect the location of errors in the user provided Jacobian . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result . -snes_compare_operator - Make the comparison options above use the operator instead of the matrix used to construct the preconditioner . -snes_compare_coloring - Compute the finite difference Jacobian using coloring and display norms of difference . -snes_compare_coloring_display - Compute the finite difference Jacobian using coloring and display verbose differences . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by `-snes_compare_coloring_threshold` . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by `-snes_compare_coloring_threshold` . -snes_compare_coloring_draw - Compute the finite difference Jacobian using coloring and draw differences - -snes_compare_coloring_draw_contour - Compute the finite difference Jacobian using coloring and show contours of matrices and differences Level: developer Note: Most users should not need to explicitly call this routine, as it is used internally within the nonlinear solvers. Developer Note: This has duplicative ways of checking the accuracy of the user provided Jacobian (see the options above). This is for historical reasons, the routine `SNESTestJacobian()` use to used with the `SNESType` of test that has been removed. .seealso: [](ch_snes), `SNESSetJacobian()`, `KSPSetOperators()`, `MatStructure`, `SNESSetLagPreconditioner()`, `SNESSetLagJacobian()` @*/ PetscErrorCode SNESComputeJacobian(SNES snes, Vec X, Mat A, Mat B) { PetscBool flag; DM dm; DMSNES sdm; KSP ksp; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(X, VEC_CLASSID, 2); PetscCheckSameComm(snes, 1, X, 2); PetscCall(VecValidValues_Internal(X, 2, PETSC_TRUE)); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix-free */ if (snes->lagjacobian == -2) { snes->lagjacobian = -1; PetscCall(PetscInfo(snes, "Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n")); } else if (snes->lagjacobian == -1) { PetscCall(PetscInfo(snes, "Reusing Jacobian/preconditioner because lag is -1\n")); PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flag)); if (flag) { PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); } PetscFunctionReturn(PETSC_SUCCESS); } else if (snes->lagjacobian > 1 && (snes->iter + snes->jac_iter) % snes->lagjacobian) { PetscCall(PetscInfo(snes, "Reusing Jacobian/preconditioner because lag is %" PetscInt_FMT " and SNES iteration is %" PetscInt_FMT "\n", snes->lagjacobian, snes->iter)); PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flag)); if (flag) { PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); } PetscFunctionReturn(PETSC_SUCCESS); } if (snes->npc && snes->npcside == PC_LEFT) { PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscLogEventBegin(SNES_JacobianEval, snes, X, A, B)); PetscCall(VecLockReadPush(X)); { void *ctx; SNESJacobianFn *J; PetscCall(DMSNESGetJacobian(dm, &J, &ctx)); PetscCallBack("SNES callback Jacobian", (*J)(snes, X, A, B, ctx)); } PetscCall(VecLockReadPop(X)); PetscCall(PetscLogEventEnd(SNES_JacobianEval, snes, X, A, B)); /* attach latest linearization point to the matrix used to construct the preconditioner */ PetscCall(PetscObjectCompose((PetscObject)B, "__SNES_latest_X", (PetscObject)X)); /* the next line ensures that snes->ksp exists */ PetscCall(SNESGetKSP(snes, &ksp)); if (snes->lagpreconditioner == -2) { PetscCall(PetscInfo(snes, "Rebuilding preconditioner exactly once since lag is -2\n")); PetscCall(KSPSetReusePreconditioner(snes->ksp, PETSC_FALSE)); snes->lagpreconditioner = -1; } else if (snes->lagpreconditioner == -1) { PetscCall(PetscInfo(snes, "Reusing preconditioner because lag is -1\n")); PetscCall(KSPSetReusePreconditioner(snes->ksp, PETSC_TRUE)); } else if (snes->lagpreconditioner > 1 && (snes->iter + snes->pre_iter) % snes->lagpreconditioner) { PetscCall(PetscInfo(snes, "Reusing preconditioner because lag is %" PetscInt_FMT " and SNES iteration is %" PetscInt_FMT "\n", snes->lagpreconditioner, snes->iter)); PetscCall(KSPSetReusePreconditioner(snes->ksp, PETSC_TRUE)); } else { PetscCall(PetscInfo(snes, "Rebuilding preconditioner\n")); PetscCall(KSPSetReusePreconditioner(snes->ksp, PETSC_FALSE)); } /* monkey business to allow testing Jacobians in multilevel solvers. This is needed because the SNESTestXXX interface does not accept vectors and matrices */ { Vec xsave = snes->vec_sol; Mat jacobiansave = snes->jacobian; Mat jacobian_presave = snes->jacobian_pre; snes->vec_sol = X; snes->jacobian = A; snes->jacobian_pre = B; if (snes->testFunc) PetscCall(SNESTestFunction(snes)); if (snes->testJac) PetscCall(SNESTestJacobian(snes, NULL, NULL)); snes->vec_sol = xsave; snes->jacobian = jacobiansave; snes->jacobian_pre = jacobian_presave; } { PetscBool flag = PETSC_FALSE, flag_draw = PETSC_FALSE, flag_contour = PETSC_FALSE, flag_operator = PETSC_FALSE; PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_explicit", NULL, NULL, &flag)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_explicit_draw", NULL, NULL, &flag_draw)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_explicit_draw_contour", NULL, NULL, &flag_contour)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_operator", NULL, NULL, &flag_operator)); if (flag || flag_draw || flag_contour) { Mat Bexp_mine = NULL, Bexp, FDexp; PetscViewer vdraw, vstdout; PetscBool flg; if (flag_operator) { PetscCall(MatComputeOperator(A, MATAIJ, &Bexp_mine)); Bexp = Bexp_mine; } else { /* See if the matrix used to construct the preconditioner can be viewed and added directly */ PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)B, &flg, MATSEQAIJ, MATMPIAIJ, MATSEQDENSE, MATMPIDENSE, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, "")); if (flg) Bexp = B; else { /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ PetscCall(MatComputeOperator(B, MATAIJ, &Bexp_mine)); Bexp = Bexp_mine; } } PetscCall(MatConvert(Bexp, MATSAME, MAT_INITIAL_MATRIX, &FDexp)); PetscCall(SNESComputeJacobianDefault(snes, X, FDexp, FDexp, NULL)); PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &vstdout)); if (flag_draw || flag_contour) { PetscCall(PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes), NULL, "Explicit Jacobians", PETSC_DECIDE, PETSC_DECIDE, 300, 300, &vdraw)); if (flag_contour) PetscCall(PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR)); } else vdraw = NULL; PetscCall(PetscViewerASCIIPrintf(vstdout, "Explicit %s\n", flag_operator ? "Jacobian" : "preconditioning Jacobian")); if (flag) PetscCall(MatView(Bexp, vstdout)); if (vdraw) PetscCall(MatView(Bexp, vdraw)); PetscCall(PetscViewerASCIIPrintf(vstdout, "Finite difference Jacobian\n")); if (flag) PetscCall(MatView(FDexp, vstdout)); if (vdraw) PetscCall(MatView(FDexp, vdraw)); PetscCall(MatAYPX(FDexp, -1.0, Bexp, SAME_NONZERO_PATTERN)); PetscCall(PetscViewerASCIIPrintf(vstdout, "User-provided matrix minus finite difference Jacobian\n")); if (flag) PetscCall(MatView(FDexp, vstdout)); if (vdraw) { /* Always use contour for the difference */ PetscCall(PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR)); PetscCall(MatView(FDexp, vdraw)); PetscCall(PetscViewerPopFormat(vdraw)); } if (flag_contour) PetscCall(PetscViewerPopFormat(vdraw)); PetscCall(PetscViewerDestroy(&vdraw)); PetscCall(MatDestroy(&Bexp_mine)); PetscCall(MatDestroy(&FDexp)); } } { PetscBool flag = PETSC_FALSE, flag_display = PETSC_FALSE, flag_draw = PETSC_FALSE, flag_contour = PETSC_FALSE, flag_threshold = PETSC_FALSE; PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON, threshold_rtol = 10 * PETSC_SQRT_MACHINE_EPSILON; PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring", NULL, NULL, &flag)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_display", NULL, NULL, &flag_display)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_draw", NULL, NULL, &flag_draw)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_draw_contour", NULL, NULL, &flag_contour)); PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_threshold", NULL, NULL, &flag_threshold)); if (flag_threshold) { PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_threshold_rtol", &threshold_rtol, NULL)); PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_compare_coloring_threshold_atol", &threshold_atol, NULL)); } if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { Mat Bfd; PetscViewer vdraw, vstdout; MatColoring coloring; ISColoring iscoloring; MatFDColoring matfdcoloring; SNESFunctionFn *func; void *funcctx; PetscReal norm1, norm2, normmax; PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &Bfd)); PetscCall(MatColoringCreate(Bfd, &coloring)); PetscCall(MatColoringSetType(coloring, MATCOLORINGSL)); PetscCall(MatColoringSetFromOptions(coloring)); PetscCall(MatColoringApply(coloring, &iscoloring)); PetscCall(MatColoringDestroy(&coloring)); PetscCall(MatFDColoringCreate(Bfd, iscoloring, &matfdcoloring)); PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); PetscCall(MatFDColoringSetUp(Bfd, iscoloring, matfdcoloring)); PetscCall(ISColoringDestroy(&iscoloring)); /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ PetscCall(SNESGetFunction(snes, NULL, &func, &funcctx)); PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)func, funcctx)); PetscCall(PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring, ((PetscObject)snes)->prefix)); PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring, "coloring_")); PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); PetscCall(MatFDColoringApply(Bfd, matfdcoloring, X, snes)); PetscCall(MatFDColoringDestroy(&matfdcoloring)); PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &vstdout)); if (flag_draw || flag_contour) { PetscCall(PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes), NULL, "Colored Jacobians", PETSC_DECIDE, PETSC_DECIDE, 300, 300, &vdraw)); if (flag_contour) PetscCall(PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR)); } else vdraw = NULL; PetscCall(PetscViewerASCIIPrintf(vstdout, "Explicit preconditioning Jacobian\n")); if (flag_display) PetscCall(MatView(B, vstdout)); if (vdraw) PetscCall(MatView(B, vdraw)); PetscCall(PetscViewerASCIIPrintf(vstdout, "Colored Finite difference Jacobian\n")); if (flag_display) PetscCall(MatView(Bfd, vstdout)); if (vdraw) PetscCall(MatView(Bfd, vdraw)); PetscCall(MatAYPX(Bfd, -1.0, B, SAME_NONZERO_PATTERN)); PetscCall(MatNorm(Bfd, NORM_1, &norm1)); PetscCall(MatNorm(Bfd, NORM_FROBENIUS, &norm2)); PetscCall(MatNorm(Bfd, NORM_MAX, &normmax)); PetscCall(PetscViewerASCIIPrintf(vstdout, "User-provided matrix minus finite difference Jacobian, norm1=%g normFrob=%g normmax=%g\n", (double)norm1, (double)norm2, (double)normmax)); if (flag_display) PetscCall(MatView(Bfd, vstdout)); if (vdraw) { /* Always use contour for the difference */ PetscCall(PetscViewerPushFormat(vdraw, PETSC_VIEWER_DRAW_CONTOUR)); PetscCall(MatView(Bfd, vdraw)); PetscCall(PetscViewerPopFormat(vdraw)); } if (flag_contour) PetscCall(PetscViewerPopFormat(vdraw)); if (flag_threshold) { PetscInt bs, rstart, rend, i; PetscCall(MatGetBlockSize(B, &bs)); PetscCall(MatGetOwnershipRange(B, &rstart, &rend)); for (i = rstart; i < rend; i++) { const PetscScalar *ba, *ca; const PetscInt *bj, *cj; PetscInt bn, cn, j, maxentrycol = -1, maxdiffcol = -1, maxrdiffcol = -1; PetscReal maxentry = 0, maxdiff = 0, maxrdiff = 0; PetscCall(MatGetRow(B, i, &bn, &bj, &ba)); PetscCall(MatGetRow(Bfd, i, &cn, &cj, &ca)); PetscCheck(bn == cn, ((PetscObject)A)->comm, PETSC_ERR_PLIB, "Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); for (j = 0; j < bn; j++) { PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol * PetscAbsScalar(ba[j])); if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { maxentrycol = bj[j]; maxentry = PetscRealPart(ba[j]); } if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { maxdiffcol = bj[j]; maxdiff = PetscRealPart(ca[j]); } if (rdiff > maxrdiff) { maxrdiffcol = bj[j]; maxrdiff = rdiff; } } if (maxrdiff > 1) { PetscCall(PetscViewerASCIIPrintf(vstdout, "row %" PetscInt_FMT " (maxentry=%g at %" PetscInt_FMT ", maxdiff=%g at %" PetscInt_FMT ", maxrdiff=%g at %" PetscInt_FMT "):", i, (double)maxentry, maxentrycol, (double)maxdiff, maxdiffcol, (double)maxrdiff, maxrdiffcol)); for (j = 0; j < bn; j++) { PetscReal rdiff; rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol * PetscAbsScalar(ba[j])); if (rdiff > 1) PetscCall(PetscViewerASCIIPrintf(vstdout, " (%" PetscInt_FMT ",%g:%g)", bj[j], (double)PetscRealPart(ba[j]), (double)PetscRealPart(ca[j]))); } PetscCall(PetscViewerASCIIPrintf(vstdout, "\n")); } PetscCall(MatRestoreRow(B, i, &bn, &bj, &ba)); PetscCall(MatRestoreRow(Bfd, i, &cn, &cj, &ca)); } } PetscCall(PetscViewerDestroy(&vdraw)); PetscCall(MatDestroy(&Bfd)); } } PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetJacobian - Sets the function to compute Jacobian as well as the location to store the matrix. Logically Collective Input Parameters: + snes - the `SNES` context . Amat - the matrix that defines the (approximate) Jacobian . Pmat - the matrix to be used in constructing the preconditioner, usually the same as `Amat`. . J - Jacobian evaluation routine (if `NULL` then `SNES` retains any previously set value), see `SNESJacobianFn` for details - ctx - [optional] user-defined context for private data for the Jacobian evaluation routine (may be `NULL`) (if `NULL` then `SNES` retains any previously set value) Level: beginner Notes: If the `Amat` matrix and `Pmat` matrix are different you must call `MatAssemblyBegin()`/`MatAssemblyEnd()` on each matrix. If you know the operator `Amat` has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process. If using `SNESComputeJacobianDefaultColor()` to assemble a Jacobian, the `ctx` argument must be a `MatFDColoring`. Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common example is to use the "Picard linearization" which only differentiates through the highest order parts of each term using `SNESSetPicard()` .seealso: [](ch_snes), `SNES`, `KSPSetOperators()`, `SNESSetFunction()`, `MatMFFDComputeJacobian()`, `SNESComputeJacobianDefaultColor()`, `MatStructure`, `SNESSetPicard()`, `SNESJacobianFn`, `SNESFunctionFn` @*/ PetscErrorCode SNESSetJacobian(SNES snes, Mat Amat, Mat Pmat, SNESJacobianFn *J, void *ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (Amat) PetscValidHeaderSpecific(Amat, MAT_CLASSID, 2); if (Pmat) PetscValidHeaderSpecific(Pmat, MAT_CLASSID, 3); if (Amat) PetscCheckSameComm(snes, 1, Amat, 2); if (Pmat) PetscCheckSameComm(snes, 1, Pmat, 3); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESSetJacobian(dm, J, ctx)); if (Amat) { PetscCall(PetscObjectReference((PetscObject)Amat)); PetscCall(MatDestroy(&snes->jacobian)); snes->jacobian = Amat; } if (Pmat) { PetscCall(PetscObjectReference((PetscObject)Pmat)); PetscCall(MatDestroy(&snes->jacobian_pre)); snes->jacobian_pre = Pmat; } PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESGetJacobian - Returns the Jacobian matrix and optionally the user provided context for evaluating the Jacobian. Not Collective, but `Mat` object will be parallel if `SNES` is Input Parameter: . snes - the nonlinear solver context Output Parameters: + Amat - location to stash (approximate) Jacobian matrix (or `NULL`) . Pmat - location to stash matrix used to compute the preconditioner (or `NULL`) . J - location to put Jacobian function (or `NULL`), for calling sequence see `SNESJacobianFn` - ctx - location to stash Jacobian ctx (or `NULL`) Level: advanced .seealso: [](ch_snes), `SNES`, `Mat`, `SNESSetJacobian()`, `SNESComputeJacobian()`, `SNESJacobianFn`, `SNESGetFunction()` @*/ PetscErrorCode SNESGetJacobian(SNES snes, Mat *Amat, Mat *Pmat, SNESJacobianFn **J, void **ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (Amat) *Amat = snes->jacobian; if (Pmat) *Pmat = snes->jacobian_pre; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESGetJacobian(dm, J, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode SNESSetDefaultComputeJacobian(SNES snes) { DM dm; DMSNES sdm; PetscFunctionBegin; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); if (!sdm->ops->computejacobian && snes->jacobian_pre) { DM dm; PetscBool isdense, ismf; PetscCall(SNESGetDM(snes, &dm)); PetscCall(PetscObjectTypeCompareAny((PetscObject)snes->jacobian_pre, &isdense, MATSEQDENSE, MATMPIDENSE, MATDENSE, NULL)); PetscCall(PetscObjectTypeCompareAny((PetscObject)snes->jacobian_pre, &ismf, MATMFFD, MATSHELL, NULL)); if (isdense) { PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobianDefault, NULL)); } else if (!ismf) { PetscCall(DMSNESSetJacobian(dm, SNESComputeJacobianDefaultColor, NULL)); } } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetUp - Sets up the internal data structures for the later use of a nonlinear solver `SNESSolve()`. Collective Input Parameter: . snes - the `SNES` context Level: advanced Note: For basic use of the `SNES` solvers the user does not need to explicitly call `SNESSetUp()`, since these actions will automatically occur during the call to `SNESSolve()`. However, if one wishes to control this phase separately, `SNESSetUp()` should be called after `SNESCreate()` and optional routines of the form SNESSetXXX(), but before `SNESSolve()`. .seealso: [](ch_snes), `SNES`, `SNESCreate()`, `SNESSolve()`, `SNESDestroy()`, `SNESSetFromOptions()` @*/ PetscErrorCode SNESSetUp(SNES snes) { DM dm; DMSNES sdm; SNESLineSearch linesearch, pclinesearch; void *lsprectx, *lspostctx; PetscBool mf_operator, mf; Vec f, fpc; void *funcctx; void *jacctx, *appctx; Mat j, jpre; PetscErrorCode (*precheck)(SNESLineSearch, Vec, Vec, PetscBool *, void *); PetscErrorCode (*postcheck)(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *); SNESFunctionFn *func; SNESJacobianFn *jac; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (snes->setupcalled) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(PetscLogEventBegin(SNES_SetUp, snes, 0, 0, 0)); if (!((PetscObject)snes)->type_name) PetscCall(SNESSetType(snes, SNESNEWTONLS)); PetscCall(SNESGetFunction(snes, &snes->vec_func, NULL, NULL)); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetDMSNES(dm, &sdm)); PetscCall(SNESSetDefaultComputeJacobian(snes)); if (!snes->vec_func) PetscCall(DMCreateGlobalVector(dm, &snes->vec_func)); if (!snes->ksp) PetscCall(SNESGetKSP(snes, &snes->ksp)); if (snes->linesearch) { PetscCall(SNESGetLineSearch(snes, &snes->linesearch)); PetscCall(SNESLineSearchSetFunction(snes->linesearch, SNESComputeFunction)); } PetscCall(SNESGetUseMatrixFree(snes, &mf_operator, &mf)); if (snes->npc && snes->npcside == PC_LEFT) { snes->mf = PETSC_TRUE; snes->mf_operator = PETSC_FALSE; } if (snes->npc) { /* copy the DM over */ PetscCall(SNESGetDM(snes, &dm)); PetscCall(SNESSetDM(snes->npc, dm)); PetscCall(SNESGetFunction(snes, &f, &func, &funcctx)); PetscCall(VecDuplicate(f, &fpc)); PetscCall(SNESSetFunction(snes->npc, fpc, func, funcctx)); PetscCall(SNESGetJacobian(snes, &j, &jpre, &jac, &jacctx)); PetscCall(SNESSetJacobian(snes->npc, j, jpre, jac, jacctx)); PetscCall(SNESGetApplicationContext(snes, &appctx)); PetscCall(SNESSetApplicationContext(snes->npc, appctx)); PetscCall(SNESSetUseMatrixFree(snes->npc, mf_operator, mf)); PetscCall(VecDestroy(&fpc)); /* copy the function pointers over */ PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)snes, (PetscObject)snes->npc)); /* default to 1 iteration */ PetscCall(SNESSetTolerances(snes->npc, 0.0, 0.0, 0.0, 1, snes->npc->max_funcs)); if (snes->npcside == PC_RIGHT) { PetscCall(SNESSetNormSchedule(snes->npc, SNES_NORM_FINAL_ONLY)); } else { PetscCall(SNESSetNormSchedule(snes->npc, SNES_NORM_NONE)); } PetscCall(SNESSetFromOptions(snes->npc)); /* copy the line search context over */ if (snes->linesearch && snes->npc->linesearch) { PetscCall(SNESGetLineSearch(snes, &linesearch)); PetscCall(SNESGetLineSearch(snes->npc, &pclinesearch)); PetscCall(SNESLineSearchGetPreCheck(linesearch, &precheck, &lsprectx)); PetscCall(SNESLineSearchGetPostCheck(linesearch, &postcheck, &lspostctx)); PetscCall(SNESLineSearchSetPreCheck(pclinesearch, precheck, lsprectx)); PetscCall(SNESLineSearchSetPostCheck(pclinesearch, postcheck, lspostctx)); PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)linesearch, (PetscObject)pclinesearch)); } } if (snes->mf) PetscCall(SNESSetUpMatrixFree_Private(snes, snes->mf_operator, snes->mf_version)); if (snes->ops->usercompute && !snes->ctx) PetscCallBack("SNES callback compute application context", (*snes->ops->usercompute)(snes, &snes->ctx)); snes->jac_iter = 0; snes->pre_iter = 0; PetscTryTypeMethod(snes, setup); PetscCall(SNESSetDefaultComputeJacobian(snes)); if (snes->npc && snes->npcside == PC_LEFT) { if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { if (snes->linesearch) { PetscCall(SNESGetLineSearch(snes, &linesearch)); PetscCall(SNESLineSearchSetFunction(linesearch, SNESComputeFunctionDefaultNPC)); } } } PetscCall(PetscLogEventEnd(SNES_SetUp, snes, 0, 0, 0)); snes->setupcalled = PETSC_TRUE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESReset - Resets a `SNES` context to the state it was in before `SNESSetUp()` was called and removes any allocated `Vec` and `Mat` from its data structures Collective Input Parameter: . snes - the nonlinear iterative solver context obtained from `SNESCreate()` Level: intermediate Notes: Any options set on the `SNES` object, including those set with `SNESSetFromOptions()` remain. Call this if you wish to reuse a `SNES` but with different size vectors Also calls the application context destroy routine set with `SNESSetComputeApplicationContext()` .seealso: [](ch_snes), `SNES`, `SNESDestroy()`, `SNESCreate()`, `SNESSetUp()`, `SNESSolve()` @*/ PetscErrorCode SNESReset(SNES snes) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (snes->ops->ctxdestroy && snes->ctx) { PetscCallBack("SNES callback destroy application context", (*snes->ops->ctxdestroy)(&snes->ctx)); snes->ctx = NULL; } if (snes->npc) PetscCall(SNESReset(snes->npc)); PetscTryTypeMethod(snes, reset); if (snes->ksp) PetscCall(KSPReset(snes->ksp)); if (snes->linesearch) PetscCall(SNESLineSearchReset(snes->linesearch)); PetscCall(VecDestroy(&snes->vec_rhs)); PetscCall(VecDestroy(&snes->vec_sol)); PetscCall(VecDestroy(&snes->vec_sol_update)); PetscCall(VecDestroy(&snes->vec_func)); PetscCall(MatDestroy(&snes->jacobian)); PetscCall(MatDestroy(&snes->jacobian_pre)); PetscCall(MatDestroy(&snes->picard)); PetscCall(VecDestroyVecs(snes->nwork, &snes->work)); PetscCall(VecDestroyVecs(snes->nvwork, &snes->vwork)); snes->alwayscomputesfinalresidual = PETSC_FALSE; snes->nwork = snes->nvwork = 0; snes->setupcalled = PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESConvergedReasonViewCancel - Clears all the reason view functions for a `SNES` object provided with `SNESConvergedReasonViewSet()` also removes the default viewer. Collective Input Parameter: . snes - the nonlinear iterative solver context obtained from `SNESCreate()` Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESCreate()`, `SNESDestroy()`, `SNESReset()`, `SNESConvergedReasonViewSet()` @*/ PetscErrorCode SNESConvergedReasonViewCancel(SNES snes) { PetscInt i; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); for (i = 0; i < snes->numberreasonviews; i++) { if (snes->reasonviewdestroy[i]) PetscCall((*snes->reasonviewdestroy[i])(&snes->reasonviewcontext[i])); } snes->numberreasonviews = 0; PetscCall(PetscViewerDestroy(&snes->convergedreasonviewer)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESDestroy - Destroys the nonlinear solver context that was created with `SNESCreate()`. Collective Input Parameter: . snes - the `SNES` context Level: beginner .seealso: [](ch_snes), `SNES`, `SNESCreate()`, `SNESSolve()` @*/ PetscErrorCode SNESDestroy(SNES *snes) { DM dm; PetscFunctionBegin; if (!*snes) PetscFunctionReturn(PETSC_SUCCESS); PetscValidHeaderSpecific(*snes, SNES_CLASSID, 1); if (--((PetscObject)*snes)->refct > 0) { *snes = NULL; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(SNESReset(*snes)); PetscCall(SNESDestroy(&(*snes)->npc)); /* if memory was published with SAWs then destroy it */ PetscCall(PetscObjectSAWsViewOff((PetscObject)*snes)); PetscTryTypeMethod(*snes, destroy); dm = (*snes)->dm; while (dm) { PetscCall(DMCoarsenHookRemove(dm, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, *snes)); PetscCall(DMGetCoarseDM(dm, &dm)); } PetscCall(DMDestroy(&(*snes)->dm)); PetscCall(KSPDestroy(&(*snes)->ksp)); PetscCall(SNESLineSearchDestroy(&(*snes)->linesearch)); PetscCall(PetscFree((*snes)->kspconvctx)); if ((*snes)->ops->convergeddestroy) PetscCall((*(*snes)->ops->convergeddestroy)(&(*snes)->cnvP)); if ((*snes)->conv_hist_alloc) PetscCall(PetscFree2((*snes)->conv_hist, (*snes)->conv_hist_its)); PetscCall(SNESMonitorCancel(*snes)); PetscCall(SNESConvergedReasonViewCancel(*snes)); PetscCall(PetscHeaderDestroy(snes)); PetscFunctionReturn(PETSC_SUCCESS); } /* ----------- Routines to set solver parameters ---------- */ /*@ SNESSetLagPreconditioner - Sets when the preconditioner is rebuilt in the nonlinear solve `SNESSolve()`. Logically Collective Input Parameters: + snes - the `SNES` context - lag - 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that Options Database Keys: + -snes_lag_jacobian_persists - sets the persistence through multiple `SNESSolve()` . -snes_lag_jacobian <-2,1,2,...> - sets the lag . -snes_lag_preconditioner_persists - sets the persistence through multiple `SNESSolve()` - -snes_lag_preconditioner <-2,1,2,...> - sets the lag Level: intermediate Notes: The default is 1 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 or `SNESSetLagPreconditionerPersists()` was called `SNESSetLagPreconditionerPersists()` allows using the same uniform lagging (for example every second linear solve) across multiple nonlinear solves. .seealso: [](ch_snes), `SNESGetLagPreconditioner()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESSetLagPreconditionerPersists()`, `SNESSetLagJacobianPersists()`, `SNES`, `SNESSolve()` @*/ PetscErrorCode SNESSetLagPreconditioner(SNES snes, PetscInt lag) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCheck(lag >= -2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lag must be -2, -1, 1 or greater"); PetscCheck(lag, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lag cannot be 0"); PetscValidLogicalCollectiveInt(snes, lag, 2); snes->lagpreconditioner = lag; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetGridSequence - sets the number of steps of grid sequencing that `SNES` will do Logically Collective Input Parameters: + snes - the `SNES` context - steps - the number of refinements to do, defaults to 0 Options Database Key: . -snes_grid_sequence - Use grid sequencing to generate initial guess Level: intermediate Notes: Once grid sequencing is turned on `SNESSolve()` will automatically perform the solve on each grid refinement. Use `SNESGetSolution()` to extract the fine grid solution after grid sequencing. .seealso: [](ch_snes), `SNES`, `SNESGetLagPreconditioner()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESGetGridSequence()`, `SNESSetDM()`, `SNESSolve()` @*/ PetscErrorCode SNESSetGridSequence(SNES snes, PetscInt steps) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveInt(snes, steps, 2); snes->gridsequence = steps; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetGridSequence - gets the number of steps of grid sequencing that `SNES` will do Logically Collective Input Parameter: . snes - the `SNES` context Output Parameter: . steps - the number of refinements to do, defaults to 0 Level: intermediate .seealso: [](ch_snes), `SNESGetLagPreconditioner()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESSetGridSequence()` @*/ PetscErrorCode SNESGetGridSequence(SNES snes, PetscInt *steps) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *steps = snes->gridsequence; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetLagPreconditioner - Return how often the preconditioner is rebuilt Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that Level: intermediate Notes: The default is 1 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 .seealso: [](ch_snes), `SNES`, `SNESSetLagPreconditioner()`, `SNESSetLagJacobianPersists()`, `SNESSetLagPreconditionerPersists()` @*/ PetscErrorCode SNESGetLagPreconditioner(SNES snes, PetscInt *lag) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *lag = snes->lagpreconditioner; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetLagJacobian - Set when the Jacobian is rebuilt in the nonlinear solve. See `SNESSetLagPreconditioner()` for determining how often the preconditioner is rebuilt. Logically Collective Input Parameters: + snes - the `SNES` context - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time the Jacobian is built etc. -2 means rebuild at next chance but then never again Options Database Keys: + -snes_lag_jacobian_persists - sets the persistence through multiple SNES solves . -snes_lag_jacobian <-2,1,2,...> - sets the lag . -snes_lag_preconditioner_persists - sets the persistence through multiple SNES solves - -snes_lag_preconditioner <-2,1,2,...> - sets the lag. Level: intermediate Notes: The default is 1 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 If -1 is used before the very first nonlinear solve the CODE WILL FAIL! because no Jacobian is used, use -2 to indicate you want it recomputed at the next Newton step but never again (unless it is reset to another value) .seealso: [](ch_snes), `SNES`, `SNESGetLagPreconditioner()`, `SNESSetLagPreconditioner()`, `SNESGetLagJacobianPersists()`, `SNESSetLagPreconditionerPersists()` @*/ PetscErrorCode SNESSetLagJacobian(SNES snes, PetscInt lag) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCheck(lag >= -2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lag must be -2, -1, 1 or greater"); PetscCheck(lag, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lag cannot be 0"); PetscValidLogicalCollectiveInt(snes, lag, 2); snes->lagjacobian = lag; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetLagJacobian - Get how often the Jacobian is rebuilt. See `SNESGetLagPreconditioner()` to determine when the preconditioner is rebuilt Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time the Jacobian is built etc. Level: intermediate Notes: The default is 1 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 or `SNESSetLagJacobianPersists()` was called. .seealso: [](ch_snes), `SNES`, `SNESSetLagJacobian()`, `SNESSetLagPreconditioner()`, `SNESGetLagPreconditioner()`, `SNESSetLagJacobianPersists()`, `SNESSetLagPreconditionerPersists()` @*/ PetscErrorCode SNESGetLagJacobian(SNES snes, PetscInt *lag) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *lag = snes->lagjacobian; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetLagJacobianPersists - Set whether or not the Jacobian lagging persists through multiple nonlinear solves Logically collective Input Parameters: + snes - the `SNES` context - flg - jacobian lagging persists if true Options Database Keys: + -snes_lag_jacobian_persists - sets the persistence through multiple SNES solves . -snes_lag_jacobian <-2,1,2,...> - sets the lag . -snes_lag_preconditioner_persists - sets the persistence through multiple SNES solves - -snes_lag_preconditioner <-2,1,2,...> - sets the lag Level: advanced Notes: Normally when `SNESSetLagJacobian()` is used, the Jacobian is always rebuilt at the beginning of each new nonlinear solve, this removes that behavior This is useful both for nonlinear preconditioning, where it's appropriate to have the Jacobian be stale by several solves, and for implicit time-stepping, where Jacobian lagging in the inner nonlinear solve over several timesteps may present huge efficiency gains. .seealso: [](ch_snes), `SNES`, `SNESSetLagPreconditionerPersists()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESGetNPC()` @*/ PetscErrorCode SNESSetLagJacobianPersists(SNES snes, PetscBool flg) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveBool(snes, flg, 2); snes->lagjac_persist = flg; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetLagPreconditionerPersists - Set whether or not the preconditioner lagging persists through multiple nonlinear solves Logically Collective Input Parameters: + snes - the `SNES` context - flg - preconditioner lagging persists if true Options Database Keys: + -snes_lag_jacobian_persists - sets the persistence through multiple SNES solves . -snes_lag_jacobian <-2,1,2,...> - sets the lag . -snes_lag_preconditioner_persists - sets the persistence through multiple SNES solves - -snes_lag_preconditioner <-2,1,2,...> - sets the lag Level: developer Notes: Normally when `SNESSetLagPreconditioner()` is used, the preconditioner is always rebuilt at the beginning of each new nonlinear solve, this removes that behavior This is useful both for nonlinear preconditioning, where it's appropriate to have the preconditioner be stale by several solves, and for implicit time-stepping, where preconditioner lagging in the inner nonlinear solve over several timesteps may present huge efficiency gains. .seealso: [](ch_snes), `SNES`, `SNESSetLagJacobianPersists()`, `SNESSetLagJacobian()`, `SNESGetLagJacobian()`, `SNESGetNPC()`, `SNESSetLagPreconditioner()` @*/ PetscErrorCode SNESSetLagPreconditionerPersists(SNES snes, PetscBool flg) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveBool(snes, flg, 2); snes->lagpre_persist = flg; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetForceIteration - force `SNESSolve()` to take at least one iteration regardless of the initial residual norm Logically Collective Input Parameters: + snes - the `SNES` context - force - `PETSC_TRUE` require at least one iteration Options Database Key: . -snes_force_iteration - Sets forcing an iteration Level: intermediate Note: This is used sometimes with `TS` to prevent `TS` from detecting a false steady state solution .seealso: [](ch_snes), `SNES`, `TS`, `SNESSetDivergenceTolerance()` @*/ PetscErrorCode SNESSetForceIteration(SNES snes, PetscBool force) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->forceiteration = force; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetForceIteration - Check whether or not `SNESSolve()` take at least one iteration regardless of the initial residual norm Logically Collective Input Parameter: . snes - the `SNES` context Output Parameter: . force - `PETSC_TRUE` requires at least one iteration. Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESSetForceIteration()`, `SNESSetDivergenceTolerance()` @*/ PetscErrorCode SNESGetForceIteration(SNES snes, PetscBool *force) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); *force = snes->forceiteration; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetTolerances - Sets various parameters used in `SNES` convergence tests. Logically Collective Input Parameters: + snes - the `SNES` context . abstol - the absolute convergence tolerance, $ F(x^n) \le abstol $ . rtol - the relative convergence tolerance, $ F(x^n) \le reltol * F(x^0) $ . stol - convergence tolerance in terms of the norm of the change in the solution between steps, || delta x || < stol*|| x || . maxit - the maximum number of iterations allowed in the solver, default 50. - maxf - the maximum number of function evaluations allowed in the solver (use `PETSC_UNLIMITED` indicates no limit), default 10,000 Options Database Keys: + -snes_atol - Sets `abstol` . -snes_rtol - Sets `rtol` . -snes_stol - Sets `stol` . -snes_max_it - Sets `maxit` - -snes_max_funcs - Sets `maxf` (use `unlimited` to have no maximum) Level: intermediate Note: All parameters must be non-negative Use `PETSC_CURRENT` to retain the current value of any parameter and `PETSC_DETERMINE` to use the default value for the given `SNES`. The default value is the value in the object when its type is set. Use `PETSC_UNLIMITED` on `maxit` or `maxf` to indicate there is no bound on the number of iterations or number of function evaluations. Fortran Note: Use `PETSC_CURRENT_INTEGER`, `PETSC_CURRENT_REAL`, `PETSC_UNLIMITED_INTEGER`, `PETSC_DETERMINE_INTEGER`, or `PETSC_DETERMINE_REAL` .seealso: [](ch_snes), `SNESSolve()`, `SNES`, `SNESSetDivergenceTolerance()`, `SNESSetForceIteration()` @*/ PetscErrorCode SNESSetTolerances(SNES snes, PetscReal abstol, PetscReal rtol, PetscReal stol, PetscInt maxit, PetscInt maxf) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveReal(snes, abstol, 2); PetscValidLogicalCollectiveReal(snes, rtol, 3); PetscValidLogicalCollectiveReal(snes, stol, 4); PetscValidLogicalCollectiveInt(snes, maxit, 5); PetscValidLogicalCollectiveInt(snes, maxf, 6); if (abstol == (PetscReal)PETSC_DETERMINE) { snes->abstol = snes->default_abstol; } else if (abstol != (PetscReal)PETSC_CURRENT) { PetscCheck(abstol >= 0.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Absolute tolerance %g must be non-negative", (double)abstol); snes->abstol = abstol; } if (rtol == (PetscReal)PETSC_DETERMINE) { snes->rtol = snes->default_rtol; } else if (rtol != (PetscReal)PETSC_CURRENT) { PetscCheck(rtol >= 0.0 && 1.0 > rtol, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Relative tolerance %g must be non-negative and less than 1.0", (double)rtol); snes->rtol = rtol; } if (stol == (PetscReal)PETSC_DETERMINE) { snes->stol = snes->default_stol; } else if (stol != (PetscReal)PETSC_CURRENT) { PetscCheck(stol >= 0.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Step tolerance %g must be non-negative", (double)stol); snes->stol = stol; } if (maxit == PETSC_DETERMINE) { snes->max_its = snes->default_max_its; } else if (maxit == PETSC_UNLIMITED) { snes->max_its = PETSC_INT_MAX; } else if (maxit != PETSC_CURRENT) { PetscCheck(maxit >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of iterations %" PetscInt_FMT " must be non-negative", maxit); snes->max_its = maxit; } if (maxf == PETSC_DETERMINE) { snes->max_funcs = snes->default_max_funcs; } else if (maxf == PETSC_UNLIMITED || maxf == -1) { snes->max_funcs = PETSC_UNLIMITED; } else if (maxf != PETSC_CURRENT) { PetscCheck(maxf >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Maximum number of function evaluations %" PetscInt_FMT " must be nonnegative", maxf); snes->max_funcs = maxf; } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetDivergenceTolerance - Sets the divergence tolerance used for the `SNES` divergence test. Logically Collective Input Parameters: + snes - the `SNES` context - divtol - the divergence tolerance. Use `PETSC_UNLIMITED` to deactivate the test. If the residual norm $ F(x^n) \ge divtol * F(x^0) $ the solver is stopped due to divergence. Options Database Key: . -snes_divergence_tolerance - Sets `divtol` Level: intermediate Notes: Use `PETSC_DETERMINE` to use the default value from when the object's type was set. Fortran Note: Use ``PETSC_DETERMINE_REAL` or `PETSC_UNLIMITED_REAL` .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESSetTolerances()`, `SNESGetDivergenceTolerance()` @*/ PetscErrorCode SNESSetDivergenceTolerance(SNES snes, PetscReal divtol) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveReal(snes, divtol, 2); if (divtol == (PetscReal)PETSC_DETERMINE) { snes->divtol = snes->default_divtol; } else if (divtol == (PetscReal)PETSC_UNLIMITED || divtol == -1) { snes->divtol = PETSC_UNLIMITED; } else if (divtol != (PetscReal)PETSC_CURRENT) { PetscCheck(divtol >= 1.0, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_OUTOFRANGE, "Divergence tolerance %g must be greater than 1.0", (double)divtol); snes->divtol = divtol; } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetTolerances - Gets various parameters used in `SNES` convergence tests. Not Collective Input Parameter: . snes - the `SNES` context Output Parameters: + atol - the absolute convergence tolerance . rtol - the relative convergence tolerance . stol - convergence tolerance in terms of the norm of the change in the solution between steps . maxit - the maximum number of iterations allowed - maxf - the maximum number of function evaluations allowed, `PETSC_UNLIMITED` indicates no bound Level: intermediate Notes: See `SNESSetTolerances()` for details on the parameters. The user can specify `NULL` for any parameter that is not needed. .seealso: [](ch_snes), `SNES`, `SNESSetTolerances()` @*/ PetscErrorCode SNESGetTolerances(SNES snes, PetscReal *atol, PetscReal *rtol, PetscReal *stol, PetscInt *maxit, PetscInt *maxf) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (atol) *atol = snes->abstol; if (rtol) *rtol = snes->rtol; if (stol) *stol = snes->stol; if (maxit) *maxit = snes->max_its; if (maxf) *maxf = snes->max_funcs; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetDivergenceTolerance - Gets divergence tolerance used in divergence test. Not Collective Input Parameters: + snes - the `SNES` context - divtol - divergence tolerance Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESSetDivergenceTolerance()` @*/ PetscErrorCode SNESGetDivergenceTolerance(SNES snes, PetscReal *divtol) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (divtol) *divtol = snes->divtol; PetscFunctionReturn(PETSC_SUCCESS); } PETSC_INTERN PetscErrorCode SNESMonitorRange_Private(SNES, PetscInt, PetscReal *); PetscErrorCode SNESMonitorLGRange(SNES snes, PetscInt n, PetscReal rnorm, void *monctx) { PetscDrawLG lg; PetscReal x, y, per; PetscViewer v = (PetscViewer)monctx; static PetscReal prev; /* should be in the context */ PetscDraw draw; PetscFunctionBegin; PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 4); PetscCall(PetscViewerDrawGetDrawLG(v, 0, &lg)); if (!n) PetscCall(PetscDrawLGReset(lg)); PetscCall(PetscDrawLGGetDraw(lg, &draw)); PetscCall(PetscDrawSetTitle(draw, "Residual norm")); x = (PetscReal)n; if (rnorm > 0.0) y = PetscLog10Real(rnorm); else y = -15.0; PetscCall(PetscDrawLGAddPoint(lg, &x, &y)); if (n < 20 || !(n % 5) || snes->reason) { PetscCall(PetscDrawLGDraw(lg)); PetscCall(PetscDrawLGSave(lg)); } PetscCall(PetscViewerDrawGetDrawLG(v, 1, &lg)); if (!n) PetscCall(PetscDrawLGReset(lg)); PetscCall(PetscDrawLGGetDraw(lg, &draw)); PetscCall(PetscDrawSetTitle(draw, "% elements > .2*max element")); PetscCall(SNESMonitorRange_Private(snes, n, &per)); x = (PetscReal)n; y = 100.0 * per; PetscCall(PetscDrawLGAddPoint(lg, &x, &y)); if (n < 20 || !(n % 5) || snes->reason) { PetscCall(PetscDrawLGDraw(lg)); PetscCall(PetscDrawLGSave(lg)); } PetscCall(PetscViewerDrawGetDrawLG(v, 2, &lg)); if (!n) { prev = rnorm; PetscCall(PetscDrawLGReset(lg)); } PetscCall(PetscDrawLGGetDraw(lg, &draw)); PetscCall(PetscDrawSetTitle(draw, "(norm -oldnorm)/oldnorm")); x = (PetscReal)n; y = (prev - rnorm) / prev; PetscCall(PetscDrawLGAddPoint(lg, &x, &y)); if (n < 20 || !(n % 5) || snes->reason) { PetscCall(PetscDrawLGDraw(lg)); PetscCall(PetscDrawLGSave(lg)); } PetscCall(PetscViewerDrawGetDrawLG(v, 3, &lg)); if (!n) PetscCall(PetscDrawLGReset(lg)); PetscCall(PetscDrawLGGetDraw(lg, &draw)); PetscCall(PetscDrawSetTitle(draw, "(norm -oldnorm)/oldnorm*(% > .2 max)")); x = (PetscReal)n; y = (prev - rnorm) / (prev * per); if (n > 2) { /*skip initial crazy value */ PetscCall(PetscDrawLGAddPoint(lg, &x, &y)); } if (n < 20 || !(n % 5) || snes->reason) { PetscCall(PetscDrawLGDraw(lg)); PetscCall(PetscDrawLGSave(lg)); } prev = rnorm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESConverged - Run the convergence test and update the `SNESConvergedReason`. Collective Input Parameters: + snes - the `SNES` context . it - current iteration . xnorm - 2-norm of current iterate . snorm - 2-norm of current step - fnorm - 2-norm of function Level: developer Note: This routine is called by the `SNESSolve()` implementations. It does not typically need to be called by the user. .seealso: [](ch_snes), `SNES`, `SNESSolve`, `SNESSetConvergenceTest()` @*/ PetscErrorCode SNESConverged(SNES snes, PetscInt it, PetscReal xnorm, PetscReal snorm, PetscReal fnorm) { PetscFunctionBegin; if (!snes->reason) { if (snes->normschedule == SNES_NORM_ALWAYS) PetscUseTypeMethod(snes, converged, it, xnorm, snorm, fnorm, &snes->reason, snes->cnvP); if (it == snes->max_its && !snes->reason) { if (snes->normschedule == SNES_NORM_ALWAYS) { PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its)); snes->reason = SNES_DIVERGED_MAX_IT; } else snes->reason = SNES_CONVERGED_ITS; } } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESMonitor - runs any `SNES` monitor routines provided with `SNESMonitor()` or the options database Collective Input Parameters: + snes - nonlinear solver context obtained from `SNESCreate()` . iter - current iteration number - rnorm - current relative norm of the residual Level: developer Note: This routine is called by the `SNESSolve()` implementations. It does not typically need to be called by the user. .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()` @*/ PetscErrorCode SNESMonitor(SNES snes, PetscInt iter, PetscReal rnorm) { PetscInt i, n = snes->numbermonitors; PetscFunctionBegin; PetscCall(VecLockReadPush(snes->vec_sol)); for (i = 0; i < n; i++) PetscCall((*snes->monitor[i])(snes, iter, rnorm, snes->monitorcontext[i])); PetscCall(VecLockReadPop(snes->vec_sol)); PetscFunctionReturn(PETSC_SUCCESS); } /* ------------ Routines to set performance monitoring options ----------- */ /*MC SNESMonitorFunction - functional form passed to `SNESMonitorSet()` to monitor convergence of nonlinear solver Synopsis: #include PetscErrorCode SNESMonitorFunction(SNES snes, PetscInt its, PetscReal norm, void *mctx) Collective Input Parameters: + snes - the `SNES` context . its - iteration number . norm - 2-norm function value (may be estimated) - mctx - [optional] monitoring context Level: advanced .seealso: [](ch_snes), `SNESMonitorSet()` M*/ /*@C SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every iteration of the `SNES` nonlinear solver to display the iteration's progress. Logically Collective Input Parameters: + snes - the `SNES` context . f - the monitor function, for the calling sequence see `SNESMonitorFunction` . mctx - [optional] user-defined context for private data for the monitor routine (use `NULL` if no context is desired) - monitordestroy - [optional] routine that frees monitor context (may be `NULL`), see `PetscCtxDestroyFn` for the calling sequence Options Database Keys: + -snes_monitor - sets `SNESMonitorDefault()` . -snes_monitor draw::draw_lg - sets line graph monitor, - -snes_monitor_cancel - cancels all monitors that have been hardwired into a code by calls to `SNESMonitorSet()`, but does not cancel those set via the options database. Level: intermediate Note: Several different monitoring routines may be set by calling `SNESMonitorSet()` multiple times; all will be called in the order in which they were set. Fortran Note: Only a single monitor function can be set for each `SNES` object .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESMonitorDefault()`, `SNESMonitorCancel()`, `SNESMonitorFunction`, `PetscCtxDestroyFn` @*/ PetscErrorCode SNESMonitorSet(SNES snes, PetscErrorCode (*f)(SNES, PetscInt, PetscReal, void *), void *mctx, PetscCtxDestroyFn *monitordestroy) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); for (PetscInt i = 0; i < snes->numbermonitors; i++) { PetscBool identical; PetscCall(PetscMonitorCompare((PetscErrorCode (*)(void))(PetscVoidFn *)f, mctx, monitordestroy, (PetscErrorCode (*)(void))(PetscVoidFn *)snes->monitor[i], snes->monitorcontext[i], snes->monitordestroy[i], &identical)); if (identical) PetscFunctionReturn(PETSC_SUCCESS); } PetscCheck(snes->numbermonitors < MAXSNESMONITORS, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many monitors set"); snes->monitor[snes->numbermonitors] = f; snes->monitordestroy[snes->numbermonitors] = monitordestroy; snes->monitorcontext[snes->numbermonitors++] = mctx; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESMonitorCancel - Clears all the monitor functions for a `SNES` object. Logically Collective Input Parameter: . snes - the `SNES` context Options Database Key: . -snes_monitor_cancel - cancels all monitors that have been hardwired into a code by calls to `SNESMonitorSet()`, but does not cancel those set via the options database Level: intermediate Note: There is no way to clear one specific monitor from a `SNES` object. .seealso: [](ch_snes), `SNES`, `SNESMonitorDefault()`, `SNESMonitorSet()` @*/ PetscErrorCode SNESMonitorCancel(SNES snes) { PetscInt i; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); for (i = 0; i < snes->numbermonitors; i++) { if (snes->monitordestroy[i]) PetscCall((*snes->monitordestroy[i])(&snes->monitorcontext[i])); } snes->numbermonitors = 0; PetscFunctionReturn(PETSC_SUCCESS); } /*MC SNESConvergenceTestFunction - functional form used for testing of convergence of nonlinear solver Synopsis: #include PetscErrorCode SNESConvergenceTest(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *cctx) Collective Input Parameters: + snes - the `SNES` context . it - current iteration (0 is the first and is before any Newton step) . xnorm - 2-norm of current iterate . gnorm - 2-norm of current step . f - 2-norm of function - cctx - [optional] convergence context Output Parameter: . reason - reason for convergence/divergence, only needs to be set when convergence or divergence is detected Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESSolve`, `SNESSetConvergenceTest()` M*/ /*@C SNESSetConvergenceTest - Sets the function that is to be used to test for convergence of the nonlinear iterative solution. Logically Collective Input Parameters: + snes - the `SNES` context . SNESConvergenceTestFunction - routine to test for convergence . cctx - [optional] context for private data for the convergence routine (may be `NULL`) - destroy - [optional] destructor for the context (may be `NULL`; `PETSC_NULL_FUNCTION` in Fortran) Level: advanced .seealso: [](ch_snes), `SNES`, `SNESConvergedDefault()`, `SNESConvergedSkip()`, `SNESConvergenceTestFunction` @*/ PetscErrorCode SNESSetConvergenceTest(SNES snes, PetscErrorCode (*SNESConvergenceTestFunction)(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *), void *cctx, PetscCtxDestroyFn *destroy) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (!SNESConvergenceTestFunction) SNESConvergenceTestFunction = SNESConvergedSkip; if (snes->ops->convergeddestroy) PetscCall((*snes->ops->convergeddestroy)(&snes->cnvP)); snes->ops->converged = SNESConvergenceTestFunction; snes->ops->convergeddestroy = destroy; snes->cnvP = cctx; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetConvergedReason - Gets the reason the `SNES` iteration was stopped, which may be due to convergence, divergence, or stagnation Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . reason - negative value indicates diverged, positive value converged, see `SNESConvergedReason` for the individual convergence tests for complete lists Options Database Key: . -snes_converged_reason - prints the reason to standard out Level: intermediate Note: Should only be called after the call the `SNESSolve()` is complete, if it is called earlier it returns the value `SNES__CONVERGED_ITERATING`. .seealso: [](ch_snes), `SNESSolve()`, `SNESSetConvergenceTest()`, `SNESSetConvergedReason()`, `SNESConvergedReason`, `SNESGetConvergedReasonString()` @*/ PetscErrorCode SNESGetConvergedReason(SNES snes, SNESConvergedReason *reason) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(reason, 2); *reason = snes->reason; PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESGetConvergedReasonString - Return a human readable string for `SNESConvergedReason` Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . strreason - a human readable string that describes `SNES` converged reason Level: beginner .seealso: [](ch_snes), `SNES`, `SNESGetConvergedReason()` @*/ PetscErrorCode SNESGetConvergedReasonString(SNES snes, const char **strreason) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(strreason, 2); *strreason = SNESConvergedReasons[snes->reason]; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetConvergedReason - Sets the reason the `SNES` iteration was stopped. Not Collective Input Parameters: + snes - the `SNES` context - reason - negative value indicates diverged, positive value converged, see `SNESConvergedReason` or the manual pages for the individual convergence tests for complete lists Level: developer Developer Note: Called inside the various `SNESSolve()` implementations .seealso: [](ch_snes), `SNESGetConvergedReason()`, `SNESSetConvergenceTest()`, `SNESConvergedReason` @*/ PetscErrorCode SNESSetConvergedReason(SNES snes, SNESConvergedReason reason) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCheck(!snes->errorifnotconverged || reason > 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "SNES code should have previously errored due to negative reason"); snes->reason = reason; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetConvergenceHistory - Sets the arrays used to hold the convergence history. Logically Collective Input Parameters: + snes - iterative context obtained from `SNESCreate()` . a - array to hold history, this array will contain the function norms computed at each step . its - integer array holds the number of linear iterations for each solve. . na - size of `a` and `its` - reset - `PETSC_TRUE` indicates each new nonlinear solve resets the history counter to zero, else it continues storing new values for new nonlinear solves after the old ones Level: intermediate Notes: If 'a' and 'its' are `NULL` then space is allocated for the history. If 'na' is `PETSC_DECIDE` (or, deprecated, `PETSC_DEFAULT`) then a default array of length 1,000 is allocated. This routine is useful, e.g., when running a code for purposes of accurate performance monitoring, when no I/O should be done during the section of code that is being timed. If the arrays run out of space after a number of iterations then the later values are not saved in the history .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESGetConvergenceHistory()` @*/ PetscErrorCode SNESSetConvergenceHistory(SNES snes, PetscReal a[], PetscInt its[], PetscInt na, PetscBool reset) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (a) PetscAssertPointer(a, 2); if (its) PetscAssertPointer(its, 3); if (!a) { if (na == PETSC_DECIDE) na = 1000; PetscCall(PetscCalloc2(na, &a, na, &its)); snes->conv_hist_alloc = PETSC_TRUE; } snes->conv_hist = a; snes->conv_hist_its = its; snes->conv_hist_max = (size_t)na; snes->conv_hist_len = 0; snes->conv_hist_reset = reset; PetscFunctionReturn(PETSC_SUCCESS); } #if defined(PETSC_HAVE_MATLAB) #include /* MATLAB include file */ #include /* MATLAB include file */ PETSC_EXTERN mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) { mxArray *mat; PetscInt i; PetscReal *ar; mat = mxCreateDoubleMatrix(snes->conv_hist_len, 1, mxREAL); ar = (PetscReal *)mxGetData(mat); for (i = 0; i < snes->conv_hist_len; i++) ar[i] = snes->conv_hist[i]; return mat; } #endif /*@C SNESGetConvergenceHistory - Gets the arrays used to hold the convergence history. Not Collective Input Parameter: . snes - iterative context obtained from `SNESCreate()` Output Parameters: + a - array to hold history, usually was set with `SNESSetConvergenceHistory()` . its - integer array holds the number of linear iterations (or negative if not converged) for each solve. - na - size of `a` and `its` Level: intermediate Note: This routine is useful, e.g., when running a code for purposes of accurate performance monitoring, when no I/O should be done during the section of code that is being timed. Fortran Notes: Return the arrays with ``SNESRestoreConvergenceHistory()` Use the arguments .vb PetscReal, pointer :: a(:) PetscInt, pointer :: its(:) .ve .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESSetConvergenceHistory()` @*/ PetscErrorCode SNESGetConvergenceHistory(SNES snes, PetscReal *a[], PetscInt *its[], PetscInt *na) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (a) *a = snes->conv_hist; if (its) *its = snes->conv_hist_its; if (na) *na = (PetscInt)snes->conv_hist_len; PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESSetUpdate - Sets the general-purpose update function called at the beginning of every iteration of the nonlinear solve. Specifically it is called just before the Jacobian is "evaluated" and after the function evaluation. Logically Collective Input Parameters: + snes - The nonlinear solver context - func - The update function; for calling sequence see `SNESUpdateFn` Level: advanced Notes: This is NOT what one uses to update the ghost points before a function evaluation, that should be done at the beginning of your function provided to `SNESSetFunction()`, or `SNESSetPicard()` This is not used by most users, and it is intended to provide a general hook that is run right before the direction step is computed. Users are free to modify the current residual vector, the current linearization point, or any other vector associated to the specific solver used. If such modifications take place, it is the user responsibility to update all the relevant vectors. For example, if one is adjusting the model parameters at each Newton step their code may look like .vb PetscErrorCode update(SNES snes, PetscInt iteration) { PetscFunctionBeginUser; if (iteration > 0) { // update the model parameters here Vec x,f; PetscCall(SNESGetSolution(snes,&x)); PetcCall(SNESGetFunction(snes,&f,NULL,NULL)); PetscCall(SNESComputeFunction(snes,x,f)); } PetscFunctionReturn(PETSC_SUCCESS); } .ve There are a variety of function hooks one many set that are called at different stages of the nonlinear solution process, see the functions listed below. .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESSetJacobian()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchSetPostCheck()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESMonitorSet()` @*/ PetscErrorCode SNESSetUpdate(SNES snes, SNESUpdateFn *func) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); snes->ops->update = func; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESConvergedReasonView - Displays the reason a `SNES` solve converged or diverged to a viewer Collective Input Parameters: + snes - iterative context obtained from `SNESCreate()` - viewer - the viewer to display the reason Options Database Keys: + -snes_converged_reason - print reason for converged or diverged, also prints number of iterations - -snes_converged_reason ::failed - only print reason and number of iterations when diverged Level: beginner Note: To change the format of the output call `PetscViewerPushFormat`(viewer,format) before this call. Use `PETSC_VIEWER_DEFAULT` for the default, use `PETSC_VIEWER_FAILED` to only display a reason if it fails. .seealso: [](ch_snes), `SNESConvergedReason`, `PetscViewer`, `SNES`, `SNESCreate()`, `SNESSetUp()`, `SNESDestroy()`, `SNESSetTolerances()`, `SNESConvergedDefault()`, `SNESGetConvergedReason()`, `SNESConvergedReasonViewFromOptions()`, `PetscViewerPushFormat()`, `PetscViewerPopFormat()` @*/ PetscErrorCode SNESConvergedReasonView(SNES snes, PetscViewer viewer) { PetscViewerFormat format; PetscBool isAscii; PetscFunctionBegin; if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)); PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isAscii)); if (isAscii) { PetscCall(PetscViewerGetFormat(viewer, &format)); PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel + 1)); if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { DM dm; Vec u; PetscDS prob; PetscInt Nf, f; PetscErrorCode (**exactSol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *); void **exactCtx; PetscReal error; PetscCall(SNESGetDM(snes, &dm)); PetscCall(SNESGetSolution(snes, &u)); PetscCall(DMGetDS(dm, &prob)); PetscCall(PetscDSGetNumFields(prob, &Nf)); PetscCall(PetscMalloc2(Nf, &exactSol, Nf, &exactCtx)); for (f = 0; f < Nf; ++f) PetscCall(PetscDSGetExactSolution(prob, f, &exactSol[f], &exactCtx[f])); PetscCall(DMComputeL2Diff(dm, 0.0, exactSol, exactCtx, u, &error)); PetscCall(PetscFree2(exactSol, exactCtx)); if (error < 1.0e-11) PetscCall(PetscViewerASCIIPrintf(viewer, "L_2 Error: < 1.0e-11\n")); else PetscCall(PetscViewerASCIIPrintf(viewer, "L_2 Error: %g\n", (double)error)); } if (snes->reason > 0 && format != PETSC_VIEWER_FAILED) { if (((PetscObject)snes)->prefix) { PetscCall(PetscViewerASCIIPrintf(viewer, "Nonlinear %s solve converged due to %s iterations %" PetscInt_FMT "\n", ((PetscObject)snes)->prefix, SNESConvergedReasons[snes->reason], snes->iter)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, "Nonlinear solve converged due to %s iterations %" PetscInt_FMT "\n", SNESConvergedReasons[snes->reason], snes->iter)); } } else if (snes->reason <= 0) { if (((PetscObject)snes)->prefix) { PetscCall(PetscViewerASCIIPrintf(viewer, "Nonlinear %s solve did not converge due to %s iterations %" PetscInt_FMT "\n", ((PetscObject)snes)->prefix, SNESConvergedReasons[snes->reason], snes->iter)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, "Nonlinear solve did not converge due to %s iterations %" PetscInt_FMT "\n", SNESConvergedReasons[snes->reason], snes->iter)); } } PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel + 1)); } PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESConvergedReasonViewSet - Sets an ADDITIONAL function that is to be used at the end of the nonlinear solver to display the convergence reason of the nonlinear solver. Logically Collective Input Parameters: + snes - the `SNES` context . f - the `SNESConvergedReason` view function . vctx - [optional] user-defined context for private data for the `SNESConvergedReason` view function (use `NULL` if no context is desired) - reasonviewdestroy - [optional] routine that frees the context (may be `NULL`), see `PetscCtxDestroyFn` for the calling sequence Calling sequence of `f`: + snes - the `SNES` context - vctx - [optional] context for private data for the function Options Database Keys: + -snes_converged_reason - sets a default `SNESConvergedReasonView()` - -snes_converged_reason_view_cancel - cancels all converged reason viewers that have been hardwired into a code by calls to `SNESConvergedReasonViewSet()`, but does not cancel those set via the options database. Level: intermediate Note: Several different converged reason view routines may be set by calling `SNESConvergedReasonViewSet()` multiple times; all will be called in the order in which they were set. .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESConvergedReason`, `SNESGetConvergedReason()`, `SNESConvergedReasonView()`, `SNESConvergedReasonViewCancel()`, `PetscCtxDestroyFn` @*/ PetscErrorCode SNESConvergedReasonViewSet(SNES snes, PetscErrorCode (*f)(SNES snes, void *vctx), void *vctx, PetscCtxDestroyFn *reasonviewdestroy) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); for (PetscInt i = 0; i < snes->numberreasonviews; i++) { PetscBool identical; PetscCall(PetscMonitorCompare((PetscErrorCode (*)(void))(PetscVoidFn *)f, vctx, reasonviewdestroy, (PetscErrorCode (*)(void))(PetscVoidFn *)snes->reasonview[i], snes->reasonviewcontext[i], snes->reasonviewdestroy[i], &identical)); if (identical) PetscFunctionReturn(PETSC_SUCCESS); } PetscCheck(snes->numberreasonviews < MAXSNESREASONVIEWS, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many SNES reasonview set"); snes->reasonview[snes->numberreasonviews] = f; snes->reasonviewdestroy[snes->numberreasonviews] = reasonviewdestroy; snes->reasonviewcontext[snes->numberreasonviews++] = vctx; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESConvergedReasonViewFromOptions - Processes command line options to determine if/how a `SNESConvergedReason` is to be viewed at the end of `SNESSolve()` All the user-provided viewer routines set with `SNESConvergedReasonViewSet()` will be called, if they exist. Collective Input Parameter: . snes - the `SNES` object Level: advanced .seealso: [](ch_snes), `SNES`, `SNESConvergedReason`, `SNESConvergedReasonViewSet()`, `SNESCreate()`, `SNESSetUp()`, `SNESDestroy()`, `SNESSetTolerances()`, `SNESConvergedDefault()`, `SNESGetConvergedReason()`, `SNESConvergedReasonView()` @*/ PetscErrorCode SNESConvergedReasonViewFromOptions(SNES snes) { static PetscBool incall = PETSC_FALSE; PetscFunctionBegin; if (incall) PetscFunctionReturn(PETSC_SUCCESS); incall = PETSC_TRUE; /* All user-provided viewers are called first, if they exist. */ for (PetscInt i = 0; i < snes->numberreasonviews; i++) PetscCall((*snes->reasonview[i])(snes, snes->reasonviewcontext[i])); /* Call PETSc default routine if users ask for it */ if (snes->convergedreasonviewer) { PetscCall(PetscViewerPushFormat(snes->convergedreasonviewer, snes->convergedreasonformat)); PetscCall(SNESConvergedReasonView(snes, snes->convergedreasonviewer)); PetscCall(PetscViewerPopFormat(snes->convergedreasonviewer)); } incall = PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSolve - Solves a nonlinear system $F(x) = b $ associated with a `SNES` object Collective Input Parameters: + snes - the `SNES` context . b - the constant part of the equation $F(x) = b$, or `NULL` to use zero. - x - the solution vector. Level: beginner Note: The user should initialize the vector, `x`, with the initial guess for the nonlinear solve prior to calling `SNESSolve()` . .seealso: [](ch_snes), `SNES`, `SNESCreate()`, `SNESDestroy()`, `SNESSetFunction()`, `SNESSetJacobian()`, `SNESSetGridSequence()`, `SNESGetSolution()`, `SNESNewtonTRSetPreCheck()`, `SNESNewtonTRGetPreCheck()`, `SNESNewtonTRSetPostCheck()`, `SNESNewtonTRGetPostCheck()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchGetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESLineSearchGetPreCheck()` @*/ PetscErrorCode SNESSolve(SNES snes, Vec b, Vec x) { PetscBool flg; PetscInt grid; Vec xcreated = NULL; DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (x) PetscValidHeaderSpecific(x, VEC_CLASSID, 3); if (x) PetscCheckSameComm(snes, 1, x, 3); if (b) PetscValidHeaderSpecific(b, VEC_CLASSID, 2); if (b) PetscCheckSameComm(snes, 1, b, 2); /* High level operations using the nonlinear solver */ { PetscViewer viewer; PetscViewerFormat format; PetscInt num; PetscBool flg; static PetscBool incall = PETSC_FALSE; if (!incall) { /* Estimate the convergence rate of the discretization */ PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_convergence_estimate", &viewer, &format, &flg)); if (flg) { PetscConvEst conv; DM dm; PetscReal *alpha; /* Convergence rate of the solution error for each field in the L_2 norm */ PetscInt Nf; incall = PETSC_TRUE; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMGetNumFields(dm, &Nf)); PetscCall(PetscCalloc1(Nf, &alpha)); PetscCall(PetscConvEstCreate(PetscObjectComm((PetscObject)snes), &conv)); PetscCall(PetscConvEstSetSolver(conv, (PetscObject)snes)); PetscCall(PetscConvEstSetFromOptions(conv)); PetscCall(PetscConvEstSetUp(conv)); PetscCall(PetscConvEstGetConvRate(conv, alpha)); PetscCall(PetscViewerPushFormat(viewer, format)); PetscCall(PetscConvEstRateView(conv, alpha, viewer)); PetscCall(PetscViewerPopFormat(viewer)); PetscCall(PetscViewerDestroy(&viewer)); PetscCall(PetscConvEstDestroy(&conv)); PetscCall(PetscFree(alpha)); incall = PETSC_FALSE; } /* Adaptively refine the initial grid */ num = 1; PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)snes)->prefix, "-snes_adapt_initial", &num, &flg)); if (flg) { DMAdaptor adaptor; incall = PETSC_TRUE; PetscCall(DMAdaptorCreate(PetscObjectComm((PetscObject)snes), &adaptor)); PetscCall(DMAdaptorSetSolver(adaptor, snes)); PetscCall(DMAdaptorSetSequenceLength(adaptor, num)); PetscCall(DMAdaptorSetFromOptions(adaptor)); PetscCall(DMAdaptorSetUp(adaptor)); PetscCall(DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_INITIAL, &dm, &x)); PetscCall(DMAdaptorDestroy(&adaptor)); incall = PETSC_FALSE; } /* Use grid sequencing to adapt */ num = 0; PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)snes)->prefix, "-snes_adapt_sequence", &num, NULL)); if (num) { DMAdaptor adaptor; const char *prefix; incall = PETSC_TRUE; PetscCall(DMAdaptorCreate(PetscObjectComm((PetscObject)snes), &adaptor)); PetscCall(SNESGetOptionsPrefix(snes, &prefix)); PetscCall(DMAdaptorSetOptionsPrefix(adaptor, prefix)); PetscCall(DMAdaptorSetSolver(adaptor, snes)); PetscCall(DMAdaptorSetSequenceLength(adaptor, num)); PetscCall(DMAdaptorSetFromOptions(adaptor)); PetscCall(DMAdaptorSetUp(adaptor)); PetscCall(PetscObjectViewFromOptions((PetscObject)adaptor, NULL, "-snes_adapt_view")); PetscCall(DMAdaptorAdapt(adaptor, x, DM_ADAPTATION_SEQUENTIAL, &dm, &x)); PetscCall(DMAdaptorDestroy(&adaptor)); incall = PETSC_FALSE; } } } if (!x) x = snes->vec_sol; if (!x) { PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMCreateGlobalVector(dm, &xcreated)); x = xcreated; } PetscCall(SNESViewFromOptions(snes, NULL, "-snes_view_pre")); for (grid = 0; grid < snes->gridsequence; grid++) PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)))); for (grid = 0; grid < snes->gridsequence + 1; grid++) { /* set solution vector */ if (!grid) PetscCall(PetscObjectReference((PetscObject)x)); PetscCall(VecDestroy(&snes->vec_sol)); snes->vec_sol = x; PetscCall(SNESGetDM(snes, &dm)); /* set affine vector if provided */ if (b) PetscCall(PetscObjectReference((PetscObject)b)); PetscCall(VecDestroy(&snes->vec_rhs)); snes->vec_rhs = b; if (snes->vec_rhs) PetscCheck(snes->vec_func != snes->vec_rhs, PETSC_COMM_SELF, PETSC_ERR_ARG_IDN, "Right hand side vector cannot be function vector"); PetscCheck(snes->vec_func != snes->vec_sol, PETSC_COMM_SELF, PETSC_ERR_ARG_IDN, "Solution vector cannot be function vector"); PetscCheck(snes->vec_rhs != snes->vec_sol, PETSC_COMM_SELF, PETSC_ERR_ARG_IDN, "Solution vector cannot be right-hand side vector"); if (!snes->vec_sol_update /* && snes->vec_sol */) PetscCall(VecDuplicate(snes->vec_sol, &snes->vec_sol_update)); PetscCall(DMShellSetGlobalVector(dm, snes->vec_sol)); PetscCall(SNESSetUp(snes)); if (!grid) { if (snes->ops->computeinitialguess) PetscCallBack("SNES callback compute initial guess", (*snes->ops->computeinitialguess)(snes, snes->vec_sol, snes->initialguessP)); } if (snes->conv_hist_reset) snes->conv_hist_len = 0; PetscCall(SNESResetCounters(snes)); snes->reason = SNES_CONVERGED_ITERATING; PetscCall(PetscLogEventBegin(SNES_Solve, snes, 0, 0, 0)); PetscUseTypeMethod(snes, solve); PetscCall(PetscLogEventEnd(SNES_Solve, snes, 0, 0, 0)); PetscCheck(snes->reason, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Internal error, solver %s returned without setting converged reason", ((PetscObject)snes)->type_name); snes->domainerror = PETSC_FALSE; /* clear the flag if it has been set */ if (snes->lagjac_persist) snes->jac_iter += snes->iter; if (snes->lagpre_persist) snes->pre_iter += snes->iter; PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_test_local_min", NULL, NULL, &flg)); if (flg && !PetscPreLoadingOn) PetscCall(SNESTestLocalMin(snes)); /* Call converged reason views. This may involve user-provided viewers as well */ PetscCall(SNESConvergedReasonViewFromOptions(snes)); if (snes->errorifnotconverged) { if (snes->reason < 0) PetscCall(SNESMonitorCancel(snes)); PetscCheck(snes->reason >= 0, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNESSolve has not converged"); } if (snes->reason < 0) break; if (grid < snes->gridsequence) { DM fine; Vec xnew; Mat interp; PetscCall(DMRefine(snes->dm, PetscObjectComm((PetscObject)snes), &fine)); PetscCheck(fine, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing"); PetscCall(DMGetCoordinatesLocalSetUp(fine)); PetscCall(DMCreateInterpolation(snes->dm, fine, &interp, NULL)); PetscCall(DMCreateGlobalVector(fine, &xnew)); PetscCall(MatInterpolate(interp, x, xnew)); PetscCall(DMInterpolate(snes->dm, interp, fine)); PetscCall(MatDestroy(&interp)); x = xnew; PetscCall(SNESReset(snes)); PetscCall(SNESSetDM(snes, fine)); PetscCall(SNESResetFromOptions(snes)); PetscCall(DMDestroy(&fine)); PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)))); } } PetscCall(SNESViewFromOptions(snes, NULL, "-snes_view")); PetscCall(VecViewFromOptions(snes->vec_sol, (PetscObject)snes, "-snes_view_solution")); PetscCall(DMMonitor(snes->dm)); PetscCall(SNESMonitorPauseFinal_Internal(snes)); PetscCall(VecDestroy(&xcreated)); PetscCall(PetscObjectSAWsBlock((PetscObject)snes)); PetscFunctionReturn(PETSC_SUCCESS); } /* --------- Internal routines for SNES Package --------- */ /*@ SNESSetType - Sets the algorithm/method to be used to solve the nonlinear system with the given `SNES` Collective Input Parameters: + snes - the `SNES` context - type - a known method Options Database Key: . -snes_type - Sets the method; use -help for a list of available methods (for instance, newtonls or newtontr) Level: intermediate Notes: See `SNESType` for available methods (for instance) + `SNESNEWTONLS` - Newton's method with line search (systems of nonlinear equations) - `SNESNEWTONTR` - Newton's method with trust region (systems of nonlinear equations) Normally, it is best to use the `SNESSetFromOptions()` command and then set the `SNES` solver type from the options database rather than by using this routine. Using the options database provides the user with maximum flexibility in evaluating the many nonlinear solvers. The `SNESSetType()` routine is provided for those situations where it is necessary to set the nonlinear solver independently of the command line or options database. This might be the case, for example, when the choice of solver changes during the execution of the program, and the user's application is taking responsibility for choosing the appropriate method. Developer Note: `SNESRegister()` adds a constructor for a new `SNESType` to `SNESList`, `SNESSetType()` locates the constructor in that list and calls it to create the specific object. .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESType`, `SNESCreate()`, `SNESDestroy()`, `SNESGetType()`, `SNESSetFromOptions()` @*/ PetscErrorCode SNESSetType(SNES snes, SNESType type) { PetscBool match; PetscErrorCode (*r)(SNES); PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(type, 2); PetscCall(PetscObjectTypeCompare((PetscObject)snes, type, &match)); if (match) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(PetscFunctionListFind(SNESList, type, &r)); PetscCheck(r, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested SNES type %s", type); /* Destroy the previous private SNES context */ PetscTryTypeMethod(snes, destroy); /* Reinitialize type-specific function pointers in SNESOps structure */ snes->ops->reset = NULL; snes->ops->setup = NULL; snes->ops->solve = NULL; snes->ops->view = NULL; snes->ops->setfromoptions = NULL; snes->ops->destroy = NULL; /* It may happen the user has customized the line search before calling SNESSetType */ if (((PetscObject)snes)->type_name) PetscCall(SNESLineSearchDestroy(&snes->linesearch)); /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ snes->setupcalled = PETSC_FALSE; PetscCall(PetscObjectChangeTypeName((PetscObject)snes, type)); PetscCall((*r)(snes)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetType - Gets the `SNES` method type and name (as a string). Not Collective Input Parameter: . snes - nonlinear solver context Output Parameter: . type - `SNES` method (a character string) Level: intermediate .seealso: [](ch_snes), `SNESSetType()`, `SNESType`, `SNESSetFromOptions()`, `SNES` @*/ PetscErrorCode SNESGetType(SNES snes, SNESType *type) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(type, 2); *type = ((PetscObject)snes)->type_name; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetSolution - Sets the solution vector for use by the `SNES` routines. Logically Collective Input Parameters: + snes - the `SNES` context obtained from `SNESCreate()` - u - the solution vector Level: beginner .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESGetSolution()`, `Vec` @*/ PetscErrorCode SNESSetSolution(SNES snes, Vec u) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(u, VEC_CLASSID, 2); PetscCall(PetscObjectReference((PetscObject)u)); PetscCall(VecDestroy(&snes->vec_sol)); snes->vec_sol = u; PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMShellSetGlobalVector(dm, u)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetSolution - Returns the vector where the approximate solution is stored. This is the fine grid solution when using `SNESSetGridSequence()`. Not Collective, but `x` is parallel if `snes` is parallel Input Parameter: . snes - the `SNES` context Output Parameter: . x - the solution Level: intermediate .seealso: [](ch_snes), `SNESSetSolution()`, `SNESSolve()`, `SNES`, `SNESGetSolutionUpdate()`, `SNESGetFunction()` @*/ PetscErrorCode SNESGetSolution(SNES snes, Vec *x) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(x, 2); *x = snes->vec_sol; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetSolutionUpdate - Returns the vector where the solution update is stored. Not Collective, but `x` is parallel if `snes` is parallel Input Parameter: . snes - the `SNES` context Output Parameter: . x - the solution update Level: advanced .seealso: [](ch_snes), `SNES`, `SNESGetSolution()`, `SNESGetFunction()` @*/ PetscErrorCode SNESGetSolutionUpdate(SNES snes, Vec *x) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(x, 2); *x = snes->vec_sol_update; PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESGetFunction - Returns the function that defines the nonlinear system set with `SNESSetFunction()` Not Collective, but `r` is parallel if `snes` is parallel. Collective if `r` is requested, but has not been created yet. Input Parameter: . snes - the `SNES` context Output Parameters: + r - the vector that is used to store residuals (or `NULL` if you don't want it) . f - the function (or `NULL` if you don't want it); for calling sequence see `SNESFunctionFn` - ctx - the function context (or `NULL` if you don't want it) Level: advanced Note: The vector `r` DOES NOT, in general, contain the current value of the `SNES` nonlinear function .seealso: [](ch_snes), `SNES`, `SNESSolve()`, `SNESSetFunction()`, `SNESGetSolution()`, `SNESFunctionFn` @*/ PetscErrorCode SNESGetFunction(SNES snes, Vec *r, SNESFunctionFn **f, void **ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (r) { if (!snes->vec_func) { if (snes->vec_rhs) { PetscCall(VecDuplicate(snes->vec_rhs, &snes->vec_func)); } else if (snes->vec_sol) { PetscCall(VecDuplicate(snes->vec_sol, &snes->vec_func)); } else if (snes->dm) { PetscCall(DMCreateGlobalVector(snes->dm, &snes->vec_func)); } } *r = snes->vec_func; } PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESGetFunction(dm, f, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESGetNGS - Returns the function and context set with `SNESSetNGS()` Input Parameter: . snes - the `SNES` context Output Parameters: + f - the function (or `NULL`) see `SNESNGSFn` for calling sequence - ctx - the function context (or `NULL`) Level: advanced .seealso: [](ch_snes), `SNESSetNGS()`, `SNESGetFunction()`, `SNESNGSFn` @*/ PetscErrorCode SNESGetNGS(SNES snes, SNESNGSFn **f, void **ctx) { DM dm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(SNESGetDM(snes, &dm)); PetscCall(DMSNESGetNGS(dm, f, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetOptionsPrefix - Sets the prefix used for searching for all `SNES` options in the database. Logically Collective Input Parameters: + snes - the `SNES` context - prefix - the prefix to prepend to all option names Level: advanced Note: A hyphen (-) must NOT be given at the beginning of the prefix name. The first character of all runtime options is AUTOMATICALLY the hyphen. .seealso: [](ch_snes), `SNES`, `SNESSetFromOptions()`, `SNESAppendOptionsPrefix()` @*/ PetscErrorCode SNESSetOptionsPrefix(SNES snes, const char prefix[]) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(PetscObjectSetOptionsPrefix((PetscObject)snes, prefix)); if (!snes->ksp) PetscCall(SNESGetKSP(snes, &snes->ksp)); if (snes->linesearch) { PetscCall(SNESGetLineSearch(snes, &snes->linesearch)); PetscCall(PetscObjectSetOptionsPrefix((PetscObject)snes->linesearch, prefix)); } PetscCall(KSPSetOptionsPrefix(snes->ksp, prefix)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESAppendOptionsPrefix - Appends to the prefix used for searching for all `SNES` options in the database. Logically Collective Input Parameters: + snes - the `SNES` context - prefix - the prefix to prepend to all option names Level: advanced Note: A hyphen (-) must NOT be given at the beginning of the prefix name. The first character of all runtime options is AUTOMATICALLY the hyphen. .seealso: [](ch_snes), `SNESGetOptionsPrefix()`, `SNESSetOptionsPrefix()` @*/ PetscErrorCode SNESAppendOptionsPrefix(SNES snes, const char prefix[]) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)snes, prefix)); if (!snes->ksp) PetscCall(SNESGetKSP(snes, &snes->ksp)); if (snes->linesearch) { PetscCall(SNESGetLineSearch(snes, &snes->linesearch)); PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)snes->linesearch, prefix)); } PetscCall(KSPAppendOptionsPrefix(snes->ksp, prefix)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetOptionsPrefix - Gets the prefix used for searching for all `SNES` options in the database. Not Collective Input Parameter: . snes - the `SNES` context Output Parameter: . prefix - pointer to the prefix string used Level: advanced .seealso: [](ch_snes), `SNES`, `SNESSetOptionsPrefix()`, `SNESAppendOptionsPrefix()` @*/ PetscErrorCode SNESGetOptionsPrefix(SNES snes, const char *prefix[]) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscCall(PetscObjectGetOptionsPrefix((PetscObject)snes, prefix)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C SNESRegister - Adds a method to the nonlinear solver package. Not Collective Input Parameters: + sname - name of a new user-defined solver - function - routine to create method context Level: advanced Note: `SNESRegister()` may be called multiple times to add several user-defined solvers. Example Usage: .vb SNESRegister("my_solver", MySolverCreate); .ve Then, your solver can be chosen with the procedural interface via .vb SNESSetType(snes, "my_solver") .ve or at runtime via the option .vb -snes_type my_solver .ve .seealso: [](ch_snes), `SNESRegisterAll()`, `SNESRegisterDestroy()` @*/ PetscErrorCode SNESRegister(const char sname[], PetscErrorCode (*function)(SNES)) { PetscFunctionBegin; PetscCall(SNESInitializePackage()); PetscCall(PetscFunctionListAdd(&SNESList, sname, function)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode SNESTestLocalMin(SNES snes) { PetscInt N, i, j; Vec u, uh, fh; PetscScalar value; PetscReal norm; PetscFunctionBegin; PetscCall(SNESGetSolution(snes, &u)); PetscCall(VecDuplicate(u, &uh)); PetscCall(VecDuplicate(u, &fh)); /* currently only works for sequential */ PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), "Testing FormFunction() for local min\n")); PetscCall(VecGetSize(u, &N)); for (i = 0; i < N; i++) { PetscCall(VecCopy(u, uh)); PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), "i = %" PetscInt_FMT "\n", i)); for (j = -10; j < 11; j++) { value = PetscSign(j) * PetscExpReal(PetscAbs(j) - 10.0); PetscCall(VecSetValue(uh, i, value, ADD_VALUES)); PetscCall(SNESComputeFunction(snes, uh, fh)); PetscCall(VecNorm(fh, NORM_2, &norm)); PetscCall(PetscPrintf(PetscObjectComm((PetscObject)snes), " j norm %" PetscInt_FMT " %18.16e\n", j, (double)norm)); value = -value; PetscCall(VecSetValue(uh, i, value, ADD_VALUES)); } } PetscCall(VecDestroy(&uh)); PetscCall(VecDestroy(&fh)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESKSPSetUseEW - Sets `SNES` to the use Eisenstat-Walker method for computing relative tolerance for linear solvers within an inexact Newton method. Logically Collective Input Parameters: + snes - `SNES` context - flag - `PETSC_TRUE` or `PETSC_FALSE` Options Database Keys: + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence . -snes_ksp_ew_version ver - version of Eisenstat-Walker method . -snes_ksp_ew_rtol0 - Sets rtol0 . -snes_ksp_ew_rtolmax - Sets rtolmax . -snes_ksp_ew_gamma - Sets gamma . -snes_ksp_ew_alpha - Sets alpha . -snes_ksp_ew_alpha2 - Sets alpha2 - -snes_ksp_ew_threshold - Sets threshold Level: advanced Note: The default is to use a constant relative tolerance for the inner linear solvers. Alternatively, one can use the Eisenstat-Walker method {cite}`ew96`, where the relative convergence tolerance is reset at each Newton iteration according progress of the nonlinear solver. .seealso: [](ch_snes), `KSP`, `SNES`, `SNESKSPGetUseEW()`, `SNESKSPGetParametersEW()`, `SNESKSPSetParametersEW()` @*/ PetscErrorCode SNESKSPSetUseEW(SNES snes, PetscBool flag) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveBool(snes, flag, 2); snes->ksp_ewconv = flag; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESKSPGetUseEW - Gets if `SNES` is using Eisenstat-Walker method for computing relative tolerance for linear solvers within an inexact Newton method. Not Collective Input Parameter: . snes - `SNES` context Output Parameter: . flag - `PETSC_TRUE` or `PETSC_FALSE` Level: advanced .seealso: [](ch_snes), `SNESKSPSetUseEW()`, `SNESKSPGetParametersEW()`, `SNESKSPSetParametersEW()` @*/ PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(flag, 2); *flag = snes->ksp_ewconv; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker convergence criteria for the linear solvers within an inexact Newton method. Logically Collective Input Parameters: + snes - `SNES` context . version - version 1, 2 (default is 2), 3 or 4 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) . alpha - power for version 2 rtol computation (1 < alpha <= 2) . alpha2 - power for safeguard - threshold - threshold for imposing safeguard (0 < threshold < 1) Level: advanced Notes: Version 3 was contributed by Luis Chacon, June 2006. Use `PETSC_CURRENT` to retain the default for any of the parameters. .seealso: [](ch_snes), `SNES`, `SNESKSPSetUseEW()`, `SNESKSPGetUseEW()`, `SNESKSPGetParametersEW()` @*/ PetscErrorCode SNESKSPSetParametersEW(SNES snes, PetscInt version, PetscReal rtol_0, PetscReal rtol_max, PetscReal gamma, PetscReal alpha, PetscReal alpha2, PetscReal threshold) { SNESKSPEW *kctx; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); kctx = (SNESKSPEW *)snes->kspconvctx; PetscCheck(kctx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "No Eisenstat-Walker context existing"); PetscValidLogicalCollectiveInt(snes, version, 2); PetscValidLogicalCollectiveReal(snes, rtol_0, 3); PetscValidLogicalCollectiveReal(snes, rtol_max, 4); PetscValidLogicalCollectiveReal(snes, gamma, 5); PetscValidLogicalCollectiveReal(snes, alpha, 6); PetscValidLogicalCollectiveReal(snes, alpha2, 7); PetscValidLogicalCollectiveReal(snes, threshold, 8); if (version != PETSC_CURRENT) kctx->version = version; if (rtol_0 != (PetscReal)PETSC_CURRENT) kctx->rtol_0 = rtol_0; if (rtol_max != (PetscReal)PETSC_CURRENT) kctx->rtol_max = rtol_max; if (gamma != (PetscReal)PETSC_CURRENT) kctx->gamma = gamma; if (alpha != (PetscReal)PETSC_CURRENT) kctx->alpha = alpha; if (alpha2 != (PetscReal)PETSC_CURRENT) kctx->alpha2 = alpha2; if (threshold != (PetscReal)PETSC_CURRENT) kctx->threshold = threshold; PetscCheck(kctx->version >= 1 && kctx->version <= 4, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only versions 1 to 4 are supported: %" PetscInt_FMT, kctx->version); PetscCheck(kctx->rtol_0 >= 0.0 && kctx->rtol_0 < 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "0.0 <= rtol_0 < 1.0: %g", (double)kctx->rtol_0); PetscCheck(kctx->rtol_max >= 0.0 && kctx->rtol_max < 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "0.0 <= rtol_max (%g) < 1.0", (double)kctx->rtol_max); PetscCheck(kctx->gamma >= 0.0 && kctx->gamma <= 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "0.0 <= gamma (%g) <= 1.0", (double)kctx->gamma); PetscCheck(kctx->alpha > 1.0 && kctx->alpha <= 2.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "1.0 < alpha (%g) <= 2.0", (double)kctx->alpha); PetscCheck(kctx->threshold > 0.0 && kctx->threshold < 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "0.0 < threshold (%g) < 1.0", (double)kctx->threshold); PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker convergence criteria for the linear solvers within an inexact Newton method. Not Collective Input Parameter: . snes - `SNES` context Output Parameters: + version - version 1, 2 (default is 2), 3 or 4 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) . gamma - multiplicative factor for version 2 rtol computation (0 <= gamma2 <= 1) . alpha - power for version 2 rtol computation (1 < alpha <= 2) . alpha2 - power for safeguard - threshold - threshold for imposing safeguard (0 < threshold < 1) Level: advanced .seealso: [](ch_snes), `SNES`, `SNESKSPSetUseEW()`, `SNESKSPGetUseEW()`, `SNESKSPSetParametersEW()` @*/ PetscErrorCode SNESKSPGetParametersEW(SNES snes, PetscInt *version, PetscReal *rtol_0, PetscReal *rtol_max, PetscReal *gamma, PetscReal *alpha, PetscReal *alpha2, PetscReal *threshold) { SNESKSPEW *kctx; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); kctx = (SNESKSPEW *)snes->kspconvctx; PetscCheck(kctx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "No Eisenstat-Walker context existing"); if (version) *version = kctx->version; if (rtol_0) *rtol_0 = kctx->rtol_0; if (rtol_max) *rtol_max = kctx->rtol_max; if (gamma) *gamma = kctx->gamma; if (alpha) *alpha = kctx->alpha; if (alpha2) *alpha2 = kctx->alpha2; if (threshold) *threshold = kctx->threshold; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode KSPPreSolve_SNESEW(KSP ksp, Vec b, Vec x, void *ctx) { SNES snes = (SNES)ctx; SNESKSPEW *kctx = (SNESKSPEW *)snes->kspconvctx; PetscReal rtol = PETSC_CURRENT, stol; PetscFunctionBegin; if (!snes->ksp_ewconv) PetscFunctionReturn(PETSC_SUCCESS); if (!snes->iter) { rtol = kctx->rtol_0; /* first time in, so use the original user rtol */ PetscCall(VecNorm(snes->vec_func, NORM_2, &kctx->norm_first)); } else { PetscCheck(kctx->version >= 1 && kctx->version <= 4, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only versions 1-4 are supported: %" PetscInt_FMT, kctx->version); if (kctx->version == 1) { rtol = PetscAbsReal(snes->norm - kctx->lresid_last) / kctx->norm_last; stol = PetscPowReal(kctx->rtol_last, kctx->alpha2); if (stol > kctx->threshold) rtol = PetscMax(rtol, stol); } else if (kctx->version == 2) { rtol = kctx->gamma * PetscPowReal(snes->norm / kctx->norm_last, kctx->alpha); stol = kctx->gamma * PetscPowReal(kctx->rtol_last, kctx->alpha); if (stol > kctx->threshold) rtol = PetscMax(rtol, stol); } else if (kctx->version == 3) { /* contributed by Luis Chacon, June 2006. */ rtol = kctx->gamma * PetscPowReal(snes->norm / kctx->norm_last, kctx->alpha); /* safeguard: avoid sharp decrease of rtol */ stol = kctx->gamma * PetscPowReal(kctx->rtol_last, kctx->alpha); stol = PetscMax(rtol, stol); rtol = PetscMin(kctx->rtol_0, stol); /* safeguard: avoid oversolving */ stol = kctx->gamma * (kctx->norm_first * snes->rtol) / snes->norm; stol = PetscMax(rtol, stol); rtol = PetscMin(kctx->rtol_0, stol); } else /* if (kctx->version == 4) */ { /* H.-B. An et al. Journal of Computational and Applied Mathematics 200 (2007) 47-60 */ PetscReal ared = PetscAbsReal(kctx->norm_last - snes->norm); PetscReal pred = PetscAbsReal(kctx->norm_last - kctx->lresid_last); PetscReal rk = ared / pred; if (rk < kctx->v4_p1) rtol = 1. - 2. * kctx->v4_p1; else if (rk < kctx->v4_p2) rtol = kctx->rtol_last; else if (rk < kctx->v4_p3) rtol = kctx->v4_m1 * kctx->rtol_last; else rtol = kctx->v4_m2 * kctx->rtol_last; if (kctx->rtol_last_2 > kctx->v4_m3 && kctx->rtol_last > kctx->v4_m3 && kctx->rk_last_2 < kctx->v4_p1 && kctx->rk_last < kctx->v4_p1) rtol = kctx->v4_m4 * kctx->rtol_last; kctx->rtol_last_2 = kctx->rtol_last; kctx->rk_last_2 = kctx->rk_last; kctx->rk_last = rk; } } /* safeguard: avoid rtol greater than rtol_max */ rtol = PetscMin(rtol, kctx->rtol_max); PetscCall(KSPSetTolerances(ksp, rtol, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT)); PetscCall(PetscInfo(snes, "iter %" PetscInt_FMT ", Eisenstat-Walker (version %" PetscInt_FMT ") KSP rtol=%g\n", snes->iter, kctx->version, (double)rtol)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode KSPPostSolve_SNESEW(KSP ksp, Vec b, Vec x, void *ctx) { SNES snes = (SNES)ctx; SNESKSPEW *kctx = (SNESKSPEW *)snes->kspconvctx; PCSide pcside; Vec lres; PetscFunctionBegin; if (!snes->ksp_ewconv) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(KSPGetTolerances(ksp, &kctx->rtol_last, NULL, NULL, NULL)); kctx->norm_last = snes->norm; if (kctx->version == 1 || kctx->version == 4) { PC pc; PetscBool getRes; PetscCall(KSPGetPC(ksp, &pc)); PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCNONE, &getRes)); if (!getRes) { KSPNormType normtype; PetscCall(KSPGetNormType(ksp, &normtype)); getRes = (PetscBool)(normtype == KSP_NORM_UNPRECONDITIONED); } PetscCall(KSPGetPCSide(ksp, &pcside)); if (pcside == PC_RIGHT || getRes) { /* KSP residual is true linear residual */ PetscCall(KSPGetResidualNorm(ksp, &kctx->lresid_last)); } else { /* KSP residual is preconditioned residual */ /* compute true linear residual norm */ Mat J; PetscCall(KSPGetOperators(ksp, &J, NULL)); PetscCall(VecDuplicate(b, &lres)); PetscCall(MatMult(J, x, lres)); PetscCall(VecAYPX(lres, -1.0, b)); PetscCall(VecNorm(lres, NORM_2, &kctx->lresid_last)); PetscCall(VecDestroy(&lres)); } } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetKSP - Returns the `KSP` context for a `SNES` solver. Not Collective, but if `snes` is parallel, then `ksp` is parallel Input Parameter: . snes - the `SNES` context Output Parameter: . ksp - the `KSP` context Level: beginner Notes: The user can then directly manipulate the `KSP` context to set various options, etc. Likewise, the user can then extract and manipulate the `PC` contexts as well. Some `SNESType`s do not use a `KSP` but a `KSP` is still returned by this function, changes to that `KSP` will have no effect. .seealso: [](ch_snes), `SNES`, `KSP`, `PC`, `KSPGetPC()`, `SNESCreate()`, `KSPCreate()`, `SNESSetKSP()` @*/ PetscErrorCode SNESGetKSP(SNES snes, KSP *ksp) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(ksp, 2); if (!snes->ksp) { PetscCall(KSPCreate(PetscObjectComm((PetscObject)snes), &snes->ksp)); PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes->ksp, (PetscObject)snes, 1)); PetscCall(KSPSetPreSolve(snes->ksp, KSPPreSolve_SNESEW, snes)); PetscCall(KSPSetPostSolve(snes->ksp, KSPPostSolve_SNESEW, snes)); PetscCall(KSPMonitorSetFromOptions(snes->ksp, "-snes_monitor_ksp", "snes_preconditioned_residual", snes)); PetscCall(PetscObjectSetOptions((PetscObject)snes->ksp, ((PetscObject)snes)->options)); } *ksp = snes->ksp; PetscFunctionReturn(PETSC_SUCCESS); } #include /*@ SNESSetDM - Sets the `DM` that may be used by some `SNES` nonlinear solvers or their underlying preconditioners Logically Collective Input Parameters: + snes - the nonlinear solver context - dm - the `DM`, cannot be `NULL` Level: intermediate Note: A `DM` can only be used for solving one problem at a time because information about the problem is stored on the `DM`, even when not using interfaces like `DMSNESSetFunction()`. Use `DMClone()` to get a distinct `DM` when solving different problems using the same function space. .seealso: [](ch_snes), `DM`, `SNES`, `SNESGetDM()`, `KSPSetDM()`, `KSPGetDM()` @*/ PetscErrorCode SNESSetDM(SNES snes, DM dm) { KSP ksp; DMSNES sdm; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(dm, DM_CLASSID, 2); PetscCall(PetscObjectReference((PetscObject)dm)); if (snes->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */ if (snes->dm->dmsnes && !dm->dmsnes) { PetscCall(DMCopyDMSNES(snes->dm, dm)); PetscCall(DMGetDMSNES(snes->dm, &sdm)); if (sdm->originaldm == snes->dm) sdm->originaldm = dm; /* Grant write privileges to the replacement DM */ } PetscCall(DMCoarsenHookRemove(snes->dm, DMCoarsenHook_SNESVecSol, DMRestrictHook_SNESVecSol, snes)); PetscCall(DMDestroy(&snes->dm)); } snes->dm = dm; snes->dmAuto = PETSC_FALSE; PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(KSPSetDM(ksp, dm)); PetscCall(KSPSetDMActive(ksp, PETSC_FALSE)); if (snes->npc) { PetscCall(SNESSetDM(snes->npc, snes->dm)); PetscCall(SNESSetNPCSide(snes, snes->npcside)); } PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetDM - Gets the `DM` that may be used by some `SNES` nonlinear solvers/preconditioners Not Collective but `dm` obtained is parallel on `snes` Input Parameter: . snes - the `SNES` context Output Parameter: . dm - the `DM` Level: intermediate .seealso: [](ch_snes), `DM`, `SNES`, `SNESSetDM()`, `KSPSetDM()`, `KSPGetDM()` @*/ PetscErrorCode SNESGetDM(SNES snes, DM *dm) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); if (!snes->dm) { PetscCall(DMShellCreate(PetscObjectComm((PetscObject)snes), &snes->dm)); snes->dmAuto = PETSC_TRUE; } *dm = snes->dm; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetNPC - Sets the nonlinear preconditioner to be used. Collective Input Parameters: + snes - iterative context obtained from `SNESCreate()` - npc - the `SNES` nonlinear preconditioner object Options Database Key: . -npc_snes_type - set the type of the `SNES` to use as the nonlinear preconditioner Level: developer Notes: This is rarely used, rather use `SNESGetNPC()` to retrieve the preconditioner and configure it using the API. Only some `SNESType` can use a nonlinear preconditioner .seealso: [](ch_snes), `SNES`, `SNESNGS`, `SNESFAS`, `SNESGetNPC()`, `SNESHasNPC()` @*/ PetscErrorCode SNESSetNPC(SNES snes, SNES npc) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(npc, SNES_CLASSID, 2); PetscCheckSameComm(snes, 1, npc, 2); PetscCall(PetscObjectReference((PetscObject)npc)); PetscCall(SNESDestroy(&snes->npc)); snes->npc = npc; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetNPC - Gets a nonlinear preconditioning solver SNES` to be used to precondition the original nonlinear solver. Not Collective; but any changes to the obtained the `pc` object must be applied collectively Input Parameter: . snes - iterative context obtained from `SNESCreate()` Output Parameter: . pc - the `SNES` preconditioner context Options Database Key: . -npc_snes_type - set the type of the `SNES` to use as the nonlinear preconditioner Level: advanced Notes: If a `SNES` was previously set with `SNESSetNPC()` then that value is returned, otherwise a new `SNES` object is created that will be used as the nonlinear preconditioner for the current `SNES`. The (preconditioner) `SNES` returned automatically inherits the same nonlinear function and Jacobian supplied to the original `SNES`. These may be overwritten if needed. Use the options database prefixes `-npc_snes`, `-npc_ksp`, etc., to control the configuration of the nonlinear preconditioner .seealso: [](ch_snes), `SNESSetNPC()`, `SNESHasNPC()`, `SNES`, `SNESCreate()` @*/ PetscErrorCode SNESGetNPC(SNES snes, SNES *pc) { const char *optionsprefix; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(pc, 2); if (!snes->npc) { void *ctx; PetscCall(SNESCreate(PetscObjectComm((PetscObject)snes), &snes->npc)); PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes->npc, (PetscObject)snes, 1)); PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix)); PetscCall(SNESSetOptionsPrefix(snes->npc, optionsprefix)); PetscCall(SNESAppendOptionsPrefix(snes->npc, "npc_")); if (snes->ops->usercompute) { PetscCall(SNESSetComputeApplicationContext(snes, snes->ops->usercompute, snes->ops->ctxdestroy)); } else { PetscCall(SNESGetApplicationContext(snes, &ctx)); PetscCall(SNESSetApplicationContext(snes->npc, ctx)); } PetscCall(SNESSetCountersReset(snes->npc, PETSC_FALSE)); } *pc = snes->npc; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESHasNPC - Returns whether a nonlinear preconditioner is associated with the given `SNES` Not Collective Input Parameter: . snes - iterative context obtained from `SNESCreate()` Output Parameter: . has_npc - whether the `SNES` has a nonlinear preconditioner or not Level: developer .seealso: [](ch_snes), `SNESSetNPC()`, `SNESGetNPC()` @*/ PetscErrorCode SNESHasNPC(SNES snes, PetscBool *has_npc) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(has_npc, 2); *has_npc = snes->npc ? PETSC_TRUE : PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetNPCSide - Sets the nonlinear preconditioning side used by the nonlinear preconditioner inside `SNES`. Logically Collective Input Parameter: . snes - iterative context obtained from `SNESCreate()` Output Parameter: . side - the preconditioning side, where side is one of .vb PC_LEFT - left preconditioning PC_RIGHT - right preconditioning (default for most nonlinear solvers) .ve Options Database Key: . -snes_npc_side - nonlinear preconditioner side Level: intermediate Note: `SNESNRICHARDSON` and `SNESNCG` only support left preconditioning. .seealso: [](ch_snes), `SNES`, `SNESGetNPC()`, `SNESNRICHARDSON`, `SNESNCG`, `SNESType`, `SNESGetNPCSide()`, `KSPSetPCSide()`, `PC_LEFT`, `PC_RIGHT`, `PCSide` @*/ PetscErrorCode SNESSetNPCSide(SNES snes, PCSide side) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidLogicalCollectiveEnum(snes, side, 2); if (side == PC_SIDE_DEFAULT) side = PC_RIGHT; PetscCheck((side == PC_LEFT) || (side == PC_RIGHT), PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONG, "Only PC_LEFT and PC_RIGHT are supported"); snes->npcside = side; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetNPCSide - Gets the preconditioning side used by the nonlinear preconditioner inside `SNES`. Not Collective Input Parameter: . snes - iterative context obtained from `SNESCreate()` Output Parameter: . side - the preconditioning side, where side is one of .vb `PC_LEFT` - left preconditioning `PC_RIGHT` - right preconditioning (default for most nonlinear solvers) .ve Level: intermediate .seealso: [](ch_snes), `SNES`, `SNESGetNPC()`, `SNESSetNPCSide()`, `KSPGetPCSide()`, `PC_LEFT`, `PC_RIGHT`, `PCSide` @*/ PetscErrorCode SNESGetNPCSide(SNES snes, PCSide *side) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(side, 2); *side = snes->npcside; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESSetLineSearch - Sets the `SNESLineSearch` to be used for a given `SNES` Collective Input Parameters: + snes - iterative context obtained from `SNESCreate()` - linesearch - the linesearch object Level: developer Note: This is almost never used, rather one uses `SNESGetLineSearch()` to retrieve the line search and set options on it to configure it using the API). .seealso: [](ch_snes), `SNES`, `SNESLineSearch`, `SNESGetLineSearch()` @*/ PetscErrorCode SNESSetLineSearch(SNES snes, SNESLineSearch linesearch) { PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscValidHeaderSpecific(linesearch, SNESLINESEARCH_CLASSID, 2); PetscCheckSameComm(snes, 1, linesearch, 2); PetscCall(PetscObjectReference((PetscObject)linesearch)); PetscCall(SNESLineSearchDestroy(&snes->linesearch)); snes->linesearch = linesearch; PetscFunctionReturn(PETSC_SUCCESS); } /*@ SNESGetLineSearch - Returns the line search associated with the `SNES`. Not Collective Input Parameter: . snes - iterative context obtained from `SNESCreate()` Output Parameter: . linesearch - linesearch context Level: beginner Notes: It creates a default line search instance which can be configured as needed in case it has not been already set with `SNESSetLineSearch()`. You can also use the options database keys `-snes_linesearch_*` to configure the line search. See `SNESLineSearchSetFromOptions()` for the possible options. .seealso: [](ch_snes), `SNESLineSearch`, `SNESSetLineSearch()`, `SNESLineSearchCreate()`, `SNESLineSearchSetFromOptions()` @*/ PetscErrorCode SNESGetLineSearch(SNES snes, SNESLineSearch *linesearch) { const char *optionsprefix; PetscFunctionBegin; PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); PetscAssertPointer(linesearch, 2); if (!snes->linesearch) { PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix)); PetscCall(SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &snes->linesearch)); PetscCall(SNESLineSearchSetSNES(snes->linesearch, snes)); PetscCall(SNESLineSearchAppendOptionsPrefix(snes->linesearch, optionsprefix)); PetscCall(PetscObjectIncrementTabLevel((PetscObject)snes->linesearch, (PetscObject)snes, 1)); } *linesearch = snes->linesearch; PetscFunctionReturn(PETSC_SUCCESS); }