1 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdmshell.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "SNESApplyNPC" 7 /*@ 8 SNESApplyNPC - Calls SNESSolve() on preconditioner for the SNES 9 10 Collective on SNES 11 12 Input Parameters: 13 + snes - the SNES context 14 . x - input vector 15 - f - optional; the function evaluation on x 16 17 Output Parameter: 18 . y - function vector, as set by SNESSetFunction() 19 20 Notes: 21 SNESComputeFunction() should be called on x before SNESApplyNPC() is called, as it is 22 with SNESComuteJacobian(). 23 24 Level: developer 25 26 .keywords: SNES, nonlinear, compute, function 27 28 .seealso: SNESGetNPC(),SNESSetNPC(),SNESComputeFunction() 29 @*/ 30 PetscErrorCode SNESApplyNPC(SNES snes,Vec x,Vec f,Vec y) 31 { 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 36 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 37 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 38 PetscCheckSameComm(snes,1,x,2); 39 PetscCheckSameComm(snes,1,y,3); 40 ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr); 41 if (snes->pc) { 42 if (f) { 43 ierr = SNESSetInitialFunction(snes->pc,f);CHKERRQ(ierr); 44 } 45 ierr = VecCopy(x,y);CHKERRQ(ierr); 46 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr); 47 ierr = SNESSolve(snes->pc,snes->vec_rhs,y);CHKERRQ(ierr); 48 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,x,y,0);CHKERRQ(ierr); 49 ierr = VecAYPX(y,-1.0,x);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "SNESComputeFunctionDefaultNPC" 57 PetscErrorCode SNESComputeFunctionDefaultNPC(SNES snes,Vec X,Vec F) 58 { 59 /* This is to be used as an argument to SNESMF -- NOT as a "function" */ 60 SNESConvergedReason reason; 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 if (snes->pc) { 65 ierr = SNESApplyNPC(snes,X,NULL,F);CHKERRQ(ierr); 66 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 67 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 68 ierr = SNESSetFunctionDomainError(snes);CHKERRQ(ierr); 69 } 70 } else { 71 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 72 } 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "SNESGetNPCFunction" 78 /*@ 79 SNESGetNPCFunction - Gets the function from a preconditioner after SNESSolve() has been called. 80 81 Collective on SNES 82 83 Input Parameters: 84 . snes - the SNES context 85 86 Output Parameter: 87 . F - function vector 88 . fnorm - the norm of F 89 90 Level: developer 91 92 .keywords: SNES, nonlinear, function 93 94 .seealso: SNESGetNPC(),SNESSetNPC(),SNESComputeFunction(),SNESApplyNPC(),SNESSolve() 95 @*/ 96 PetscErrorCode SNESGetNPCFunction(SNES snes,Vec F,PetscReal *fnorm) 97 { 98 PetscErrorCode ierr; 99 PCSide npcside; 100 SNESFunctionType functype; 101 SNESNormSchedule normschedule; 102 Vec FPC,XPC; 103 104 PetscFunctionBegin; 105 if (snes->pc) { 106 ierr = SNESGetNPCSide(snes->pc,&npcside);CHKERRQ(ierr); 107 ierr = SNESGetFunctionType(snes->pc,&functype);CHKERRQ(ierr); 108 ierr = SNESGetNormSchedule(snes->pc,&normschedule);CHKERRQ(ierr); 109 110 /* check if the function is valid based upon how the inner solver is preconditioned */ 111 if (normschedule != SNES_NORM_NONE && normschedule != SNES_NORM_INITIAL_ONLY && (npcside == PC_RIGHT || functype == SNES_FUNCTION_UNPRECONDITIONED)) { 112 ierr = SNESGetFunction(snes->pc,&FPC,NULL,NULL);CHKERRQ(ierr); 113 if (FPC) { 114 if (fnorm) {ierr = VecNorm(FPC,NORM_2,fnorm);CHKERRQ(ierr);} 115 ierr = VecCopy(FPC,F);CHKERRQ(ierr); 116 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Preconditioner has no function"); 117 } else { 118 ierr = SNESGetSolution(snes->pc,&XPC);CHKERRQ(ierr); 119 if (XPC) { 120 ierr = SNESComputeFunction(snes->pc,XPC,F);CHKERRQ(ierr); 121 if (fnorm) {ierr = VecNorm(F,NORM_2,fnorm);CHKERRQ(ierr);} 122 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Preconditioner has no solution"); 123 } 124 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No preconditioner set"); 125 PetscFunctionReturn(0); 126 } 127