Actual source code: ex7.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\
3: matrix-free techniques with user-provided explicit preconditioner matrix.\n\n";
5: #include <petscsnes.h>
7: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
8: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
9: extern PetscErrorCode OtherFunctionForDifferencing(void*,Vec,Vec);
10: extern PetscErrorCode FormInitialGuess(SNES,Vec);
11: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
13: typedef struct {
14: PetscViewer viewer;
15: } MonitorCtx;
17: typedef struct {
18: PetscBool variant;
19: } AppCtx;
21: int main(int argc,char **argv)
22: {
23: SNES snes; /* SNES context */
24: SNESType type = SNESNEWTONLS; /* default nonlinear solution method */
25: Vec x,r,F,U; /* vectors */
26: Mat J,B; /* Jacobian matrix-free, explicit preconditioner */
27: MonitorCtx monP; /* monitoring context */
28: AppCtx user; /* user-defined work context */
29: PetscScalar h,xp = 0.0,v;
30: PetscInt its,n = 5,i;
33: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
34: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
35: PetscOptionsHasName(NULL,NULL,"-variant",&user.variant);
36: h = 1.0/(n-1);
38: /* Set up data structures */
39: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&monP.viewer);
40: VecCreateSeq(PETSC_COMM_SELF,n,&x);
41: PetscObjectSetName((PetscObject)x,"Approximate Solution");
42: VecDuplicate(x,&r);
43: VecDuplicate(x,&F);
44: VecDuplicate(x,&U);
45: PetscObjectSetName((PetscObject)U,"Exact Solution");
47: /* create explict matrix preconditioner */
48: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&B);
50: /* Store right-hand-side of PDE and exact solution */
51: for (i=0; i<n; i++) {
52: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
53: VecSetValues(F,1,&i,&v,INSERT_VALUES);
54: v = xp*xp*xp;
55: VecSetValues(U,1,&i,&v,INSERT_VALUES);
56: xp += h;
57: }
59: /* Create nonlinear solver */
60: SNESCreate(PETSC_COMM_WORLD,&snes);
61: SNESSetType(snes,type);
63: /* Set various routines and options */
64: SNESSetFunction(snes,r,FormFunction,F);
65: if (user.variant) {
66: /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
67: MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);
68: MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);
69: } else {
70: /* create matrix free matrix for Jacobian */
71: MatCreateSNESMF(snes,&J);
72: /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */
73: /* note we use the same context for this function as FormFunction, the F vector */
74: MatMFFDSetFunction(J,OtherFunctionForDifferencing,F);
75: }
77: /* Set various routines and options */
78: SNESSetJacobian(snes,J,B,FormJacobian,&user);
79: SNESMonitorSet(snes,Monitor,&monP,0);
80: SNESSetFromOptions(snes);
82: /* Solve nonlinear system */
83: FormInitialGuess(snes,x);
84: SNESSolve(snes,NULL,x);
85: SNESGetIterationNumber(snes,&its);
86: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
88: /* Free data structures */
89: VecDestroy(&x); VecDestroy(&r);
90: VecDestroy(&U); VecDestroy(&F);
91: MatDestroy(&J); MatDestroy(&B);
92: SNESDestroy(&snes);
93: PetscViewerDestroy(&monP.viewer);
94: PetscFinalize();
95: return ierr;
96: }
97: /* -------------------- Evaluate Function F(x) --------------------- */
99: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
100: {
101: const PetscScalar *xx;
102: PetscScalar *ff,*FF,d;
103: PetscInt i,n;
104: PetscErrorCode ierr;
106: VecGetArrayRead(x,&xx);
107: VecGetArray(f,&ff);
108: VecGetArray((Vec) dummy,&FF);
109: VecGetSize(x,&n);
110: d = (PetscReal)(n - 1); d = d*d;
111: ff[0] = xx[0];
112: for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
113: ff[n-1] = xx[n-1] - 1.0;
114: VecRestoreArrayRead(x,&xx);
115: VecRestoreArray(f,&ff);
116: VecRestoreArray((Vec)dummy,&FF);
117: return 0;
118: }
120: /*
122: Example function that when differenced produces the same matrix free Jacobian as FormFunction()
123: this is provided to show how a user can provide a different function
124: */
125: PetscErrorCode OtherFunctionForDifferencing(void *dummy,Vec x,Vec f)
126: {
129: FormFunction(NULL,x,f,dummy);
130: VecShift(f,1.0);
131: return 0;
132: }
134: /* -------------------- Form initial approximation ----------------- */
136: PetscErrorCode FormInitialGuess(SNES snes,Vec x)
137: {
139: PetscScalar pfive = .50;
140: VecSet(x,pfive);
141: return 0;
142: }
143: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
144: /* Evaluates a matrix that is used to precondition the matrix-free
145: jacobian. In this case, the explict preconditioner matrix is
146: also EXACTLY the Jacobian. In general, it would be some lower
147: order, simplified apprioximation */
149: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
150: {
151: const PetscScalar *xx;
152: PetscScalar A[3],d;
153: PetscInt i,n,j[3];
154: PetscErrorCode ierr;
155: AppCtx *user = (AppCtx*) dummy;
157: VecGetArrayRead(x,&xx);
158: VecGetSize(x,&n);
159: d = (PetscReal)(n - 1); d = d*d;
161: i = 0; A[0] = 1.0;
162: MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
163: for (i=1; i<n-1; i++) {
164: j[0] = i - 1; j[1] = i; j[2] = i + 1;
165: A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d;
166: MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
167: }
168: i = n-1; A[0] = 1.0;
169: MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
170: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
171: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
172: VecRestoreArrayRead(x,&xx);
174: if (user->variant) {
175: MatMFFDSetBase(jac,x,NULL);
176: }
177: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
178: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
179: return 0;
180: }
181: /* -------------------- User-defined monitor ----------------------- */
183: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy)
184: {
186: MonitorCtx *monP = (MonitorCtx*) dummy;
187: Vec x;
188: MPI_Comm comm;
190: PetscObjectGetComm((PetscObject)snes,&comm);
191: PetscFPrintf(comm,stdout,"iter = %D, SNES Function norm %g \n",its,(double)fnorm);
192: SNESGetSolution(snes,&x);
193: VecView(x,monP->viewer);
194: return 0;
195: }
198: /*TEST
200: test:
201: args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_cancel -snes_monitor_short
202: requires: x
204: test:
205: suffix: 2
206: args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_cancel -snes_monitor_short
207: output_file: output/ex7_1.out
208: requires: x
210: TEST*/